UFDC Home  myUFDC Home  Help 



Full Text  
PAGE 1 1 FROM ELECTRONIC STRUCTURE OF PO INT DEFECTS TO PHYSICAL PROPERT IE S OF COMPLEX MATERIALS USING ATOMIC LEVEL SIMULATIONS By HAIXUAN XU 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 PAGE 2 2 2010 Haixuan Xu PAGE 3 3 To my family with love and gratitude PAGE 4 4 ACKNOWLEDGMENTS First I would like to express my deepest gratitude and respect to my advisor Prof. Simon R. Phillpot not only for his invaluable gu idance throughout my Ph.D study, but also h is enthusiasm towards science and patience with students with a great sense of humor w hich has made my time at UF extremely special. I admire him as a great researcher, mentor and motivator. I feel extremely fortunate to have him as my advisor I would like to expand my appreciation to my committee members Prof. Susan B. Sinnott, Prof. Fere shteh Ebrahimi Prof. Eric D. Wachsman Prof. Jacob Jones and Prof. Aravind Asthagiri for their continuous support and advice. I would like to thank Dr. Roger E. Stoller, with whom I had the fortune to work as a student intern at Oak Ridge National Laborat ory. His kindness and insights into the subject have mad e my time there enjoyable and unforgettable. I would also like to convey my thanks to my experimental collaborator Prof. Venkatraman Gopalan at PSU and Prof. Volkmar Dierolf at Lehigh University. I w ould like to express my appreciation to all the members of Computational Materials Science Focus Group (CMSFG) especially Dr. Rakesh Behera, Dr. Jun He, Dr Tao Liang, Prof Man Yao (Dalian University of Technology), Dr. Alex Chernatinsky, and Dr. Wen Dung Hsu The CMSFG ha d a positive impact on my research and life outside the school. Finally, I would like to express my gratitude to my wife, Luyi Yin, for her unconditional love and invaluable support. She has been the constant inspiration to me. My heartf elt appreciation goes out to my parents, to whom I owe the most. Their love and belief has shape me to who I am today. This dissertation is dedicated to them. PAGE 5 5 TABLE OF CONTENTS page ACKNOWLEDGMENTS .................................................................................................................... 4 LIST OF TABLES ................................................................................................................................ 8 LIST OF FIGURES ............................................................................................................................ 10 ABSTRACT ........................................................................................................................................ 14 CHAPTER 1 INTRODUCTION ....................................................................................................................... 16 1.1 Point Defects in Crystalline Materials ............................................................................. 16 1.1.1 Energetics and Diffusion of Point Defects .......................................................... 17 1.1.2 Electronic Structure of Point Defect .................................................................... 20 1.1.3 The Effects of Point Defects on Properties of Materials .................................... 23 1.2 Defect Stability in Lithium Niobate ................................................................................. 27 1.3 Mechanical Properties of Ceria based Materials ............................................................ 29 1.4 Radiation Damage in Titanium ........................................................................................ 31 1.5 Overview of This Dissertation .......................................................................................... 32 2 SIMULATION METHODS ....................................................................................................... 34 2.1 Density Functional Theory ............................................................................................... 35 2. 1.1 Schrodinger Equation ........................................................................................... 35 2.1.2 Born Oppenheimer Approximation ..................................................................... 35 2.1.3 Hohenberg Kohn Theorem .................................................................................. 36 2.1.4 Kohn Sham Equation ........................................................................................... 37 2.1.5 Exchange Correlation Functionals ...................................................................... 38 2.1.6 Pseudopotential Approximation .......................................................................... 39 2.1.7 Limitations ............................................................................................................ 40 2.1.8 DFT+U ................................................................................................................. 41 2.2 Molecular Dynamics (MD) Simulation ........................................................................... 41 2.2.1 Interatomic Interactions ........................................................................................ 42 2.2.2 Integration Scheme ............................................................................................... 44 2.2.3 Thermostat ............................................................................................................. 46 2.2.4 Barostat ................................................................................................................. 47 2.2.5 Ensembles .............................................................................................................. 48 3 INTRINSIC DEFE CTS IN LITHIUM NIOBATE ................................................................... 50 3.1 Background of Intrinsic Defects ....................................................................................... 50 3.2 Structure of LiNbO3 and Constitute Binary Oxides ........................................................ 52 3.3 Computational Details ...................................................................................................... 54 PAGE 6 6 3.4 Thermodynamic Framework ............................................................................................ 58 3.5 Defect Formation Energy ................................................................................................. 64 3.5. 1 Formalism .............................................................................................................. 66 3.5.2 Point Defects ......................................................................................................... 67 3.5.3 Defect Clusters ...................................................................................................... 70 3.5.4 Def ect Stability Range .......................................................................................... 73 3.6 Association Effects and Defect Configurations .............................................................. 74 3.7 Diffusion Mechanisms ...................................................................................................... 82 3.8 D iscussion .......................................................................................................................... 86 4 EXTRINSIC DEFECTS IN LITHIUM NIOBATE .................................................................. 91 4.1 Background of Extrinsic Defects ..................................................................................... 91 4.2 Structure of Reference Binary Oxides ............................................................................. 92 4.2.1 Crystal Structure of MgO and FeO ...................................................................... 93 4.2.2 Crystal Structure of Er2O3 Nb2O3 and Fe2O3 .................................................. 94 4.3 Computational Details ...................................................................................................... 96 4.4 Formation Energies and Site Preference ........................................................................ 101 4.4.1 Single Defects ..................................................................................................... 102 4.4.2 Defect Clusters .................................................................................................... 106 4.5 Charge Transfer Levels ................................................................................................... 114 4.6 Discussions ...................................................................................................................... 117 4.6.1 Er Doped LiNbO3 ............................................................................................... 117 4.6.2 Fe Doped LiNbO3 ............................................................................................... 121 5 POINT DEFECTS IN CERIA BASED ELECTROLYTES .................................................. 124 5.1 Solid Oxide Fuel Cell Electrolytes ................................................................................. 124 5.2 Structure of Ceria Based Materials ................................................................................ 125 5.3 Computational Details .................................................................................................... 128 5.4 Interatomic Potentials Evaluation .................................................................................. 130 5.4.1 Lattice Constant and Thermal Expansion ......................................................... 130 5.4.2 Chemical Expansion ........................................................................................... 132 5.4.3 Dielectric Properties and Oxygen Migration Energy ....................................... 134 5.4.4 Elastic Constants ................................................................................................. 135 5.4. 5 Potential Summary .............................................................................................. 137 5.5 Mechanical Softening Effects ......................................................................................... 138 5.6 Oxygen Diffusion ............................................................................................................ 143 5.7 Summa ry .......................................................................................................................... 146 6 CASCADE SIMULATIONS IN TITANIUM ........................................................................ 148 6.1 Background ...................................................................................................................... 148 6.2 Structure o f Titanium ...................................................................................................... 150 6.3 Computational Details .................................................................................................... 152 6.3.1 MEAM Potential ................................................................................................. 154 6.3.2 ZBL Potential ...................................................................................................... 159 PAGE 7 7 6.4 Defect Analysis Method ................................................................................................. 160 6.4.1 Common Neighbor Analysis .............................................................................. 161 6.4.2 Lattice Matching Analysis ................................................................................. 162 6.5 Single Crystal Results and Discussions ......................................................................... 163 6.6 Effects of Grain Boundaries ........................................................................................... 172 7 SUMMARY AND FUTURE WORK ..................................................................................... 175 7.1 Lithium Niobate .............................................................................................................. 175 7.2 Ceria based Materials ..................................................................................................... 178 7.3 Titanium ........................................................................................................................... 179 LIST OF REFERENCES ................................................................................................................. 180 BIOGRAPHICAL SKETCH ........................................................................................................... 194 PAGE 8 8 LIST OF TABLES Table page 3 1 The comparison of bond lengths and coordination numbers between Li2O, Nb2O5, and LiNbO3. ............................................................................................................................ 53 3 2 Lattice parameters calculated by PAW LDA and PAW GGA, compared with previous results.. ..................................................................................................................... 54 3 3 Di fference in DFE, and relative energy difference, between 2 x 2 x 1 system and 2 x 2 x 2 system for various defects. ........................................................................................... 58 3 4 The heats of formation of Li2O, Nb2O5, and LiNbO3 calculated from the GGA calculations, com pared with experimental results ................................................................ 62 3 5 The en ergy of lithium metal, niobium metal, oxygen gas, Li2O, Nb2O5, and LiNbO3 as calculated by GGA. ........................................................................................................... 63 3 6 Change in oxygen chemical po tential with respect to 0 K value from experiment 130. ..... 64 3 7 Defect formation energies for various defects, compared with previous electronic structure and atomistic simulations. ...................................................................................... 72 3 8 Association energy of Li2O Pseudo Schottky, Nb2O5 PseudoSchottky, 5 NbLi +4 VNb and NbLi + 4VLi .................................................................................................... 75 3 9 Defect formation energies and polarization of various A3 configurations.. ...................... 79 3 10 Defect formation energies and polarizatio n of various A4 configurations ......................... 81 3 11 Activation Energy (Ea) obtained from first principles calculations, compared with the experimental results. .............................................................................................................. 85 4 1 Wyckoff position and fraction coordinates of Er2O3 ........................................................... 95 4 2 Wyckoff position and fraction coordinates of Nd2O3 .......................................................... 94 4 3 Fundamental properties comparison between GGA, GGA+U and experiments for FeO. ......................................................................................................................................... 98 4 4 Fundamental properties comparison between GGA, GGA+U and experiments for Fe2O3. ...................................................................................................................................... 98 4 5 Fundamental properties comparison between GGA, GGA+U and experiments for Nd2O3. ..................................................................................................................................... 99 PAGE 9 9 4 6 Defect formation energies for defect reactions under different reference states. The chemical potential of Er is assumed to be the same as the value in Er2O3. ...................... 107 4 6 Charge transfer level of impurities calculated by GGA with comparison of experiments. .......................................................................................................................... 116 5 1 Potential Fluorite structure data from International Tables of Crystallography, Volume A: space Group Symmetry, Fourth Edition. ........................................................ 127 5 2 Atomic position of ceria from neutron diffraction ........................................................... 127 5 3 Parameters of the five Buckingham potentials for ceria based systems. .......................... 131 5 4 Parameters of the Inaba 220 Born Huggins Mayer potential for CeO2. No shell model is available for this potential. ............................................................................................... 131 5 5 Properties of CeO2 from experime nt and from various interatomic potentials. ............... 134 5 6 Comparison of elastic properties of CeO2 from experiment, DFT calculations and a s predicted by empirical potentials.. ...................................................................................... 136 5 7 Anisotropy in the Youngs modulus of CeO2 using different methods. Numbers in brackets indicate fractional value relative to that in the <001> direction. ....................... 137 5 8 Ionic Radi i for relevant ions from Shannon. ...................................................................... 141 6 1 The system size and number of simulations preformed for each PKA energy. The neutron energ y that is equivalent to average PKA energy is also listed. .......................... 154 6 2 The parameters of spline function for MEAM_II 83. ......................................................... 158 6 3 The fundamental properties calculated using MEAM_I and MEAM_II, compared with EAM, FS, DFT and experimental data. ...................................................................... 159 6 4 The fundamental properties calculated using MEAM_I and MEAM_II, compared with EAM, FS, DFT and experimental data. ...................................................................... 159 6 5 The predefined parameters of the ZBL potential 246. ......................................................... 160 6 6 Coordination number distributions provided by CNA with distinguishing HPC and FCC stacking sequence. ....................................................................................................... 162 PAGE 10 10 LIST OF FIGURES Figure page 1 1 The Gibbs free energy as fun ctional of defect concentration ............................................. 18 1 2 Diekes Diagram. ................................................................................................................... 22 1 3 Sugano Tanabe diagram ...................................................................................................... 23 1 4 Illustration of defect polarization int eracting with bulk polarization Ps is the polarization of the bulk system, and Pd is the polarization of defects. ............................... 25 1 5 Brouwer diagram of fluorite structure m ixed ionic electronic conductor ......................... 27 2 1 Illustration of the concept of pseudopotential. ..................................................................... 40 2 2 Cubic spline function of MEAM_II, which ar e fitted for hcp titanium metal. .................. 44 3 1 Schematic of possible defect arrangement involving one niobium antisite and four lithium vacancies. ................................................................................................................... 52 3 2 The density of state (DOS) and partial DOS (PDOS) using GGA for LiNbO3. ................ 56 3 3 Difference in defect formation energy (relative to the 1 x 1 x 1 system) as a function of system size.. ....................................................................................................................... 57 3 4 Stability range of chemical potentials (in eV) of the elements in LiNbO3.. ....................... 65 3 5 Defect formation energies of various p oint defects as a function of Fermi energy. .......... 68 3 6 Thermodynamic transition levels of O related and Nb related defects. The values under the line are with respect the valence band maximum (VBM). ................................. 69 3 7 Formation energies of neutral defects and defect clusters un der Nb2O5 rich conditions (blue) and Li2O rich conditions (red). ................................................................ 71 3 8 Defect stability range: in the region AEFD, the Li Frenkel i s the dominant defect reaction; whereas in EBCF, the niobium antisite compensated by lithium vacancies is dominant.. ............................................................................................................................... 73 3 9 Defect formation energies dependence on the number of lithium vacancies in the first nearest neighbor positions of niobium antisite.. ................................................................... 77 3 10 Local structure of niobium antisite. ...................................................................................... 78 3 11 Schematics of possible diffusion paths in the LiNbO3 ........................................................ 83 PAGE 11 11 3 12 NEB calculations of diffusion barriers for lithium vacancy for several diffusion paths.. ...................................................................................................................................... 84 3 13 The DFEs of model III, Li Frenkel and model I as a function of niobium chemical potential.. ................................................................................................................................ 88 4 1 Crystal structure of MgO and FeO. ....................................................................................... 93 4 2 Crystal structure of Er2O3. ..................................................................................................... 94 4 3 Crystal structure of Nd2O3. .................................................................................................... 95 4 4 Crystal Fe2O3. ................................................................................................. 96 4 5 Electronic density of States of Er2O3 calculated using GGA and GGA+U. The parameters U J = 4 eV, U J = 7.85 eV, and U J=10.3 eV are considered. ............... 100 4 6 Defect formation energies of the lowest energy charge state of ErLi and ErNb as a function of Fermi energy ..................................................................................................... 104 4 7 Partial density of state (PDOS) for ErLi .. and ErNb .......................................................... 105 4 8 Formation energies of Er related defect clusters under Li2O rich conditions (blue) and Nb2O5 rich conditions (red). ......................................................................................... 108 4 9 Possible lithium vacancy positions in the first nearest neighbor (FNN) around Er sitting on lithium site. The oxygen sub lattice is not shown. ............................................ 109 4 10 Defect formation energy for ErLi .. + 2VLi shows the effects of association of the lithium vacancies around ErLi .. site. .................................................................................... 110 4 11 Defect formation energies of various possible defect reactions for impurities under Nb2O5 reference state. .......................................................................................................... 112 4 12 Defect formation energies of various possible defect reactions for impurities under Li2O reference state. ............................................................................................................. 113 4 13 Schematics of calculations of charge transfer levels. ........................................................ 115 4 14 The relative position of Fe2+/Fe3+ charge transfer level with regard to the band gap. ..... 123 5 1 Electrical conductivity of various electrolytes as a functional of temperature ............... 125 5 2 Schematics of fluorite structure for CeO2. .......................................................................... 126 5 3 Kagome net formed by Ce from <1 1 1> view. Only the Ce sublattice is shown. .......... 127 PAGE 12 12 5 4 Temperature dependence of the lattice parameter of CeO2 obtained from MD simulation using different potentials compared with experiment data ........................... 132 5 5 The effect of composition on the lattice constants of CeO2x obtained from MD simulations at 1073K. The values are compared with exper imental results ................... 133 5 6 Composition dependence on the room temperature Youngs modulus of CeO2x obtained from MD simulation and from experiment ........................................................ 139 5 7 The effects of temperature on the normalized bulk modulus (B), shear modulus (G), Youngs modulus (Y), and anisotropic factor (Z) for CeO2 .............................................. 140 5 8 The effect of dopant ion size on the lattice parameter of Ce0.8M0.2O1.9 (where M = In3+, Y3+, Gd3+, Ce3+, La3+) obtained fro m MD simulation, compared with the experimental results ............................................................................................................ 141 5 9 Dopant size dependence on the Youngs modulus of Ce0.8M0.2O1.9 (where M = In3+, Y3+, Gd3+, Ce3+, La3+) obtained using Grimes Potential. ................................................... 143 5 10 Diffusivity dependence on temperature for various dopedceria systems. ....................... 144 5 11 Activation energies of oxygen diffusion of various ceria based electrolytes. The results are compared with experimental values and DFT results. ..................................... 145 5 12 Migration barrier along [001], [110], and [111] directions obtained from NEB ca lculations. .......................................................................................................................... 146 6 1 Phase diagram of titanium ................................................................................................. 150 6 2 ............................................................................... 151 6 3 S ........................................................................................ 152 6 4 The evolution of number of Frenkel defect with time in cascade simulations at 100K with PKA = 1 keV using different time steps. ................................................................... 153 6 5 The schematics of lattice matching analysis (LMA). Both interstitials and vacancies are shown. ............................................................................................................................. 163 6 6 Threshold energies at various crystallographic orientations using MEAM_I, compared with previous results using FS potential. .......................................................... 164 6 7 Snapshot of cascade simulation in single crystal system at different simulation time. Only the atoms with coordination rather than 12 are shown. ............................................ 166 6 8 The number of survival Frenkel defects as a function of simulation time with various directions, positions, and energies of primary knockoff atom. ........................................ 167 6 9 The effect of PKA directions on survival Frenkel defects. ............................................... 168 PAGE 13 13 6 10 The comparison of number of survival Frenkel defects as a function of PKA energies between MEAM_I, MEAM_II, and FS potentials.. ........................................................... 169 6 11 The potential energy as a functional of nearest neighbor distance using FS, MEAM_2006 and MEAM_2008 potentials. ...................................................................... 171 6 12 Polycrystalline structure of titanium in the cascade simulations. Courtesy of D. H. Kim. ....................................................................................................................................... 172 6 13 Snapshot of cascade simulation of polycrystalline system when PKA = 10 keV at 100K. Only the atoms with coordination rather than 12 are shown. ................................ 174 PAGE 14 14 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 FROM ELECTRONIC STRUCTURE OF PO INT DEFECTS TO PHYISCAL PROPERT IE S OF COMPLEX MATERIALS USING ATOMIC LEVEL SIMULATIONS By Haixuan Xu May 2010 Chair: Simon R. Phillpot Major: Materials Science and Engineering P oint defects play a significant role in determin ing the physical and chemical properties of materials Atomic level simulation is a powerful tool to investigate and characterize the effect of these point defects In this study, various aspect of t he structure and stability of complex materials have been determined and predicted for lithium niobate, ceria based systems, and titanium The production, evolution, and dynamic behavior of defects have been explored. The focus has been on establishing the relationship between point defects and fundamental properties of bulk materia ls Lithium niobate is an important ferroelectric and non linear optical material. For lithium niobate t he dominant defects with the lowest formation energies and their equilib rium structures are predicted under various experimental relevant environments The site preferences with corresponding charge compensation mechanisms are compared with experimental observation s The diffusion mechanism and energy barrier are determined to elucidate the dynamic behavior of defect and defect clusters. The effects of point defects on the polarization of the system are also discussed. PAGE 15 15 Ceria based systems are considered as potential electrolytes of solid oxide fuel cells. In ceria based system s the effects of sub stoichiometry, temperature and ionic radii on the mechanical properties are evaluated using molecular dynamics simulation. It is observed that sub stoichiometry lead to a significant softening of the elastic constant s Similar results are predicted for doped ceria systems. Th ese softe ning effects arise from the significantly reduced strength of ionic interactions Titanium is a candidate material for cladding of fast nuclear reactor system due to its high corrosion resistance and exc ellent mechanical properties. In this study cascade simulations are carried out to investigate its radiation resistance. The effect of a high energy atom (primary knock on atom) is simulated with various energies, positions and orientations A high disordered region with a large number of point defects is observed during the initial phases of simulations (ballistic phase), followed by recombination of interstitials and vacancies (relaxation phase). The effects of primary knock on atom energies on rem nant defects are established The orientation effects of primary knock on atom and the effects of grain boundaries are also evaluated. PAGE 16 16 CHAPTER 1 INTRODUCTION 1.1 Point Defects in Crystalline Materials T he long range order of atoms is the defining featur e of crystalline materials However, most crystalline materials are not perfect: the atomic arrangement can be interrupted and altered by various types of defects. According to their geometry and shape, th e defects can be classified int o point, line and planar defects 1, 2 For instance, p oint defects can be a regular lattice point not occupied by the proper ion Linear defect s, such as dislocations, are defects that some of the atoms in the crystalline materials are misaligned. Planar defe cts are surface imperfection in polycrystalline materials, such as grain boundaries or surfaces It has been shown that properties of materi als are strongly affected by the presence of these defects 1. Therefore, it is of great importance to establish the relationship between these defects and the physical properties of materials. This study will focus on point defects. There are many types of point defects 1. Vacancy : Site usually occupied by atom or ion is not occupied. Interstitial: Atom is at a location in the crystal lattice that is usually empty. Antisite : Site usually occupied by a specific type of atoms or ions (A) is occupied by a different type of atoms or ions (B). Substitutional impurity : An atom, which is of the different type than any of the bulk atoms of the materials, replaces one of the bulk atoms in the crystal lattice. Vacancies, interstitials, and antisites are considered to be intrinsic point defects, while subsituti onal impurities are considered to be extrinsic defects. C ombination of these point defect s form s defect complex es such as Schottky defects and Frenkel pairs PAGE 17 17 A Schottky defect is a combination of different types of vacancies that are on different sublatti ce 1. For Schottky defect s, the numbers of vacancies formed on each sublattice need s to satisfy the requirement of charge balance For example, in a AB O3 oxide with A +1 charged and B +5 charged the Schottky defects can be expressed as O B AV V V Null 3' ' [1 1 ] By contrast a Frenkel pair is a combination of a vacancy and interstitial of the same species1, w hich can be formed when an atom move s from its original crystal lattice site to a nearby interstitial site, creating a vacancy behind. The Frenkel pair s can occur on both cation and ani on sublattice. The Frenkel pairs on anion sublattice are usually represented as a nti Frenkel defects. In the ABO3 oxide mentioned above the Frenkel defects can be: A iV A Null [1 2] In addition, other types of defect complexes could form in materials. For example, in a ABO3 oxide, an antisite can be formed when a B ion occupy a A site ( BA .). This antisite needs to be compensated by defects with opposite charges to maintain the overal l charge neutrality either locally or globally such as by vacancies of B cations. The possible defect reaction is: ' '4 5B AV B Null [1 3] In specific materials, all the above defect complexes are possible. However, only one will prevail, depending on the energetics of these defect reactions. 1.1.1 Energetics and Diffusion of Point Defects Energetics: The competition between the energy o f creating a point defect (enthalpy change H) and the entropy increase s leads to a total Gibbs free energy decrease in the materials. In the dilute limit, t he Gibbs free energy change ( G )as a functional of defect number can be expressed by 1: PAGE 18 18 ) ln ln ( ln N n n n n N N N kT kTn h n Gd d v d d d d [1 4 ] w here n is the number of defects, N is the number of the atoms involved, k is the Boltzmann constant, T is the temperature, hd is the enthalpy of formation v is vibration frequency when atoms are at equilibrium structure, and v is the vibration frequency when the defects are present in the nearest neighbor positions. As the equation shows, the Gibbs free energy change is a function of both temperature and nd. If T is kept constant, the dependence of Gibbs free energy on nd is shown in Figure 1 1 Figure 1 1. The Gibbs free energy as functional of defect concentration 1. The defect concentration, at which Gibbs free energy is minimum, is the equilibrium defect concentration ( neq). neq can be obtained using the following standard formula 1: PAGE 19 19 ) exp( kT S T h N nd eq [1 5 ] As mentioned previously, the equilibrium concentration of defects is also a functional of temperature Increas ing the temperature lead s to an exponential increase in the number of the defects in the system. Diffusion: The motion of these defects could happen when the energy of an atom is large enough to overcome the energy barrier of diffusion. The diffusivity of vacancy can be calculated by 1: ) exp(2T k H DB m o vac [1 6 ] w here is a geometric constant depending on crystal structure is the jump distance, 0 is the vibration frequency of the atoms, H is the migration energy, kB is the Boltzmann constant, and T is the temperature. Since the multiplication of and is unity, the vacancy diffusivity solely depends on the jump frequency and the distance of the jump. In the vacancy mechanism the atom and vacancy swap lattice site, which leads to: vac vac ion ionc D c D [1 7 ] w here Dion is the diffusivity of ions, cion is the ion concentration, Dvac is the diffusivity of vacancy, and cvac is the concentration of vacancies In general, cvac is very small, and cion can be approximate d as unity. The equation can be simplified as vac vac ionc D D [1 8 ] This indicates that the diffusivity of ions is much smaller than vacancy diffusivity. The vaca ncies move often, however, they are not numerous. By contrast atoms move less frequently, PAGE 20 20 but there are a lot of them. Dion is the self diffusivity, which can be calculated from mean square displacement (MSD) of the ions: 2  6 1o i ionr r t D [1 9 ] w here ri is the position of ion i at time t and r0 is the original position of ions. The i nterstitial mechanism involves defect jump s from one interstitial si te to another interstitial site The interstitial diffusivity can be calculated by 1: ) exp(2 int intT k H DB m o [1 10] From comparison between Equation 1 6 and 1 10, it can be seen that there is little difference between the formulas for vacan cy and interstitial diffusivity, e xcept the successful j ump frequency and jump distance. In addition, interstitial diffusivity is affected by a geometry factor, which depends on the structure of the materials. However, it is noted that the migration e nergy of interstitial is generally different from that of vacancies. 1.1.2 Electronic Structure of Point Defect The electronic structure and properties of a material are modified by the presence of point defects. The interaction between point defects and bulk materials sometimes give unique feature of the ma terial, such as the F center defects (color center) in alkali h alides 3. This is because point defects can introduce new energy levels into the band gap of the host crystal, which leads to transition s involving these levels thereby producing electronic features not present in the perfect materials. Many practical applications are based on the deliberate incorporatio n of impurities in the system, such as introducing rare earth elements and transition metals. Rare ear th elements have unfilled 4f electrons screened by fully filled 5s and 5p shells. Transition metals have unfilled 3d PAGE 21 21 outer shells. The electronic configurations of rare earth elements and transition metals give them different feature s from bulk and lead to various application s such as solid state laser, amplifiers, and superconductors 4 The rare earth element s sometimes called lanthanide ions, have an electron configuration of 5s25p64f n, where n v aries from 1 (Ce) to 13 (Yb). The valence electrons of these elements are from 4f orbitals. However, they are shielded by 5s and 5p outer electrons, which usually make them weakly affected by the host ions. Therefore, the Stark energy level s of rare earth elements in crystal are very similar to those of free ions 4, which also leads to the energy levels similar from one material to another. In fact, the energy level of rare earth elements provided by Die ke 5 based on particular ho st (lanthanum chloride) can be used for a wide range of materials. This energy level diagram, the socalled Dieke diagram 5, is shown in Figure 1 2. The diagram shows the energy of 2S+1 LJ states for rare earth elements in LaCl3, in which S is the total spin, L is the orbital angular momentum, and J is the sum of L and S The width of each state indicates the magnitude of the crystal field splitting, which the center of each multiplet gives the approximate location of its corresp onding free ions energy level 4. Transition metal ions, on the other hand, have an unshielded 3d outer shells. Therefore, there is a strong interaction between the ions and host lattice. Due to this s trong interaction, the energy levels of transition metal ions in one material may be very different from another 4, depending on the crystal field. Sugano and Tanabe have calculated the energy of the st ates for 3dn ions (n=2 to n=8) as a function of crystal field strength, generating the socalled Sugano Tanabe diagram 6, 7 which is quite useful in understanding these ions in a variety of materials (Figure 1 3). The x axis of the diagram is expresses in terms of ligand field splitting parameter, divided by the Racah parameter B. The y axis is expressed in terms of energy (E) divided by PAGE 22 22 B. Sugano Tanabe diagram shows how the 2S+1 LJ levels split up as the ration between the crystal field strength and the interelectronic interaction (Dq/B) increases, w hich allows deducing useful information about the nature of the optical bands of transition metal ions 4. Figure 1 2. Diekes Diagram 5. PAGE 23 23 Figure 1 3. Sugano Tanabe diagram 6, 7. 1.1.3 The Effects of Point Defects on Properties of Materials When defects are introduced into the system, the local elastic and electrostatic environments change. The elastic interaction comes from the mismatch in size between defects and host ions, while the electrostatic interaction is from the charge difference. Due to this change of elastic and electrostatic interactions, many properties of the system can be affected by the PAGE 24 24 presence of point defects, including mechanical properties, thermal properties, electrical and magnetic properties, ferroelectricity, and opt ical properties. Mechanical properties: The mechanical behaviors consist of elastic and plastic part s Both plastic and e lastic properties could be significantly affected by point defects. This is because that elastic modulus is a measure of the interatom ic force or bond strength. The elastic modulus in a simple model can be expressed as the first derivative of the interatomic force The presence of point defects can modify the interatomic forces, and hence the elastic modulus. In some cases, the influence of the point defects on elastic modulus could be significant. For instance, it has been shown experimentally that elastic modulus of TiNx decreases substantially as the function of the concentration of nitrogen vacancies increase 8, 9. The f racture toughness which is the function of both Youngs modulus and surface energy, will also be affected by the point defects thro ugh the effect in the Youngs modulus. For example, it has been found that the fracture strength and fracture toughness of doped ceria appear to decrease with increasing the dopant concentration when the dopant level is less than 20 mol %. Further studies on the fracture properties of rare each doped ceria showed that the fracture toughness was influence by the dopant concentration rather than the dopant species Since the oxygen vacancy concentration is directly associated with dopant concentration in dope d ceria system, the results indicates the critical role of oxygen vacancy on the fracture properties. In summary point defects play a n important role in determining the mechanical properties of the materials Ferroelectricity: A ferroelectric material pos sesses a spontaneous polarization, which can be reversed by an external electric field. In ceramic materials, several polarization mechanisms exist, such as electronic polarization, ionic displacement polarization, and space charge PAGE 25 25 polarization. The electr onic polarization is caused by the displacement of the electron cloud relative to its nucleus. For the ionic polarization, the effective cation center is separated from the effective anion center by a certain distance 1. F erroelectric materials can form polarization domains, res ulting in a domain wall between two domains with opposite polarization directions However, the presence of point defects may significantly change the local polarization, domai n wall structure, and domain wall motion. The arrangement of the point defects may not preserve the long range symmetry of the material, therefore, forming a local polarization itself A n illustration of how the defect polarization may interact with bulk polarization is given in Figure 1 4. The direction of local polarization does not necessarily have to b e the same as bulk polarization, i n most cases, it is different 10. Therefore, the interaction between polarization due to point defects and bulk polarization could change the macroscopic properties of the materials, such as coercive field for domain reversal. Figure 1 4. Illustration of defect polarizat ion interacting with bulk polarization 10. Ps is the polarization of the bulk system, and Pd is the polarization of defects. When point defects are located at domain walls, the local structure of domain wall will be changed as point defects can introduce strain and electrostatic interactions into the system PAGE 26 26 Furthermore, there could be a pinning effect on domain wall m otion from the point defects. This is because point defects have to rearrange themselves during the domain wall motion. However, this arrangement process involves the migration of the point defects. The energy barrier for migration may be substantial in s ome materials making the rearrangement process very difficult Therefore, point defects can greatly limit the motion of domain walls. Ionic Conductivity: The ionic conductivity in materials can also be tuned by point defects. In some cases, the ionic conductivity is obtained through defects that are thermally generated, such as Frenkel pairs or Schottky defects. In other cases, impurities are intentionally added to control the ionic conductivity of the materials. In the latte r, the external impu rities have to be in equilibrium with intrinsic point defects since t he defect reaction has t o satisfy the requirements of charge balance, mass balance and site balance. Therefore the concentration of one type of defect s can be affected by the concentrat ion of another T he ionic conductivity is also a function of defect concentration. For example, as the defect concentration changes, the pre factors are different. Furthermore the activation energy of ionic conductivity consists of migration energy and a ssociation energy. The association energy take s into accoun t the defect defect interaction; therefore, it strongly depends on the defect concentrations. T he ionic conductivity can be related with diffusivity through Nernst Einstein equation, which is expr essed as 1: kT D C e Zion ion ion 2 2 [1 1 1 ] w here ion is the ionic conductivity, Z is valence of ion, e is the charge of electron, Cion is the concentration, Dion is the diffusivity of ions, k is the Boltzmann constant, and T is the temperature. PAGE 27 27 In experiments, desired conductivity is usually achieved by controlling the temperature and other thermodynamic parameters, such as oxygen partial pressure or impurities In such system, many defect reactions are possible and multiple defect equilibrium will occur. The dependence of defect concentration of different species on the oxygen partial pressure is usually shown in a diagram, term as Brouwer diagram. The Brouwer diagram consider s the net effect when all of the possible defect reactions simultaneously, therefore, providing a convenient way to represent the variation in defect concentration with changes the activity of a component of the compound A exemplar case is given in Figure 1 5 for fluorite structure mixed ionic electronic conductors. Figure 1 5. Brouwer diagram of fluorite structure mixed ionic electronic conductor 11. 1.2 Defect Stability in Lithium Niobate LiNbO3 is one of the most technologically important oxides with an extraordinary combination of ferroelectric, piezoelectric, acoustic, and optical properties 10, 12. For a very long PAGE 28 28 time work has focused on the congruent composition of LiNbO3 1315 This is because only recently the stoichiometric LiNbO3 be synthesized u sing the double crucible Czochralski (DCCZ) 16, 17 and vapor transport equilibration (VTE) methods 18 20. This indicates that LiNbO3 has a strong tendency to form a lithium deficient structure. However, the reason why stoichiometric LiNbO3 is difficult to obtain is not yet clear The change of composition from the congruent (lithium poor) composition to the stoichiometric composition results in significant changes in the physical properties of the system, including the Curie temperature 21, the ferroelectric coerciv e field 21, 22, and the photorefractive properties 10. These changes in physical properties are thought to be due to the point defects and defect clusters. Yet, very little is known about the detailed str ucture and energetics of such clusters. Three models for intrinsic point defect structures in LiNbO3 have been proposed. In the first model (Model I), two lithium vacancies are c ompensated by an oxygen vacancy 23. In the second model (Model II), n iobi um antisites are compensated by niobium vacancies 24. In the third model (Model III) niobium antisite s are compensated by lithium vacancies 10. For each model, different configu rations of the defects are possible. Nevertheless, no fundamental understanding of either the local structure and the se effects or their effects on the material properties is available. In addition, t o achieve the desired functionalities, various dopants are added to the system. For example, Mg is introduced to increase the resistance to photorefractive damage 25, 26; Fe doping is used for applications such as holographic storage and beam coupl ing 27; Nd is added for solid state laser applications 28 31 In addition, dopant ions have also been employed as PAGE 29 29 probes to investigate the structure of domain wall s and defect/domainwall interactions 3234. However, little is known about the incorporation site preference of these impurities. To investigate the s ite selectivity of dopant ions, various experimental techniques, including electron spin resonance 4, electron nuclear double resonance (ENDOR) 35, Mossbauer spectroscopy 36, 37, Rutherford backscattering 38, X ray standing wave analysis (XSW) 39, 40 extended xray absorption fine structure (EXAFS) 41, and p roton i nduced X ray e mission (PIXE) 42 have been employed. Moreover, several optical and magnetic resonance spectroscopy studies 4347 have been used to determine the local env ironments and configurations around the dopant sites. However, a general consensus has not been reached as to what the dominant dopant configurations are. Furthermore some of these studies used congruent LiNbO3 samples grown by the Czochralski method 13, 14, others used stoichiometric LiNbO3 produced through vapor transport equilibration (VTE) 1820. This makes it very difficult to compare results of experiments on samples from made with different synthesis technique; as a result the influence of the sample stoichio metry on the site selectivity and distribution is unclear. 1.3 Mechanical Properties of Ceria based Materials The first test of a solid oxide electrolyte for solid oxide fuel cells (SOFCs) was carried out in the 1930s Since then, much research has focused on finding an appropriate electrolyte 48, 49. To be used in practical applications, solid electrolytes need to satisfy several major requirements, including a high ionic conduction, negligible electronic conduction and thermodynamic stability over a wide range of temperature s and oxygen partial pressure s In addition, solid electrolytes must have a thermal expansion compatibility with electrodes and suitable mechanical properties 50. Among various electrolytes, yttria stabilized zirconia (YSZ) 5154 has been investigated and used extensively. Compared with other solid electrolytes, YSZ exhibits very small electronic PAGE 30 30 contribution to total conductivity in the oxy gen partial pressure range most important for practical applications 50. However, YSZ electr olytes require a very high operating temperature (~1000 oC), which could cause many interface reactions such as the interaction between solid electrolytes and electrodes 55. Thus, research focus on other potential electrolyte materials is of great interest. As an alternative to YSZ, doped ceria 5659 h as been the subject of a great deal attention due to its higher ionic conductivity than that of YSZ. The electrical co nductivity of Ce0.80Gd0.20O1.9 (gadolinia stabilized ceria) is about one order of magnitude higher than that of YSZ in the intermediate temperature (IT) range (500oC 700oC). Although Bi2O3based materials showed higher ionic conductivity than YSZ and doped ceria, these materials are unstable at high temperatures and low partial pressures 60, which severely limits its application. The higher ionic conductivity of ceria based electrolytes with respect to YSZ is th e main advantage for use in SOFCs. However, doped ceria systems undergo reduction under low oxygen partial pressures, and the magnitude of electrical conductivity and stability are greatly dependent on the type and concentration of dopants 55. Overall, ceria based electrolyte s show significant potential for practical application. Due to the presence of oxygen vacancies in the doped ceria structure and important role of defect in determining physical properties of ceramics, the mechanical properties are quite different from that of the perfect structure. T he mechanical property of doped ceria electrolytes are found to be inferior to YSZ 52. Furthermore, the n anoindentation experiments indicate the elastic modulus decrease s with increasing the oxygen vacancy concentration However, the reason for this is yet fully understood. PAGE 31 31 1.4 Radiation Damage in Titanium Titanium metal is widely used in industry, recognized for its high co rrosion resistance and high strength with relatively less weight Titanium is c apable of withstanding acids and moist chlorine Moreover, it is a very strong material with high ductility. The tensile strength of titaniu m is compa rable to that of low strength martensitic stainless steel s and is better than that of austenitic or ferritic stainless steels Due to t he relatively high strength and exceptional corrosion resistance, titanium has been considered as a potential material f or nuclear reactor cladding. The degradation of thermo mechanical behavior of cladding materials in the nuclear industry is a longstanding issue. The cur rent cladding material is zircaloy. However, their corrosion resistance and ductility are poor compared with titanium. The radiation process creates numerous point defects within a very short time, generating a dramatic thermal gradient and pressure gradient The local structure is significantly modified and may be very unlike of the bulk structure. In a ddition, the defect generation has significant effects on the subsequent defect evolution and diffusion processes, substantially modifying the properties of materials. Therefore, it is critical to understand the effects of radiation on the structure and m echanical properties of titanium. In addition, t here are some indications that reducing the grain size will lead to an increase of resistance to the radiation damage 61. However, the role of grain boundaries during the radiation process is poorly understood. The presence of grain boundaries may change the defect evolution process and subsequently change the mechanical behavior of the materials. Consequently, establishing the coupling between radiation damage, defect evolution, and mechanical properties is of great importance. PAGE 32 32 1.5 Overview of This Dissertation This work is focused on the point defects energetics, dynamics, and their effects on the physical properties of lithium niobate, ceria based materials, and titanium using density functional theory and molecular dynamic simulations. In lithium niobate, the objective is to determine the dominant intrinsic defects/clusters under various experimental relevant conditions such as for congruent and stoichiometric composition, as well as their local structures and configurations. Furthermore, the effects of such point defects and defect clusters on the polarization and ferroelectricity of the system will be explored. To determine whether th e dominant defects are static or dynamic stable, the energy barrier of diffusion and corresponding mechanism will be elucidated. The details of intrinsic defects in lithium niobate are described in Chapter 3 In addition to intrinsic point defects, various impurities have been introduced into the material to achieve functionalities. Therefore, the understanding of the incorpora tion site preference for different dopants and their charge compensation mechanism s is requisite. Furthermore, as impurities may introduce additional levels within the band gap of the bulk materials, the charge transfer levels of impurities need to be pred icted to demonstrate the effects of impurities on the optical properties of the materials The details of extrinsic defects in lithium niobate are given in Chapter 4. In ceria based materials, the objective of this work is to understand the effects of oxy gen vacancies on the mechanical properties using molecular dynamic simulation. To achieve this, several empirical potentials have been evaluated first. Moreover, the effects of the dopant ionic radii on the mechanical properties will be established in orde r to search the optim al ionic radius of the dopants ( Chapter 5 ). PAGE 33 33 For titanium metal, radiation effects on the materials will be investigated using cascade simulation s in order to understand the microstructure evolution of the materials under radiation da mage. The survival Frenkel defects will be analyzed and their relationship with the primary knock on atom energy will be investigated and compared with previous studies Furthermore, the role of grain boundaries during the radiation process will be elucida ted. The details of radiation damage in titanium are discussed in Chapter 6. The summary of this work and envisioned fut ure work are presented in Chapter 7. PAGE 34 34 CHAPTER 2 SIMULATION METHODS Electronic and atomic level simulations have been used extensively to study various properties and phenomena of materials, as they can provide information with atomic resolution that sometime is difficult to obtain from experiments. In the present study, electronic structure calculation at the level of density functional theory ( DFT) 62, 63 and atomistic simulation at the level of molecula r dynamics (MD) 64 are employed DFT is one of the most powerful and successful approaches to investigate the electronic structure in condensed matter materials In numerous cases, the results predicted from DFT agree quite satisfactorily with the experimental results and observations 63. Furthermore in principle, DFT is an approach that does not require empirical parameter s In practice, certain approximations are needed. Therefore, due to this high fidelity, DFT is considered as a powerful t ool in solid state physics, chemistry and materials scienc e In general, DFT currently can only deal with up to hundreds of a toms and thousands of electrons. This prohibit s many phenomena and properties at the larger scale s from being considered, such as dislocations phase transformation s MD, on the other hand, can easily deal with millions of atoms at finite temperature using todays computational power Moreover MD can be used to study the nonequilibrium dynamic process of the system, which significantly broadens its applications. U nlike DFT, no electronic information is needed or provided by MD Rather a classical interatomic potential 64 with a set of empirical parameters is used in order to describe the interactions between ions. In general, the se empirical parameters are fitted to specific system s which lead s to poor transferability of the parameter and less material fidelity compared with DFT As both DFT and MD have their PAGE 35 35 strength s and limitations, discretion is needed in choosing the appropriate method for a specific problem or process. 2 .1 Density Functional Theory 2.1.1 Schrodinger Equation The quantum state of a system consisting of electrons and nuclei is described using Schrdinger Equation. The time independent, nonrealistic Schrdinger Equation 65 is: E H [2 1 ] where H is the Hamiltonian operator, is the eigenstate wave function and E is the eigenvalue For atoms with electrons and nuclei, H can be expressed as: J I J I J I I I I j i j i I i I i I i i eR R e Z Z M r r e R r e Z m H2 2 2 2 2 2 22 2 1 2 [2 2 ] where, me is the mass of electron MI is the mass of nuclei. The H consists of five terms, which are kinetic energy of electron s the electron nuclei interactions, the electronic electron interactions, kinetic energy of nuclei, and interactions between nuclei and nuclei, r espectively. However, analytical solutions are only available for very few simple cases, such as the one dimensional potential well or the hydrogen atom 62. For complex systems, only approximate solution s are possible. To obtain such approximate solutions, several assumptions have to be made. 2.1.2 Born Oppenheimer Approximation The Born Oppenheimer approximation 66, is based on the observation that electrons move much faster and weigh much less than the nuclei. It assumes the electronic motion and nuclear motion can be decoupled and that the electrons are in equilibrium with nuclei Thi s approximation is typically made although not all the methods employ it. In the Born  PAGE 36 36 Oppenheimer approximation, the electronic wavefunction depends on only the nuclear position but not the nuclear velocity. The Schrdinger Equation for electrons can then be written as: E r r e R r e Z mj i j i I i I i I i i e 2 2 2 22 1 2 [2 3 ] w here only the first three terms in equation 2 2 remain which are the kinetic energy of electrons the electron nuclei interactions, and the electronic electron interactions. The kinetic energy of electrons and electron nuclei term can be easily calculated. However, it is a tremendous challenge to describe the interaction between electrons and electrons, which is an interacting many body problem. All the above equations ( 2 1, 2 2, and 2 3) are ex pressed in term of wavefunctions. From wavefunctions, the probability of finding an electron at a given place or the electron density at a particular coordinate can be calculated by i i ir r r n* [2 4 ] where, the i indicates the complex conjugate of the wave function. The link between wavefunction and electron density makes the electron density a fundamental variable of the system 2.1.3 Hohenberg Kohn Theorem Two fundamental theorems in term s of electron density were established by Hohenberg and Kohn 67. Theorem I states that For any system of electrons in an external potential Vext(r), that potential is determined uniquely, except for a constant, by the ground state of electron density n(r) which established a one to one mapping between the external potential and th e ground state electron density and allow ed all properties of the system to be determined from t he ground state electron density. Theorem II states that A universal functional for the energy E[n] PAGE 37 37 of the density n(r) can be defined for all electron systems. The exact ground state energy is the global minimum for a given Vext(r), and the density n(r ) which minimizes this functional is the exact ground state density which is a direct consequence of variational principle. Although the Hohenber g Kohn theorem are very powerful, neither theorem define d the functional of the electron density or indicate d a systematic way of finding the functional 2.1.4 KohnSham Equation A route to actual calculations is provided by the Kohn Sham equations 68, which replace d the difficult the ma nybody problem with an auxiliary independent electron problem that can be exactly solv ed 62. Since there is no unique way to define the auxiliary system, it is usually rephrased as an ansatz The KohnSham equations assume the exact ground state density can be represented by the ground state electron density of a non interacting auxiliary system ; the auxiliary Hamiltonian is chosen to have the usual kinetic operator and an effect ive local p otential. Therefore, the KohnSham equation can be expressed as: r r r V m hi i i eff 2 22 [2 5 ] r V r V r V r Vxc Hartree ext eff [2 6 ] w here Vext(r) is the external potential due to the nuclei and any other external potential, Vext(r) is the Hartree potential, Vxc(r) is the exchange and correlation potential, i is the wave function of state i i is the Kohn Sham eigenvalue. The Hartree potential is defined as: 3 2) ( r d r r r n e r VH [2 7 ] The exchange and correlation potential is given as: PAGE 38 38 r n r E r VXC XC [2 8 ] where, EXC(r) is the exchange correlation energy. This exchangecorrelation functional is the only unknown term in the KohnSham equations. Therefore, reasonable approximation s such as local density ap proximation (LDA) and generalized gradient approximation (GGA), are requi red to solve the equations 62. 2.1.5 Exchange Correlation Functionals In LDA 62, the exchange correlation energy per electron at a given point r is assumed to be equal to the exchange correlation energy per e lectron in a homogeneous electron gas that has the same density at r Therefore, LDA does not consider the exchange correlation energy due to the inhomogeneities in system. Despite this, LDA work s remarkably well for many systems, even in very inhomogeneo us cases. The GGA on the other hand, includes not only the local density but also a density gradient term, which account s for the inhomogeneities of the system. There are many forms of GGA, such as Perdew and Wang (PW9 2 ) 69, Perdew, Burke, and E rnzerhof (PBE) 70, and Lee Yang Pa rr (LYP). In general, GGA predict s lower binding energies than LDA, which improves the agreement with experiments, as LDA is well known to overbind the system However, GGA does not lead to a consistent improvement over LDA, as it violates the sum rul e and other relevant constraints Sometimes, GGA even leads to worse results compared with LDA, making the selection of the appropriate functional a challenging task The wave functions ( ) in E quation 2 5 can be generally represented using either plane waves or localized orbitals. Each method would be mo re appropriate for a particulate range of problems and applications. Localized orbitals provide an intuitive picture that captures the essence of atomic like feature in solids, which is widely used in qu antum chemistry. PAGE 39 39 Furthermore, localized orbital s allow, in principle, linear scaling of DFT calculations with system size. However, the localized orbitals basis set depends on the positions of the ions, which needs a Pulay correction added to the Hellman Feynman forces 62. The completeness and superposition errors of localized orbitals basis set s are very difficult to control. In general plane waves provide general solutions of differential equations, such as the Schrdinger Equation. Therefore, plane waves are a natural c hoice for system with periodic boundary conditions (PBC) 64. They provide important information about the band structure s of many materials, including metal, semiconductors, and insulators. Moreover, the conversion between real and reciprocal space is easy and no Pulay correction is needed for the force calculations on atoms. However, t he converge nce of plane waves is very slow. To overcome this, the electron ion interaction are usually represented by pseudopotentials 2.1.6 Pseudopotential Approximation The key idea of pseudopotentials 71 is to replace the strong Coulomb potential nucleus and the effects of tightly bounded core electrons with an effective ionic potential that acts on the valence electrons and has the same scattering properties as the all electron potential beyond a given cut off radius (Figure 2 1) The approximation is based on the fact that most of the pr operties of the materials are only determined by the valence electrons. Compared with all electron method, the pseudopotential approach can significantly reduce the number of basis sets. The advent of norm conserving pseudopotentials (NCPPs) 72 and ultrasoft pseudopotentials (USPPs) 73, 74 has led to accurate electronic s tructure calculations that are the basis of modern research. The PAW 75 method introduces a projector acting on the valence electrons and auxiliary localize functions that like USPP method. It is a combination of both augmented wave method (APM) and the pseudopotential approach 76. Furthermore, it can be shown tha t NCPPs and PAGE 40 40 USPPs are approximations to the PAW method 77. In general, the PAW method is expected to be as accurate as all electron methods and as efficient as pseudopotentials approach. Figure 2 1. Illustration of the concept of pseudopotential. Picture taken with permission from Wikipedia: http://en.wikipedia.org/wiki/Pseudopotential 2 .1.7 Limitations In general, DFT is remarkably successful in determining a wide range of properties of the materials. However, there are certain limitations of the method. The best known problem of DFT is the underestimation of the band gap. The band gap of s emiconductor and insulator is always underestimated by DFT. This is because DFT is the ground state theory while calculating PAGE 41 41 the band gap involves excited states. In addition LDA usually predict s small lattice constants than experiments which result s in too large cohesive energies and too large bulk moduli. GGA, on the other hand, largely correct s this overbinding problem 63. However, GGA has a tendency to overshoot this problem for the heavy elements. Furthermore, DFT neglects strong correlation s and van der Waals i nteractions. Neglecting strong correlation cause s DFT fail to predict the correct band structure of Mott insulator with d or f orbitals 78. 2 .1.8 DFT +U To overcome these problems, many approaches have been developed, such as DFT+U 78, 79. DFT+U attempts to improve the band structure of a Mott insulator by adding an onsite repulsion term to the DFT Hamiltonian which result s in a localized feature of d or f electrons that DFT has failed to describe. Many version s of the DFT+U have been published in the literature. Among them, Dudarevs method 79 is probably the most widely used one due to easy implementation and lower computational effect requirement. In Dudarevs method, the total energy of the system is expressed as: 2 ,2m m DFT U DFTn n J U E E [2 9 ] where n is the occupation number of the m th d or f state. The U and J are empirical parameters. It is noted that in Dedarevs method, only the difference between U and J (U J ) is meaningful. 2 .2 Molecular Dynamics (MD) Simulation MD is routinely used to investigate the structure, dynamics, and thermodynamics of various materials. In MD, the motion of atoms in the system is predicted by solving Newtons law of motion, which is PAGE 42 42 i i im a F [2 1 0 ] where, Fi is the force acting on an atom i ; mi is the mass of atom i and ai is the acceleration of atom i The force on each atom can be calculated from the gradient of the potential energy: Vi i F [2 1 1 ] where, V is the potential energy. Therefore, to describe the interaction between atoms, an appropriate potential is needed. For different types of materials, the potentials that have been used to describe the system may be in different form s For any specific mate rial, even using the same form of potential, different parameterizations may be used. Therefore, it is of great importance to evaluate the potential parameter before any practical calculations have been performed. 2.2.1 Interatomic Interactions For ionic system s one of the most widely used potential is the Bu c kingham potential combined with Coulombic interac tions. The Buckingham potential is expressed as: 6expij ij ij ij ij ij Buckr C r A r V [2 12] where Aij, ij and Cij are the empirical parameters that are fitted to the system, and rij is the separation between two ions i and j Th e first term in Equation 2 12 is a repulsive interaction, which decreases exponentially with increase in distance between atoms i and j The second term in Equation 2 12 is an at tractive contribution due to the dipole interactions. The Coulombic interaction is a long range interaction and can be described as: N i N i j ij j i Coulr q q E1 1 ,2 1 [2 13] PAGE 43 43 where, N is the total number of ions in the system, qi and qj are the charges on atom i and j respe ctively The r1 dependence makes the sum conditionally convergent. The summation of the Coulombic term over the system can be carri ed out using either Ewald sum 64 or direct sum 80 The embedded atom method (EAM) was developed for metallic system 81. The potential energy of EAM consists of two parts: embedding energy and a twobody potential term. The embedding term, while the two body term depends on the interatomic bond distances. The total energy is expressed as 81: i i j ij ij i i iR F E) (2 1 [2 14] w here Fi is the embedding function, i is the background electron density at the site i ij is the pair function between atom i and atom j Rij is the distance between atom i and j The modified embedded atom method (MEAM) is a extension of the EAM to include angular contributions 82. In the present study, two types of MEAM potentials have been us ed. One is the traditional MEAM developed by Baskes 82. The details of this MEAM potential are given in section 6.3.1. Another type of MEAM potential is in the function form of cubic spline 83, t hereafter named MEAM_II. MEAM_II removes the constraints of fixed angular feature and allows additional flexibility of the potential. The total energy in MEAM_II is calculat ed by 83, 84 ij i i ijU r E ) ( ) ( [2 15] w here ij) is the pair interaction between atom i and j is the self energy function depending on the electron density of atom i The electron density of atom i is calculated using j jk jik ik ij ij ig r f r f r )] [cos( ) ( ) ( ) ( [2 16] w here jik is the angle between atom j i and k centered on atom i All five functions ij) f(r) and g [ ] are represented by cubic spline function s which are fitted to the PAGE 44 44 fundamental properties of the system. MEAM_II can be related to the Stillinger Weber potential [U(x)=x and =0 ] and EAM potential [ f=0 or g=0 ] 83. E xample s of five cubic spline functions are given in Figure 2 2 Figure 2 2 Cubic spline function of MEAM_II, which are fitted for hcp titanium metal 83. 2.2.2 Integration Scheme The integration scheme of Newtons equation of motion is a central part of the MD simulations. At a small increment of time ( t ), referred as the time step, the position, velocity and accelerations of the atoms are calculated base d on integration methods. There are several integration schemes avai lable, such as Verlet algorithm 64, which is one of the most widely used integration algorithms. In the current study, 5thorder Gear predictor corrector method 85 is employed. This is because the G ear algorithm is faster than the Verlet algorithm when the PAGE 45 45 number of the system is large 64. Moreover, the Gear algorithm is more accurate than the Verlet algorithm with a longer time step, therefore, achieving a higher degree of energy conservation 64. The central idea of the predictor corrector is that a Taylor expansion is used to predict the positions, velocities, accelerations, and higher order derivatives of position at time t + t then a corr ection term is added based on the predicted acceleration. For a fifth order predictor corrector, the p redictor formulas are given as 64: t d t t d t d t t c t t c t d t t c t t b t t b t d t t c t t b t t a t t a t d t t c t t b t t a t t v t t v t d t t c t t b t t a t t v t t r t t rp p p p p p 2 3 2 4 3 2 5 4 3 22 1 6 1 2 1 24 1 6 1 2 1 120 1 24 1 6 1 2 1 [2 17] w here r is the position v is the velocity and a is the acceleration of each atom in the system. b c and d are third fourth and fifth derivative s of the position with respect to time. T he superscript p represents predicted values. The total energy of the system can be calculated based on the predicted positions of atoms. The interatomic forces at time t + t are then calculated using the derivatives of potential energies with respect to the predicted atoms positions With information on interatomic force s accelerations are easily obtained. However, the se predicted acceleration values are based on Taylor expansion rather than phy sics Therefore, correction terms are needed in order to describe the physical behavior of atoms. T he corrected positions and other derivatives are calculated by 64: PAGE 46 46 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 18] T he superscript c corresponds to the correct ed values The coefficients for a fifth order predictor corrector are c0 = 3 / 16, c1 = 251/ 3 6 0 c2 = 1, c3 = 11/ 18, c4 = 1/6 and c5 = 1/60 86. These corrected values are then used to predict the positions and other higher order derivatives for the next iterati on 2.2.3 Thermostat T o simulate materials under realistic conditions, it is typically necessary to control the temperature. There are several commonly used approaches to achieve this. The simplest is probably the velocity rescaling thermostat 64. Within this method, the velocity of each atom, which is associated with the kinetic energy, is adjusted: ) ( ) (arg argt v T T t vcurrent current et t et t [2 19] w here Tcurrent is the current temperature, Ttarget is the target temperature, vcurrent is the current velocity, an d vtarget is the velocity after adjusting. However, this method has several problems 87. First the energy and momentum conservation are not obeyed. Secondly, the rate of temperature adjustment depends on the MD time step. Furthermore, the ergodic theorem will not hold, which could lead to a failure of prediction of physical properties based on the ensemble average 88. Due to these limitations, great caution is needed when velocity rescaling method is used. PAGE 47 47 In the current study of cascade simulations in titanium Berendsen thermostat is employed The Berendsen th ermostat 89 is found to be more reasonable in both the physics and statistical mechanics than velocityrescaling Within Berendsen thermostat, the velocity is adjusted by: ) ( 1 1 ) (arg argt v T T dt t vcurrent current et t et t [2 20] w here dt is the time step in the MD simulation, and is the adjustable parameter for the thermal equilibration. It is can be seen that the Berendsen thermostat is equivalent to velocity rescaling thermostat when = dt When is infinity, the velocity will not be adjusted, which corresponds to the microcannonical ensemble. In general, the Berendsen thermostat is gentler than velocity rescaling. The Berendsen thermostat can be considered as the first order approximation of t he Langevin thermostat 64. Ther efore, the dynamics of the system can be close to the realistic evolution. 2.2.4 Barostat In addition to temperature, controlling the supercell volume and shape is also requisite during the simulation under realistic conditions, such as under an external pressure. To achieve this, barostat schemes are need, such as Anders e n barostat 90 and Parrinello Rahman barostat 91. The idea of the Andersen barostat is to assume the system in a container that can be compressed by piston with a mass of M. Then the Lagrangian can be written as: N i N j i ij i i N NpV V M V U m V V V L1 1 2 3 / 1 2 3 / 22 1 ) ( 2 1 [2 21] w here i is the scaled coordinates and def ined as Cartesian coordinated divided by cubic root of the volume( V ), p is the external pressure acting on the piston. Then the motions of the atoms derived by Andersen can be expressed as: PAGE 48 48 dt V d r m P dt dri i i iln 3 1 [2 22] dt V d P r U r dt dPi ij ij iln 3 1 ) (' [2 23] V r U r m P p dt V Mdij ij i i/ ) ( 31 2 3 2' 2 2 2 [2 24] It is noted that Andersen barostat is isotropic, which means the external pressure is hydrostatic. In many cases, this is not a good approximation. The Parinello Rahman barostat extend s the Andersen method by allowing the shape of the simulation cell to change. This is achieved by introducing a 3x3 matrix H so that i iHs r H V ) det( [2 25] The Largrangian of Parrinello Rahman barostat can be expressed as 91: N i T ij i T i ipV h h MTr r U S G S m L1) ( 2 1 ) ( 2 1 [2 26] Both Andersen and Parrinello Rahman barostat s are used in the current study. 2.2.5 Ensembles By controlling temperature, pressure, and other thermodynamic properties, such as volume, total number of atom s and internal energy, simulations can be carried out under well defined thermodynamic conditions, such as in a well defined thermodynamic ensemble. S imulation of a system with fixed N, V, and E corresponds to the isochoric microcannonical ensemble (NVE) 64. Similarly, NVT is canonical ensemble and NPT is isothermalisobaric ensemble. These three are the most widely used ensembles in MD simulations, and have been employed in the current study. In addition, there are other types of ensembles, such as grand canonical ensemble, in PAGE 49 49 which the total number of the atoms (N) in the system can be changed while a fixed uniform chemical potential has to be maintained. PAGE 50 50 CHAPTER 3 INTRINSIC DEFECTS IN LITHIUM NIOBATE 3.1 Background of Intrinsic Defects LiNbO3 is an important ferro py ro and piezoelectric material with many promising physical properties 12. Its applications include use as a second harmonic generator, a parametric oscillator, a transducer and nonvolatile memory 18, 92. Nonstoichiometric, intrinsic defects and defect clusters have been ide ntified as the origin of substantial differences in properties between materials with slightly different Li/Nb ratios. While a large body of experimental literature exists, and some theoretical investigations have been performed, many of the key conjecture s in the literature, especially on defect clusters, remain unverified by either experiments or theory. The purpose of this work is to systematically analyze the intrinsic defects and defect clusters in lithium niobate using density functional theory (DFT) calculations combined with thermodynamic calculations. It appears that the growth process affects the type of defects produced. The typical composition of LiNbO3 grown from congruent melting is Li / (Li + Nb) = 0.485, which indicates a lithium deficient defect structure 10. More nearly stoichiometric compositions have been achieved t hrough vapor transport equilibration (VTE) 93, 94 and double crucible Czochralski (DCCZ) 16, 17 methods. The change of composition from 0.485 (congruent) to 0.5 (stoichiometric) causes large shifts in the Curie temperature 21, cohesive field for domain reversal 17, 21, 22, built in internal field 95, and other properties 10. In particular, it has been conjectured that the temperature stability and field dynamics of a defect cluster consisting of a niobium antisite surrounded by four lithium vacancies can explain much of the observed macroscale switching behavi or in congruent lithium niobate 10. However, no experimental verification or detailed PAGE 51 51 theoretical analysis of this intrinsic defect cluster has been presented to date. Thus, a fundamental understanding of this and other intrinsic defects in LiNbO3 is essential. Based on the experimental data, several defect models in congruent LiNbO3 have been proposed. Prokhorov et al. 23 proposed that oxygen vacancies surrounded by two lithium vacancies ( so called, Model I) dominate at room temperature. However, it was soon determined from experiments that the density of LiNbO3 increases with increasing Li2O deficiency 96, which is inconsistent with Model I. Schirmer et al. 24 concluded that niobium antisites compensated by niobium vacancies (Model II) are the dominant d efects, and that oxygen vacancies are present in negligible concentration. However, Donnerberg et al. 97, using atomic level simulations, sho wed that the formation of niobiu m vacancies to compensate t he niobium antisites is energetically less favorable than the formation of lithium vacancies This leads to Model III consisting of NbLi + 4 VLi (shown in Figure 3 1) Model III was supported by X ray and neutron diffraction studies 98 101. On the one hand, Schirmer et al. 24 point ed out that the niobium vacancy model (model I I) and lithium vacancy model (model II I) can be reconciled if it is assumed that there are ilmenite type stacking fault s in congruent LiNbO3. On the other hand, nuclear magnetic resonance (NMR) studies 102, 103 concluded that a combination of Model I and Model III provide both qualitative and quantitative agreement with Li NMR spectra Inconsistent with this is the assum ption that only Model III can be used to explain the temperature dependence of exp erimental Li and Nb NMR spectra 102, 103. These contradictory experimental results demonstrate that the nature of intrinsic defect arrangements in LiNbO3 is still not fully understood. PAGE 52 52 Figure 3 1.Schematic of possibl e defect arrangement involving one niobium antisite and four lithium vacancies (so called Model III). 3 .2 Structure of LiNbO3 and Constitute Binary Oxides Ferroelectric LiNbO3 has an R3c structure, which belongs to the t rigonal system. The lattice cons tant of LiNbO3 are a = 5.1474 and c = 13. 8561 Both lithium and niobium sit on the 6a Wyckoff position while oxygen occupies 18 b. This indicates that there are six lithium six niobium and eighteen oxygen in each unit cell. Compared with the paraelectric phase, which has a R 3c space group, the lithium and niobium ions are displaced along z direction, leading to spontaneous polarization of the system. Li2O h as an anti fluorite structure (space g roup: Fm 3m), with a lattice constant of 4.57 104. For Nb2O5, several different structures have been reported, including R Nb2O5 (C2/m), P PAGE 53 53 Nb2O5 (I4122), M Nb2O5 (I4/mmm), N Nb2O5 (C2/m), H Nb2O5 (P2/m), T Nb2O5 (Pbam), and B Nb2O5 (C2/c) 105 Among these T Nb2O5, and B N b2O5 are known to be highpressure phases 106; H Nb2O5 is observed to be stable at high temperatures (T > 1100oC) 105; R Nb2O5 is metastable and transform s to P Nb2O5 on slow heating and to N Nb2O5 on rapid heating 105. According to Holtzberg et al. 107, the room temperature phase of Nb2O5 is amorphous. Here, metastable R Nb2O5 is used to determine the chemical potential range, since it is the only phase for which complete structural information is available from the literature 105 Both lithium metal and niobium metal have Im 3m structure, whereas the molecular oxygen gas is simulated by putting an oxygen dimer in a vacuum box. The coordination numbers of the ions in the three systems can be determined from examination of the crystal structures. As shown in Table 3 1 both of the structurally distinct Nb ions in LiNbO3 structure are six fold coordinated. Likewise, the two structurally distinct Nb ions in R Nb2O5 are also six fold coordinated. Moreover, the bond lengths in the two systems determined from the DFT calculations are both close to each other and to the experimental values. The situation is different for Li2O. In particular, the Li in anti fluorite structured Li2O is 4 fold coordinated, whereas it is six fold coordinated in LiNbO3. Consistent with these differences in coordination, the calculat ed bond lengths are quite different. Table 3 1 The comparison of bond lengths and coordination numbers between Li2O, Nb2O5, and LiNbO3. Bond Length Coordination Number DFT Experiment LiNbO 3 Li O ( S ) 2.02 2.05 15 6 Li O ( L ) 2.25 2.27 15 6 Nb O ( S ) 1.90 1.88 15 6 Nb O ( L ) 2.14 2.13 15 6 Li 2 O Li O 1.98 1.98 104 4 Nb 2 O 5 Nb O ( S ) 1.95 2.02 105 6 Nb O ( L ) 2.15 2.19 105 6 PAGE 54 54 3 .3 Computational Details The Vienna Ab Initio Simulation Package (VASP) 108, 109 is employed to carry out all the calculations of crystallographic structures, elect ronic structures, and defect formation energies. The projected augmented wave (PAW) method 75, which combines much of the accuracy of an allelectron method 110 with the flexibility of the pseudopotential approach 111, is used. The Li 2 s1, Nb 4p64d45s1, and O 2s22p4 are treated as valence electrons. The cut off energy for plane wave basis set is 400 eV 74. The allowed error in energ y from relaxation is 0.001 eV. The electronic minimization algorithm for energy calculation is based on residual minimization scheme direct inversi on in the iterative subspace (RMM DIIS) 112. LDA vs. GGA : The results obtained from the local density approximation (LDA) and the general gradient approximation (GGA) for the lattice parameters of the ferroelectric phase of LiNbO3 (R3c) are compared with experimental values in Table 3 2 for a supercell containing 2 x 2 x 2 conventional unit cells (240 atoms with a total of 1440 electrons). As is typical GGA predicts larger lattice parameters than LDA. Table 3 2 Lattice parameters calculated by PAW LDA and PAW GGA, compared with previous results. USPP denotes u ltrasoft pseudopotentials; NCPP denotes norm conserving pseudopotential a () c () V ( 3 ) PAW LDA 5.049 13.686 302.198 PAW GGA 5.138 13.914 318.119 USPP LDA 74 5.064 13.667 303.523 USPP LDA 73 5.086 13.723 307.420 NCPP LDA 72 5.125 13.548 308.172 NCPP GGA 72 5.255 13.791 329.816 Experiment 5.147 15 5.151 113 13.853 15 13.876 113 317.873 15 318.844 72 As Table 3 2 shows, GGA gives excellent agreement with the experimental values, with the in plane lattice parameter, a being underestimated by 0.2% and the c lattice parameter being overestimated by 0.6%; as a result the calculated volume is within 0 .1% of the experi mental PAGE 55 55 value. The GGAPAW calculations give considerably better agreement with experimental values than the LDA and GGA calculations with the norm conserving pseudopotentials (NCPP) 72 or LDA calculations with ultrasoft pseudopotentials (USPP) 73, 74. The electronic densities of states (DOS) have been calculated using PAW LDA and PAW GGA. LDA yields a band gap of 3.42 eV, whereas the GGA value is 3.5 eV. These results are consistent with previous calculations using NCPPs of 3.48 and 3.50 eV for LDA and GGA respectively 72. The calculated band gap is also in good agreement with the value of 3.5 eV, previously calculated using the orthogonalized linear co mbination of atomic orbital method (OLCAO) 114. The band gap calculated by GGA is 8% less than the experimental value of 3.78 eV 115, 116; this agreement can be considered as excellent, since DFT is widely known to underestimate band gaps 62. The DOS calcula ted using GGA is given in Fig ure 3 2 The highest occupied valence band exhibits mainly O 2p features, whereas the lowest unoccupied conduction band mainly consists of Nb 4d electrons. This DOS is also consistent with calculations performed by OLCAO 114 and NCPP 72. Based on this analysis of structure and band structure, GGA is chosen for the calculations of DFEs. K point m esh : The M onkhorst Pack 117 method is used to carry out the integration in the Brillouin zone. Since for LiNbO3, the length of c is approximately double that of a for a 2 x 2 x 2 system, we chose to have half the number of kpoints in the z direction ( c axis) as in the x and y direction s. A 4 x 4 x 2 k point mesh has been used in the current study for both the perfect structure an d defect structure calculations. As discussed by Van de Walle and Neugebauer 118, of the DFEs because the interaction between a defect and its mirror image through periodic boundary conditions is at its maximum at PAGE 56 56 the k point mesh has been shifted from the 0.5). Figure 3 2. The density of state (DOS) and partial DOS (PDOS) using GGA for LiNbO3. System s ize : The supercell method 119121, containing a finite number of atoms and defects with periodic boundary conditions, is used to calculate the formation energies of both isolated defects and defect clusters. In general, the larger the supercell, the smaller the error due to the interaction between the defect and its mirror image, thereby more closely approximating the dilute limit that we wish to capture. Because larger supercells require a formidable amount of computational time, we ha ve systematically analyzed the system size convergence. To characterize the convergence, supercell sizes of 1 x 1 x 1, 2 x 1 x 1, 2 x 2 x 1, and 2 x 2 x 2, containing 30, 60, 120, and 240 atoms respectively, are considered. The d ifference in DFE as a funct ion of supercell size is calculated; the results of typical cases are shown in Fig ure 3 3 T he 120atom supercell and 240atom supercell produce similar results for various types of vacancies PAGE 57 57 and niobium antisite (see also Table 3 3 ). The larger size effec ts for the Nb vacancy and Nb antisite are presumably due to their larger charges and point to the dominant effects of electrostatic charge on artificial defect defect interactions Based on the energies of the two larger system size s the error in the DFE from supercell size effects is estimated to be ~ 0.1 0.2 eV. Figure 3 3. Difference in defect formation energy (relative to the 1 x 1 x 1 system) as a function of system size. Supercell sizes of 1 x 1 x 1, 2 x 1 x 1, 2 x 2 x 1, and 2 x 2 x 2 correspond to 1, 2, 4, and 8 in the xaxis For the neutral defect clusters, the electrostatic contribution to the system size dependence can be expected to be rather small and only due to dipole and higher multipole interactions ; however, the strain con tribution can be expected to be larger than for the point defects considered. As discussed below, from the analysis of the energy difference s between defect reactions the size dependence is expected to be weak enough that we can be confident in the PAGE 58 58 conclu sions drawn. Therefore, all of the analyses reported are based on results obtained from the 2 x 2 x 2 (240 atom) supercell. Table 3 3 Difference in DFE, and relative energy difference, between 2 x 2 x 1 system and 2 x 2 x 2 system for various defects L i Vacancy O Vacancy Nb Vacancy Nb Antisite Difference (eV) 0.051 0.005 0.121 0.014 Even for the relatively large system used here, the precise optimization method can have an effect on the results obtained. The approach used here follows the scheme of Van de Walle and Neugebauer; in particular, the defects are added to an optimized perfect crystal structure and then the positions of all atoms within 4.5 of the defect are optimized 118. This implies that within a system of 240 atoms approximately 116 atoms are rel axed for all the defect cases. The DFE calculated in this manner differs from that obtained from a full optimization of all the atomic positions by 0.052 eV, 0.050 eV, and 0.140 eV for the Li vacancy, O vacancy and Nb antisite respectively. A similar approximation has been used for defects in GaN 118 and ZnO system 122. An alternative to using large system sizes to minimize the effects of t he artificial interactions between defects and their images would be to compensate for them in an analytic manner. In particular, Makov and Payne 123 have developed a general approach to correct the calculati on for cubic system. However, this approach has been found to lead to an overestimate of the correction term for some semiconductor systems 118, 124, 125. We have not included a Makov Payne correction in our charged supercell analysis, but rather rely on our s ystem size being sufficiently large, an assumption that we justify by the data in Fig ure 3 3 3 .4 Thermodynamic Framework T he formation energy of a defect or a defect cluster depends on the chemical potential of the ions added or removed from the perfect crystal to form the defect. Below we outline a PAGE 59 59 thermodynamically consistent process to determine the physically possible range of chemical potentials This was first developed in the context of binary oxides 122 126. The extension to ternary oxides used here is substantially the same as that used previously for LiNbO3 74, BaTiO3 127 and SrTiO3 128. The stability of the system against decomposition into its constituent elements places upper limits on the chemical potential of each element in LiNbO3. ] [ ] [3Metal LiNbOLi Li [3 1 ] ] [ ] [3Metal LiNbONb Nb [3 2 ] ] [ ] [2 3O LiNbOO O [3 3 ] The total energy of a stoichiometric unit of LiNbO3 can be expressed as: ] [ ] [ 3 ] [ ] [3 3 3 3LiNbO E LiNbO LiNbO LiNbOtotal O Nb Li [3 4 ] Thus lower limits of chemical potential for each element can then be obtained: ]) [ ] [ ] [ ( 3 1] [3 3Metal Metal LiNbO E LiNbONb Li total O [3 5 ] ] [ ] [ 3 ] [ ] [2 3 3Metal O LiNbO E LiNbONb O total Li [3 6 ] ] [ ] [ 3 ] [ ] [2 3 3Metal O LiNbO E LiNbOLi O total Nb [3 7 ] Furthermore, since neither Li2O nor Nb2O5 precipitates from bulk LiNbO3, the range s of chemical potential s are subject to the additional constraints: ] [ ] [ ] [ 22 3 3O Li E LiNbO LiNbOtotal O Li [3 8 ] ] [ ] [ 5 ] [ 25 2 3 3O Nb E LiNbO LiNbOtotal O Nb [3 9 ] The above sets of equations define the ranges of chemical potential consistent with the stability of LiNbO3 against decomposition into binary oxides or into its elemental components. PAGE 60 60 Knowing the boundary of chemical potent ial for each element, the range of possible values for the DFEs of single point defects can be calculated. To characterize the defect reaction, knowledge of the range of chemical potentials is not sufficient: the chemical potential of each element has to b e known. Physically, it is possible to change the oxygen chemical potential through the temperature and/ or oxygen partial pressure. There are two limits on the oxygen chemical potential in LiNbO3. Under extremely oxidizing conditions, the chemical potenti al of oxygen reaches its maximum. The criterion for the equilibrium of oxygen in LiNbO3 and O2 is: ] [ ] [2 3O LiNbOO O [3 10] Under reducing conditions, the Li and Nb in LiNbO3 must be in equilibrium with metallic Li and Nb. This places a lower limit on the oxygen chemical potential of: ]) [ ] [ ] [ ( 3 1] [3 3Metal Metal LiNbO E LiNbONb Li total O [3 11] As we shall see, these constraints actually correspond to unrealistically high and low oxygen partial pressures respectively, with the physically attainable r ange of partial pressures being much narrower. Bringing all of these stability criteria together, we can construct the ternary chemical potential map shown in Fig ure 3 4 where the points in chemical potential space should be analyzed in a manner analogou s to that for the composition in a conventional ternary phase diagram 129. In the following, the points and lines refer to the ternary chemical potential map in Fig ure 3 4 Case 1: Li2O r eference s tate ( l ine AD) : When Li2O is chosen as a reference state, the compositionweighted sum of the lithium and oxygen chemical potential s is equal to the total energy of Li2O. Likewise the criterion for the equilibrium of Li in LiNbO3 and Li2O is: PAGE 61 61 ) ] [ ( 2 1 ] [2 32O total O Li LiO Li E LiNbO [3 12] The chemical potential of Nb in LiNbO3 can then be determined: ] [ 3 ] [ ] [3 3 32 2LiNbO LiNbO E LiNbOO Li Li O total O Li Nb [3 13] Case2: Nb2O5 r eference s tate (l ine BC) : In this case, the Nb in LiNbO3 and Nb2O5 are in equilibrium: ) 5 ] [ ( 2 1 ] [5 2 35 2O total O Nb NbO Nb E LiNbO [3 14] The chemical potentials of Li in LiNbO3 can then be determined as: ] [ 3 ] [ ] [3 3 35 2 5 2LiNbO LiNbO E LiNbOO Nb Nb O total O Nb Li [3 15] Reconciliation of the t wo r eference s tates : One might think that these separate analyses would be equivalent. To explore this, we can take Li [LiNbO3] in Equation 3 12 and 3 15 in term s of individual energy contributions. From Equation 3 12 we get the chemical potential for Li in LiNbO3, with a Li2O reference as: Li Li2O [LiNbO3] as: ) ( ] [ ] [2 2 1 32O Li H Metal LiNbOf Li O Li Li [3 16] Likewise from Equation 3 15, we can determine the chemical potential of Li in LiNbO3 with an Nb2O5 reference as: ) ( ) ( ] [ ] [5 2 2 1 3 35 2O Nb H LiNbO H Metal LiNbOf f Li O Nb Li [3 17] Thus: diff O Li Li O Nb LiH LiNbO LiNbO ] [ ] [3 32 5 2 [3 18] where )) ( ) ( ( ) (2 5 2 2 1 3O Li H O Nb H LiNbO H Hf f f diff [3 19] PAGE 62 62 As Table 3 4 shows, both experimental data and our GGA calculations show that Hdiff is non zero. As a result, the two reference states are thermodynamically different from each other. Table 3 4 The heats of formation of Li2O, Nb2O5, and LiNbO3 calculated from the GGA calculations, compared with experimental results. Both experimental and calculated values are at T=0 K GGA(eV) Experiment(eV) H(Li 2 O) 6.280 6.139 130, 131 H(Nb 2 O 5 ) 18.262 19.202 130 19.943 130 19.687 130 19.735 130 19.666 130 19.813 130 19.687 131 H (LiNbO 3 ) 14.433 14.149 132 H diff 2.162 1.236 The result of this analysis is that the chemical potentials of all three elements are constrained. The bounds that these constraints impose could be determined from experiments on the range of stability of the various oxides ; unfortunately the complete experimental data se t required is not available We therefore chose to use the results of GGA calculations to determine these constraints, as we are then able to develop a completely consistent data set of chemical potentials and DFEs Calculation of c hemical p otentials: a ppl ication to LiNbO3: Using the method outlined above, the GGA calculations allow us to explicitly define the domain of possible chemical potentials. To accurately calculate the ranges of chemical potential for each element the energy of lithium metal, niobium metal, oxygen gas, Li2O, Nb2O5 and LiNbO3 have been determined using GGA (Table 3 5 ). The heats of formation for Li2O, Nb2O5 and LiNbO3 have been calculated from the difference between energies of the oxides and compositionweighted sum of the energies of the elements. From the data shown in Table 3 4 the error associated with experimental results for PAGE 63 63 Nb2O5 is significantly larger than for Li2O and LiNbO3. In particular, the calculated Hf values for Li2O and LiNbO3 are within 2.3 % of the experimental values (Table 3 4 ), whereas for Nb2O5 the error is between 4 and 10 % of the experimental values. This larger deviation in f for Nb2O5 is not surprising due to the experimental uncertainties in the structure discussed above U sing the val ues shown in Table 3 5 and Equation 3 1 to Equation 3 10, the range s of possible chemical potential of each element are calculated. We then consider the chemical potential of oxygen. The dependence of O on temperature and partial pressure can be calculated if the dependence of O [T, Po] at one particular pressure Po is known 126: ) ln( 2 1 ] [ ] [o o o oP P kT P T P T [3 20] Table 3 5 The energy of lithium metal, niobium metal, oxygen gas Li2O, Nb2O5, and LiNbO3 as calculated by GGA Space Group Energy (eV) per Atom Li (Metal) Im 3m 1.895 Nb (Metal) Im 3m 10.049 O 2 (Gas) 4.392 Space Group Energy (eV) per Formula Unit LiNbO 3 R3c 39.552 Li 2 O Fm 3m 14.461 R Nb 2 O 5 C2/m 60.319 Thus high values of O correspond to high oxygen partial pressures and conditions that tend to favor oxidation. Similarly, low values of O correspond to low partial pressures of oxygen and conditions that favor reduction. The value of O calculated by DFT is at 0 K. Its depend ence on temperature can be derived from experimental thermodynamic data, yielding: ) ( 2 1 ] 0 [ ] [2 ] [2O P T G P K P To o O o o [3 2 1 ] The experimentally determined temperature dependence of the oxygen chemical potential at 1 atm is given in Table 3 6 130. The requirements of stability of LiNbO3 d iscussed above, place limits o n the oxygen partial pressures which, at 298 K, can lie in range of 4.910139 to 2.61018 PAGE 64 64 atm. Of course, both of these limits are physically unreasonable. Thus in Fig ure 3 4 the more physically reasonable pressure range of 1020 atm to 10 atm is also indicated. Table 3 6 Change in oxygen chemical potential with respect to 0 K value from experiment 130. Temperature (K) Change in u O (eV) 100 0.150 200 0.341 298.15 0.544 300 0.548 400 0.765 500 0.991 600 1.233 700 1.460 800 1.702 900 1.949 1000 2.199 Referring to the ternary stability map in Fig ure 3 4 at any point within the triangle, the sum of the chemical potentials of each element yields the total energy of LiNbO3. Furthermore, considering the constraints imposed by Equation 3 8 and Equation 3 9 we can further narrow down the range of chemical potentials. The region enclosed by points A, D, and PLi Nb sati sfies Equation 3 8 while the region enclosed by points B, C, PO Li, and PO Nb satisfies Equation 3 9 The in tersection of these two regions defines the thermodynamically allowable range of chemical potentials ; t his stability region is thus defined by the shaded quadrilateral enclosed by the points A, B, C, and D. Therefore, to study the general case, the DFE calculations for all points within the stability range are performed. At this point, all of the methods needed for defect energy calculations are defined. In the next Section, we present the results. 3.5 Defect Formation Energy T he current paper focuses on following intrinsic defects: Li Frenkel, Nb Frenkel, O anti Frenkel, Li2O pseudoSchottky, Nb2O5 pseudoSchottky, LiNbO3 Schottky, and three defect PAGE 65 65 clusters that have been discussed in the literature: 5 NbLi +4 VNb 2 NbLi +3 VLi + VNb and NbLi + 4 VLi 10. Figure 3 4. Stability range of chemica l potentials (in eV) of the elements in LiNbO3. The region enclosed by points A, D, and PLi Nb satisfies Equation 3 8 while the region enclosed by points B, C, PO Li, and P O Nb satisfies Equation 3 9 The intersection of these two regions defines the the rmodynamically allowable range of chemical potentials. This stability region is thus defined by the shaded quadrilateral enclosed by the points A, B, C, and D. Line AD represents using Li2O as reference state; line BC represents using Nb2O5 as reference. T he oxygen partial pressure range (in atm) is that for room temperature The chemical potential unit is in eV/atom. For computational efficiency, the DFEs of these intrinsic defect complexes are first determined from the formation energies of the individual charged defects that make up the complex. Thus, the formation energy of Li2O pseudo Schottky comes from separate calculations on the forma tion energy of lithium vacancy and oxygen vacancy. The association energy resulting from the interaction of the indiv idual charge d defects in these neutral defect clusters will be discussed in section 3.6 PAGE 66 66 Before characterizing the defect clusters, it is important to analyze the individual defects with various charge states. The following section will discuss all the res ults obtained for individual defects followed by the neutral charged clusters. 3 .5.1 Forma lism The DFE of a defect or defects, denoted as with charge state q is defined as 118 ) ( ) ( ) ( ) ( V E q n perfect E q E P T q Ev F i i i total total f [3 22] 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 t he bulk single crystal system. Therefore, as discussed by Van de Walle and Neugebauer 118, 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. Strictly speaking, the free energy rather than the tot al energy should be used in Equation 3 22 for the calculation of the D FE s The total internal energies of a supercell obtained from DFT calculations correspond to the Helmholtz free energy at zero temperature, neglecting zero point vibrations 126. 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 118. As a consequence, we do not expect this neglect of the entropy term to qualitatively change our PAGE 67 67 conclusions. Detailed analyses by He et al. 133 and Kohan et al. 122 also concluded that entropic effects can be neglec ted. 3.5.2 Point Defects T he DFEs of all the native point defects with various charge states in LiNbO3 have been calculated using Equation 3 22. The interstitial sites are assumed to be at (0.0, 0.0, 0.13936) for individual defects, which correspond s to the center of an empty oxygen octahedron. Since the DFE depends on the location of the Fermi energy, the influence of the Fermi energy on the stability of each individual defect is considered In the current study, the reference zero of Fermi energy is as signed to be the valence band maximum (VBM) of the perfect structure ; the highest Fermi energy thus corresponds to the conduction band minimum (CBM). As discussed by Van de Walle et al. 118, it is necessary to implement an additional procedure to calculate the shift of bands caus ed by the presence of defects. Thus, to calculate this shift, t he average electrostatic p otential difference ( 3 22) between the perfect structure and t he defected structure has been determined using the approach of Kohan et al. 122. For pure stoichiometric LiNbO3, which is an insulator, the Fermi energy is expected to be at the middle of the band gap. In the present paper, we are simulating the dilute limi t of pure stoichiometric system; therefore the shift of Fermi energy is considered negligible. The dependence of DFEs on the Fermi lev el is shown in Fig ure 3 5 using Nb2O5 (the BC line in Fig ure 3 4 ) as the reference state. For each individual defect, Fig ure 3 5 includes only th e charge states with lowest DFE values over the corresponding Fermi level range. The slopes of the lines in the Fig ure 3 5 represent the charge states of the defects. When Fermi energy is at mid dle of the band gap, VLi VNb VO .., NbLi ., Lii ., Nbi ., and Oi .. a re the energetically favored charge state for each defect. It is also observed that, as the Fermi energy increases (i.e. moves from the PAGE 68 68 VBM to wards the C BM), the thermodynamically stable charge of the positive defects decreases. Correspondin gly, the charge states of the negative defects decrease as the Fermi energy increases. This order of stability is broadly consistent with that determined in the DFT calculations of Li et al.74 Figure 3 5. Defect formation energies of various point defects as a function of Fermi energy; the lowest energy charge state is given in each case. The Fermi energy ranges from EF=0 (left) at the VBM to EF=3.5eV at the conduction band minimum (right) In Fig ure 3 5 the points at which the slopes change are the thermodynamic transition q1 / q2). At these transition points the Fermi level of the charge state q1 and q2 have the same energy. For Nb related defects, the stable Nb antisite changes from 3+ to 1+ directly without passing through the 2+ charge state. C orrespondingly, the Nb interstitial changes from 2+ to 0, and Nb vacancy changes from 2 to 4 . To characterize this phenomenon, the thermodynamic transition level of adjacent charge state of Nbrelated and O related defects are calculated and shown in Fig ure 3 6 PAGE 69 69 Figure 3 6. Thermodynamic transition levels of O related and Nb related defects. The values under the line are with respect the valence band maximum (VBM). The transitions enclosed by squares have a negative U character. It can be seen that there are order switching for the Nb antisite, Nb interstitial and Nb / 4 ) is lower than (2 / 3 ). These switchings of the order of thermodynamic transition levels are indicative of the so called negative U effects 134. However, as Fig ure 3 6 shows, the energies for the thermodynamic transition levels for the respective Nb related defects discussed above are extremely close to each other; indeed they are close to the estimated uncertainties in the present calculations. The electron lattice interaction result in different atomic relaxations for different charge state of the defects. Taking the Nb antisite as an example, the average bond length of 0, 1+, 2+, PAGE 70 70 3+ and 4+ are 2.130 2.092 2.063 2.038 and 2.011 respectively. These decrea ses in bond lengths as the charge state of Nb antisite increases are due to the increasing electrostatic interactions. 3 .5.3 Defect Clusters The constraint of charge neutrality means that only neutral defect complexes will be present in the system. Wh ile this leads to a need to consider a plethora of different defect clusters, it does provide the simplification of not having to consider the effects of the Fermi energy, since this term drops out in Equation 3 22. In the current section, the charged defe ct that make the neutral defect are assumed to be infinitely far away from each other. The as sociation effects are added in section 3.6 From Fig ure 3 7 it is can be seen that formation energy of the Li Frenkel, O anti Frenkel, Nb Frenkel, and LiNbO3 Sch ottky are independent of the choice of the chemical potential. This is because these defects preserve the stoichiometry of the system, obviating any chemical potential term. By contrast, the Li2O pseudoSchottky, Nb2O5 pseudoSchottky, 5 NbLi +4 VNb 2 NbLi +3 VLi + VNb and NbLi + 4 VLi defect clusters change the stoichiometry of the system. Thus, their formation energies are influenced by the reference state of chemical potential. The calculations show that the formation energy for Nbantis ite compensated purely by Nb vacancies (5 NbLi +4 VNb ) is different for the two reference states ; however it is energetically unfavorable under both conditions. The DFE for NbLi + 4 VLi clusters is positive with the Li2O (AD in Fig ure 3 4 ) reference state but negative for the Nb2O5 (BC in Fig ure 3 4 ) reference state, indicating that this defect cluster will occur spontaneously under these conditions PAGE 71 71 Figure 3 7. Formation energies of neutral defects and defect clusters under Nb2O5 rich conditions (blue) and Li2O rich conditions (red) The Frenkel defects analyzed above consist of individual fully charged defects. Frenkel defects with partial charged or neutral individual defects can also be considered. However, the formation energies o f Frenkel defects with partial charged or neutral individual defects all turn out to have higher energies than those formed from normal fully charged individual defects. For example, the formation energy of the Frenkel defects consisting of Li i x+ VLi x is 1.9 eV higher than the Frenkel defect formed from Li i .+ VLi Table 3 7 compares the results of our DFT calculations with the previous DFT calculations and with calculations using empirical potentials. For the purposes of comparison, it is instructive t o distinguish between those neutral defects that maintain the systems stoichiometry (Li Frenkel, Nb Frenkel, O Frenkel, and LiNbO3Schottkty) and those defects and defect clusters that do not ma intain the system stoichiometry. For the stoichiometric defects, all four calculations conclude PAGE 72 72 that the Li Frenkel defect has the lowest formation energy. They also all predict that the Nb Frenkel has the highest DFE. Three of the four calculations conclude that the LiNbO3 Schottky has a lower DFE than the O Frenkel, with only the older atomistic calculations of Donnerberg et al. 97 giving the opposite order. Table 3 7 Defect formation energies for various defects compared with pr evious electronic structure and atomi sti c simulations Defects (eV/defect) Jackson (Atomistic) Donnerberg (Atomistic) USPP LDA PAW GGA PAW GGA Reference State None None Nb 2 O 5 Nb 2 O 5 Li 2 O Li Frenkel 1.37 0.93 1.78 0.77 0.77 Nb Frenkel 11.72 6.26 6.43 3.80 3.80 O Frenkel 4.76 3.42 5.80 3.10 3.10 LiNbO 3 Schottky 3.95 3.91 4.84 2.25 2.25 Li 2 O Pseudo Schottky 1.81 1.94 0.89 0.07 1.51 Nb 2 O 5 Pseudo Schottky 5.09 2.85 7.29 3.19 2.57 Nb Li + 4 V Li 1.51 2.74 4.40 0.98 1.62 Nb Li +4 V Nb 10.23 5.06 2.08 2.03 3.47 With respect to the nonstoichiometric defects, it is difficult to compare the results from the atomistic simulations with the DFT calculations, since the atomistic calculations neither include an explicit reference nor take into acco unt the chemical potential. However, the results of the DFT calculations 74 can be compared, despite the concerns with the previous DFT calculations mentioned in the Introduction. In particular the USPP and PAW calculations yield the same order of stability of defec ts for the Nb2O5 reference state, with the NbLi + 4 VLi cluster having the lowest energy which, importantly, both analyses predict to be negative. This is strong evidence that this defect cluster should form spontaneously in Nb2O5 rich environments. Although they agree in the order of stability of the defects, the corresponding USPP and PAW calculations yield considerably different specific values for DFEs. Unfortunately the paper of Li et al. 74. does not break down the final DFEs into the contributions from the presence of the PAGE 73 73 defect and the contributions from the chemical potentials; it is thus not possible to analyze these differences further. 3 .5.4 Defect Stability Range The DFEs calculated so far using DFT are based on specific chemical potentials Knowledge of the defect formation dependence over the complete chemical potential stability range is also of considerable importance. The lowest energy defect as a function of chemical potential of each element is given in Fig ure 3 8 It is observed that Li Frenkel and NbLi + 4 VLi are th e two important defects The Li Frenkel is the dominant defect within region AEFD (in Fig ure 3 4), whereas the NbLi + 4 VLi cluster is the dominant defects within region EBCF (in Fig ure 3 4 ). F igure 3 8. Defect stability range: in the region AEFD, the Li Frenkel is the dominant defect reaction; whereas in EBCF, the niobium antisite compensated by lithium vacancies is dominant T he Li Frenkel is the dominant defect throughout the whole range of O along the AD line, while the defect cluster consists of niobium antisite compensated and four lithium vacancies cluster dominates along BC f or all oxygen partial pressures. PAGE 74 74 As discussed above the dominant defects might change as O transitions from its maximum (oxidizing conditions) to its minimum (reducing conditions). Remarkably, however, the Li Frenkel is the dominant defect throughout the whole range of O along the AD line, while the NbLi + 4 VLi cluster dominates along BC for all oxygen partial pressures. 3 .6 Association Effects and Defect Configurations The above defect reaction analysi s is based on the formation energies of individual point defect ; that is, the defects are assumed to be so far away from each other that there i s no interaction between them. Here, the validity of this approximation is examined by calculating the formation energies of the clustered defects. The change in energy relative to the well separated defects is the association energy: negative association energy corresponds to the preference for the defects to cluster. Calculations of such associated clusters require the choice of specific defect arrangements. For the Frenkel defects, the vacancy and interstitial are assumed to be first neighbors. For the L i2O pseudoSchottky, an oxygen atom and two lithium atoms in its first nearest neighbor have been removed from the system. For the Nb2O5pseudoSchottky, two nearest Nb atoms are removed along with five oxygen atoms that are the first nearest neighbor to t he removed Nb atoms. The arrangement for the NbLi + 4 VLi follows the model of Kim et al. 135, which assumes the defect complexes composed of a niobium antisite surrounded by three Li vacancies in the neares t neighbor plus one independent vacancy along the z direction. We analyze different possible arrangements for these defects. For 5 NbLi +4 VNb initially a Li site is chosen and replaced with an Nb atom. Then four of the first nearest Li atom sites are replaced with Nb to obtain the 5 NbLi configuration. F our Nb sites which are the nearest neighbor of the first Li atom replaced are then removed from the system to create 4 VNb PAGE 75 75 In all cases shown in Table 3 8 the defects have a tendency to cluster, as evidenced by the negative association energies. Since, this association does not change any of t he thermodynamics of the system; the se diffe rences in energy are independent of the reference state. Interestingly, although association lowers the energy of all of the defect clusters analyzed, it does not change their relative order. This gives us high confidence that the conclusions drawn from the isolated defect calculations are valid. The simulations of Li, Nb and O Frenkel pairs resulted in the perfect structure with no defects, since the vacancies and interstitials of the corresponding pairs were too close to each other. Table 3 8 Associatio n energy of Li2O Pseudo Schottky, Nb2O5 PseudoSchottky, 5 NbLi +4 VNb and NbLi + 4 VLi Association Energy Li 2 O Schottky 0.23 Nb 2 O 5 Schottky 1.20 Nb Li + 4 V Li 0.25 Nb Li +4 V Nb 1.12 Previously, the DFE of the NbLi + 4 VLi cluster without considering association effects was determined from the sum of the DFEs of its isolated constituent defects 136. That is, each defect was assumed to be at an extremely long distance from others. However, there may be a significant difference between the energy of isolated defects and any specific structural variant of the cluster arising from the electrostatic and elastic interaction among the defects. To explore the association of the lithium vacancies around the niobium antisite, several carefully chosen arrangements are considered. Since both the Nb antisite and the lithium vacancies lie on the Li sub lattice, the specific structure can be characterized in terms of this sub lattice alone; the energetics are governed by the structural relaxation of all the ions in the clu ster and in its vicinity. Because there are so many possibilities, it is not feasible to characterize all the possible arrangements of the Li vacancies. However, the most important distinction is between those in PAGE 76 76 which the vacancies are very close to the N b antisite and those in which defects are a considerable distance apart. Thus, the clusters can be categorized in terms of the number of first nearest neighbors (FNNs), with the remaining vacancies further away. (i ) 0 lithium vacancy as FNN and 4 far away (A0) (ii) 1 lithium vacancy as FNN and 3 far away (A1) (iii) 2 lithium vacancies as FNNs and 2 far away (A2) (iv ) 3 lithium vacancies as FNNs and 1 far way (A3) (v ) 4 lithium vacancies as FNNs (A4) In the A0 structure, there is no association of the vacancies; this is the case analyzed previously and found to have negative formation energy for the Nb2O5 reference state 136. For A1 and A2, only the lithium vacancies in the first nearest neighbor positions are explicitly included in the supercell calculation. Thus, taking A2 as an example, two lithium vacancies are included in the supercell; the DFE of the other two lithium vacancies are assumed to be the same as that of an isolated single lithium vacancy. The DFE of isolated lithium vacancies are 1.13 eV and 3.29 eV under Nb2O5 and Li2O reference states, respectively. The effects of Fermi level on the DFE of lithium vacancy are given by Xu et al 136. T he total formation energy of the entire defect cluster is the sum of the DFEs the partial cluster and any isolated Li vacancies. For A3 and A4, to maintain the overall c harge neutrality and characterize possible local charge effects, all the lithium vacancies are explicitly included in the supercell calculation. Thus, taking A3 as an example, three lithium vacancies are in the FNN position while the fourth lithium vacancy is far away from defect cluster. The dependence of the DFE on the number of lithium vacancies in the FNNs of niobium antisite is given in Fig ure 3 9. PAGE 77 77 As illustrated in the figure, the energy decreases with the number of FNNs. Thus, A1 and A2 cases are not considered further. The discussion is focused on the A3 and A4 cases. Since there are six possible sites for the vacancies in the first neighbor shell, there are a number of structurally distinct variants of the A3 and A4 clusters. Each data point in Fig ure 3 9 corresponds Figure 3 9. Defect formation energies dependence on the number of lithium vacancies in the first nearest neighbor positions of niobium antisite. For A3 and A4, different points represent different configurations, as discussed in Section.3.6 to a distinct configuration. Indeed, for the energy of the A4 the association energy is ~ 0.25 eV/defect, which characterizes the attractive interaction among these charged defects. Because the association energies of the A3 and A4 structure s are very similar, the energetics of a number of variants of each is examined in detail below. Configuration A3: While there are a finite number of crystallographically distinct arrangements PAGE 78 78 of the three FNN vacancies, the fourth far away vacancy could be in an almost unlimited number of distinct locations. One particular A3 model was proposed by Kim et al. 135, Figure 3 10. Local structure of niobium antisite. A. Kims model : Defect cluster 4 5 6 8 (red) is parallel to the bulk polarization direction; cluster 1 2 3 7 (black) is anti parallel to the bulk polarization direction. B. Proposed configuration Defect clusters 12 3 8 and 1 2 3 9 (red) are parallel to the bulk polariz ation direction with different value of dipole moment. Defect cluster 4 5 6 7 (black) is anti parallel to the bulk polarization direction. The bulk polarization direction is shown in blue. The oxygen sublattice is not shown for clarity. C. Configurations of A4. 1 6 represent the six possible positions of lithium vacancies in the first near neighbors (FNNs). The dipole moment direction is from negative charged defect center (lithium vacancies) to positive charged defect center (niobium antisite). The bulk po larization is set to in the figure. in which the fourth Li vacancy lies along the z direction on the opposite side of the niobium antisite, as shown in Fig ure 3 1 0 (a). Two different configurations using Kims model are considered in the present study: one in which the dipole moment of the cluster is parallel to the ferroelectric polarization (4 5 6 8 Kim parallel ), the other with the defect dipole is anti parallel to the ferroelectric polarization (1 2 3 7 Kim antiparallel ). The polarization of each syst em is PAGE 79 79 shown in Table 3 9 the Kim parallel structure in which the defect and bulk polarization are parallel has a slightly lower energy (by 0.03 eV) than Kim antiparallel structure in which defect and bulk polarization are anti parallel. Because both of th e structures preserve the rhombohedra symmetry of the system, there are no components of polarization orthogonal to the bulk uniaxial polarization; as a result the polarization remains purely uniaxial. In Fig ure 3 1 0 (b) three configurations have been proposed, which could yield much large r electric dipole mo ments than Kims model (Table 3 9 ). These correspond to vacancies on the 1 2 3 8, 12 3 9 and 4 5 6 7 sites, respectively. Case s 1 2 3 8 and 1 2 3 9 have dipole moments parallel to the bulk polarization direction while the dipole moment of 4 5 6 7 is antiparallel to the bulk polarization. These configurations are similar to Kims model, in that three lithium vacancies lie in an FNN plane, (1, 2, and 3 in Fig ure 3 1 0 (c) ). However, they differ in the position of the fourth lithium vacancy (8 or 9), which in all cases lies on the same side of the Nb antisite as the plane of three Li FNN vacancies. The polarization and DFEs calculated based on these proposed configurations are given in Table 3 9 Table 3 9 Defect formation energies and polarization of various A3 configurations. Kim parallel represents that both defect and bulk have the same direction of polarization. Kim antiparallel represents that the polarization direction of defect and bulk is opposite. The 1 2 3 8 and 4 5 6 7 configurations are illustrated in Fig ure 3 1 0 (a) and Fig ure 3 1 0 (b). All these models only possess polarization along [0001] direction. The bulk polarization is calculated to be 63.7 C/cm2, which is in excellent agr eement with experimental value 62 4 C/cm2 base on near stoichiometric sample 137 The polarization of congruent LiNbO3 is ~ 71 C/cm2 138. Defect Cluster DFE (eV /Defect) System Polarization ( C/cm2) Polarization Change Caused by Defect Clusters( C/cm2) 4 5 6 8 (Kim parallel) 135 1.222 60.5 3.2 1 2 3 7 (Kim antiparallel) 135 1.191 59.8 3.9 1 2 3 8 (parallel) 1.224 69.2 5.5 1 2 3 9 (parallel) 1.201 74.5 10.8 4 5 6 7 (antiparallel) 1.223 60.6 3.1 PAGE 80 80 The comparison of the DFEs and polarizations of both Kims model and the variants of the proposed configuration are given in Table 3 9 There is only ~0.03 eV per defect difference between the highest and lowest energy variants. According to Gopalan et al. 10, when the dipole moment of the defect and polarization of the bulk are in the same direction, the system has a lower energy. Gopala n et al. 10 considered this to be the equilibrated state. When the dipole moment of the defects is in the opposite dir ection from the bulk ferroelectric polarization, the polarization of the defected system is lower than the bulk value and the energy is higher. This may be considered to be a frustrated state. Table 3 9 does not show any clear relationship between the pola rization and DFE. For the corresponding structures, such as parallel vs. anti parallel Kims model, th e DFE decreases as the local polarization of the system (ferroelectric plus defect dipole) increases, which is consistent with the prediction of Gopalan e t al.10 However, it is not possible to compare inequivalent structures because the difference can be caused by both polarization and structure. Configuration A4: Taking four lithium atoms from six possible FNN positions yields fifteen (C6 4 = 15) possible defect clusters. However, due to the three fold symmetry of LiNbO3 basal plane, only four of these configurations are distinct; these are listed in Table 3 10 and illustrated in Fig ure 3 10 (c). All of the A4 structures have lower energies than the lowest energy A3 structure. The energy range among these four structures is only 0.047 eV, which is equivalent to the therma l energy at 545 K. This indicates that all such configurations are likely to be present at room temperature. All of these defect arrangements break the symmetry of the ferroelectric phase, thereby leading to non uniaxial contributions to the polarization ( see Table 3 10). These nonuniaxial components are not negligible, amounting to up to 7.9% of the uniaxial component (for the A4 A structure). However, since all the symmetry equivalent structures of PAGE 81 81 each of the four structures are equally likely, the syst em maintains the overall, though not local, polarization of a bulk system as uniaxial, as the nonuniaxial contributions of an ensemble of defect clusters can be expected to cancel out. Table 3 10. Defect formation energies and polarization of various A4 configurations (see Fig ure 3 10 (c ). D* represents degeneracy. Change in Polarization represents that polarization changes caused by the defect cluster in [0001] direction using bulk as reference. Defect Cluster Vacancy positions D DFE (eV /Defect) System Polarization (C/cm 2 ) [2 1 1 0] [ 121 0] [00 0 1] Change in Polarization (C/cm2) A4 A 1 3 5 6 6 1.296 3.5 3.4 61. 6 2.1 A4 B 1 2 3 5 3 1.254 1.3 2.3 58.2 5.5 A4 C 3 4 5 6 3 1.294 0.3 0.0 61. 6 2.1 A4 D 1 2 5 6 3 1.301 2.3 1.1 64. 5 0.8 A ll of these configurations considered have a negative DFE, indicating spontaneous formation of the defects. However, the current study is based on the assum ption of a dilute solution of defects; presumably the formation energies will increase with increasing concentration, such that the system remains stable. A previous study and Li et al. 74 also observed negative formation energies for this material. Furthermore, the energy difference between different configuratio ns is rather small. This indicates that all the configurations could exist in the system. Therefore, the polarization of the system, such as congruent LiNbO3, is an overall effect of different configurations. The polarization of the 1 2 3 8 system is very similar to the experimental values for the polarization. However, the comparison between simulation and experiments should be treated with caution because the dipole moments calculated in the present study represent an average over the simulation supercell T hus it represents a local polarization rather than the experimental polarization over a macroscopic region. PAGE 82 82 3 .7 Diffusion Mechanisms As the above analysis has shown, the DFEs of all of the A3 and A4 complexes are within 0.11 eV of each other. It is therefore of considerable interest to determine whether these structure s can be expected to be essentially static or whether the energy barr iers to Livacancy migration are sufficiently low that they might be highly dynamic with Li vacancies diffusing among the various nearly energetically equivalent sites. If the defects are indeed dynamic, the non uniaxial components of the polarization of e ach individual cluster may dynamically cancel out locally. To set the context for this analysis, the diffusion of an isolated Li vacancy through the Li sublattice in LiNbO3 is first considered. Based on a qualitative structural analysis, the energy barriers of s everal diffusion paths have been calculated ( Fig ure 3 1 1 ). Path A: D iffusion of a lithium vacancy directly to its first nearest neighbor in the Lisublattice Each Li ion has six first Li neighbors at a distance of 3.771 Therefore, the jump direction for an individual Li ion can be any one of the following six directions: [0 2 2 1], [ 2 0 2 1], [2 2 0 1], [0 2 2 1], [2 0 2 1], and [ 2 2 0 1]. As the diffusion directions maintain the three fold symmetry, for both the plane above and be low, the overall effect would lead to an isotropic diffusion mechanism It is noted that such diffusion is not in the basal (0001) plane. Path B: D iffusion of a Li through the neighboring vacant site then to its first nearest neighbor in the Li sublattic e. For example, in Fig ure 3 11(b), lithium vacancy first move towards [1 1 0 0] then to [0 0 0 1] to reach [2 2 0 1]. Path C: D iffusion within the basal plane. The diffusion distance for lithium vacancy within the basal plane (0001) is 5.138 along [ 2110] The coordination number of a vacancy is also six. However, this diffusion direction is perpendicular to the [0001] direction, which leads to anisotropic diffusion mechanism. PAGE 83 83 Path D: Along the [0001] direction, the stacking of cations follows the order: Li Nb vacant site Li Nb vacant site; thus, diffusion of a lithium vacancy directly along the z direction will not be favorable as the niobium ions block the diffusion path. This path is not considered further. (a ) (b) Figure 3 1 1 (a ) Schematics of p ossible diffusion paths in the LiNbO3, vertical direction is [ 0001] Diffusion path A represents the diffusion of lithium vacancy directly to its first nearest neighbor in the Li sublattice Path B represents the diffusion through the va cant site then to its first nearest neighbor in the Li sublattice Path C represents the diffusion within the basal plane. Path D is the diffusion of lithium vacancy directly along z direction, which is blocked by the niobium ions. The energy favorable pat h is highlighted with red and in solid line. (b) The diffusion direction of lithium vacancy to its six first nearest neighbor shown in three dimensions. This corresponds to the diffusion path A in the figure left (black).The diffusion path B is shown in bl ue. Only Li sublattice is shown. The energy barriers associated with each of path A, B and C are given in Fig ure 3 1 1 It is found that the migration energy for diffusion within the basal plane (Path C) is much higher than the values for diffusion between the first nearest neighbors of the Li sublattice (Paths A and B). It is also observed that the diffusion path between the first nearest neighbors in the Li sublattice PAGE 84 84 is neither in a straight line nor exactly through the vacant site, but is a compromise between the two. Furthermore, the diffusion barriers for this path is calculated to be 1.63 eV (Figure 3 12) which is in excellent agreement with the experimental values of the Li migration energy of 1.55 eV 139, 1.62 eV 140, and 1.17 eV 141 (Table 3 11). This result is c onsistent with the conjecture that the diffusion mechanism in LiNbO3 is through the migration of lithium vacancies on the Li sublattice 142. Figure 3 1 2 NEB calculations of diffusion barriers for lithium vacancy for s everal diffusion paths. Path A: FNN Direct represents diffusion of lithium vacancy directly to its first nearest neighbor in the Li sublattice Path B: FNN curve represents the diffusion through the vacant site then to its first nearest neighbor in the Li sublattice Path C: represents the diffusion within the basal plane Both path A and B converge to the energy favorable path indicated in Fig ure 3 1 1 during the NEB calculation, in which the relaxation of the ions is considered The diffusivity of lithium vacancy can be calculated using the standard diffusion formula 1 T k H DB vacexp0 2 [3 2 3 ] PAGE 85 85 w is the 0 is the vib ration frequency of the atoms, H is the migration energy, kB is the Boltzmann constant, and T is the temperature. 0 is assumed to be the same as the thermal energy kBT, where h is the Plank constant. At 1200 K, in which most diffusivity are 0 2.5 x 1013 s1. T he jump distance between the first nearest Li neighbors is 3.771 This yields a diffusivity at 1200 K of 4.2 x 109 cm2/s for an activation enthalpy of 1.63 eV, as obtained from DFT. The diffusivity calculated from experimental data at 1200 K is 1.6 x 107 cm2/s ( Mehta et al. 139), 2.8 x 108 cm2/s ( Halstead et al. 140), and 2.6 x 107 cm2/s ( Schmidt et al. 141) respectively (Table 3 11). Table 3 11. Activation Energy (Ea) obt ained from first principles calculations, compared with the experimental results. The pre factor (D0) of the calculation is based on the 0 2.5 x 1013 s1 and a jump distance between the first nearest neighbors of 3.771 The diffusivity is calculated at T=1200K. Calculation Experiment (This Work) Mehta et al 139 Halstead et al 140 Schmidt et al 141 D 0 (cm 2 /s) 3.56 10 2 5.10 10 1 1.80 10 1 2.10 10 2 E a (eV) 1.6 3 1.55 1.62 1.17 Diffusivity (cm 2 /s) 5.1 0 10 9 1.58 10 7 2.84 10 8 2.57 10 7 These variations in diffusivity mainly reflect the differences in the activation energies. In particular, the activation energy obtained by Schmidt et al. is significantly smaller than the others are, which lead to a higher diffusivity. The difference is even larger at low temperatures. At room temperature (298 K), the diffusivity is 2.4 x 1030 cm2/s, based on the DFT results. The diffusivity calculated from experimental data at room temperature is 3.2 x 1027 cm2/ s ( Mehta et al 139), 7.4 x 1029 cm2/s ( Halstead et al. 140), and 3.6 x 1022 cm2/s ( Schmidt et al. 141) respectively. These are extremely low diffusion constants, which indicate that the defect clusters are rather stable at such temperature. PAGE 86 86 The average diffusion distance can then be calculated using the standard random walk equation: Dt x 6 [3 24] w here x is the diffusion distance, D is the diffusivity and t is the time. At room temperature, for lithium vacancy to diffuse in the order of angstroms, it requires the diffusion time to be in the range of years. This indicates that each configuration of the defect complex should be very stable at room temperature, which may lead to significant local dipole moments, even including non uniaxial contributions. However, the average diffusion distance is predicted to be in the order of tens of nanometers after five hours of annealing at 250 C 143, which indicates the lithium vacancy should still be able to diffuse and change the configurations of the defect complexes. This is in agreement with the experimental observation that the relative peak emission intensities measured by combined excitation emission spectroscopy can be modified by the annealing of the sample 143. 3 .8 Discussion As discussed in the section 3.1 the X ray and NMR measurements do not provide a consistent picture of the defects in LiNbO3. The X ray data indicate that the oxygen vacancy concentration is negligible, i.e., that there are no Model I defects (VO ..+ 2 VLi ); by contrast, NMR suggests that th ere are almost equal numbers of Model I and Model III defects ( NbLi + 4 VLi ). It is of considerable interest to analyze our results to see if they can shed any light on this issue. With reference to the thermodynamics analysis, it is clear that room temperature and standard pressure, ( O [LiNbO3] = O [gas ]), are most relevant for applications; thus the following analysis is performed under these conditions. However, the DFE still cannot be PAGE 87 87 uniquely determined since the separate values of chemical pot ential of Li and Nb cannot be determined uniquely, but are coupled. Our calculations are consistent with the experiments in that they indicate that both Model I and Model III are energetically favorable relative to most other defect structures. In particular, as shown in Figure 3 13(a), Model III and the Li Frenkel are the dominant defects under different conditions, while Model I has the second lowest DFE through most of the chemical potential stability range. The error bars of 0.2 eV /defect represent our estimated uncertainty in the formation energies. The Li Frenkel defects do not, of course, change the stoichiometry of the system and are thus not relevant with regards to interpreting the X Ray and NMR data. The concentrations of each defect can be calc ulated from: ) / exp( kT G N Cf site d [3 2 5 ] where Cd is the number of the defects, Nsites is the number of sites in the crystal where the defect could occur, k is the Boltzmann constant, and T is the temperature in Kelvin. The calculated defect concentrations are shown in Figure 3 13(b) This figure must be interpreted with caution, since negative formation energy would apparently lead to a concentration of unity for the corresponding defect, w hich is not reasonable, since the calculations on which Figure 3 13 are based assume a dilute solution of defects; presumably the formation energies will increase with increasing concentration, such that the system remains stable. For low Nb chemical potentials (Lirich conditions), the Li Frenkel is the dominant defect, with a small admixture of Model I. At higher Nb chemical potentials, Model III is dominant and produces a concentration of lithium vacancies that is at least two orders of magnitude larger than produced by Model I. Since Model I defects lead to oxygen vacancies while Model III defects do not, this would be consistent with the experimental finding 24 that for Nb rich conditions, oxygen vacancies PAGE 88 88 concentration is much smaller than the lithium vacancy concentrations. However, it should be noted that although the defect concentration of Model I is much smaller than Model III, it may still reach substantial levels under these Nbrich conditions. In order to achieve the experimentally postulated ratio of 1.1: 1 for Model I and Model III, the DFEs of these tw o defect reactions would have to be much more similar in magnitude than they are calculated to be; thus this study does not support this scenario. Figure 3 13. The DFEs of M odel III Li Frenkel and M odel I as a function of niobium chemical potential. (b) The defect concentration caused by model III Li Frenkel and model I as a function of niobium chemical potential. The chemical potential unit is in eV /atom We finally address the mechanism for the vapor trans port equilibration (VTE) method 18 used widely in transforming a congruent composition of lithium niobate into a stoichiometric composition. In the vapor transport equilibration method 20, a target wafer is exposed to Li2O vapor produced by decomposition of Li2CO3 at high temperature (> 730 oC) 144. This produces a constant vapor pressure of Li2O equal to the vapor pressure of stoichiometric LiNbO3. The wafer, which is congruent LiNbO3 and deficient in Li2O, absorbs Li2O from the vapor u ntil it PAGE 89 89 reaches the stoichiometric composition. Then Li diffuses into a considerable depth due to the high diffusivity. The whole growth procedure is thus under Li2O rich conditions In the language of our calculations, this experimental VTE procedure thus corresponds to using Li2O as the reference state. Therefore, according to the analysis in section 3.5 the lithium Frenkel defect is predicted to be the most stable with a positive DFE. Using Equation 3 25, the concentration of defect is predicted to be negligible, which result in a stoichiometric composition. Our calculations thus indicate that the relative energetics of the defect clusters can change under variations of temperatures and chemical potent ials during annealing experiments which may result in different relative abundances 143. Furthermore, these results strongly suggest that the defect chemistry can be expected to b e sensitive to the thermal and chemical history of the sample mak ing direct comparison s among experimental results difficult; indeed, this sensitivity could be the origin of some of the experimental inconsistencies. In summary, the intrinsic defects and defect clusters in LiNbO3 have been characterized using DFT PAW approach. The formation energies of various types of defects and defect clusters have been reported and the dominant defects are determined over the chemical potential stability range. The ass ociation effects have also been investigated. The Li Frenkel was found to have the lowest DFE under Li2O rich conditions, while the cluster consisting of a niobium antisite compensated by lithium vacancy ( NbLi + 4 VLi ) was found to be energetically stable under Nb2O5rich conditions. These results are consistent with experimental observations. In particular, for the first time, the intrinsic defect cluster mentioned above, which was conjectured before to play a critical role in domain structure and dynamics, has been verified using first principles calculations to indeed have the lowest formation energy. The defect analysis also sheds light on the mechanism by PAGE 90 90 which a congruent composition of LiNbO3 can be transformed into a stoichiometric compositi on by the widely used VTE method today. PAGE 91 91 CHAPTER 4 EXTRINSIC DEFECTS IN LITHIUM NIOBATE 4 .1 Background of Extrinsic Defects During the past decades, there have been numerous studies of doped LiNbO3 43, 45, 143, 145147 for optical lasers, optical amplifiers, and integrated optical circuits 148 151. Va rious dopants have been added to the system for different purposes. For example, Mg has been introduced into the syst em to increase the resistance to photorefractive damage 25, 26; Er and Nd have been added for solid state laser applications 28 31 Dopant ions have also been employed as probe s to investigate the structure of domain walls and defect/domainwall interactions 3234. Various experimental techniques, including electron spin resonance 4, X ray standing wave analysis (XSW) 39, 40 Rutherford backscattering 38, extended x ray absorption fine structure (EXAFS) 41, and ion beam channeling 152 have been employed to investigate the s ite selectivity of dopant ions Moreover, several optical and magnetic resonance spectroscopy studies 43 47 have been used to determine the local environment s and configurations around the dopant sites. While some of these studies used congruent LiNbO3 samples grown by the Czochralski method 13, 14, others used stoichiometric LiNbO3 produced through vapor transport equilibration (VTE) 18 20 from congruent samples. The influence of the sam ple stoichiometry on the site selectivity and distribution is also unknown, which makes the comparison between experiments on samples from different synthesis technique problematic. In addition to understanding how the dopant s are incorporated into the s tructure, an understanding of the associated charge compensation mechanisms is also necessary. Atomic level simulations 153 using an empirical potential have been used to predict the dominant defect for various dopant ions. However, the empirical potential s did not explain the discrepancies between experiments For instance erbium sitting on lithium site compensated by erbium sitting PAGE 92 92 on niobium site is predicted to be dominant at 0 K by the empirica l potential. However, a t room temperature, the same potential predicted that Er3+ ions occupy the Nb site, compensated by NbLi .. It is not clear how such a small temperature change (300 K = 0.026 eV) could cause such a qualitative change in the dominant crystal defect. The experimental and theoretical results thus indicate that the charge compensation mechanism s for dopant ions in LiNbO3 are still not well understood In the present study, the formation energies of several dopants with different valance charges (Mg2+, Fe2+, Er3+, Nd3+, and Fe3+) have been determined using DFT to predict the site preference and corresponding charge compensation mechanisms The influences of the composition and growth condition on the defect chemistry have been assessed based on thermodynamic analysis under Nb2O5 and Li2O reference states respectively Furthermore, the cha rge transfer level of iron with regard to the band structure has been calculated to elucidate the effects of coexist ence of both 2+ and 3+ charged ions on the optical properties of LiNbO3. 4 .2 Structure of Reference Binary Oxides In order to calculate the formation energies and to determine the dominant defects, the chemical potentials of the dopants are needed However, it is extremely difficult to determine the exact value of the chemical potential, as it is a variable depending on many factors, including partial pressure, temperature, and compositions Therefore, appropriate reference states are needed, such as the metallic or oxides phase of the dopants For the metallic phase to be used it is assumed that materials are in an extremely reducing environment which is far from the condition of dopants in LiNbO3. A more general and plausible way is to use the co rresponding binary oxides as reference state, in which the bonding environments of dopant s is more similar to PAGE 93 93 that in LiNbO3. Therefore, the structures of the binary oxides need to be clarified first. In current study, the thermodynamic stable and most com mon phase for each binary oxide is used. 4 .2.1 Crystal Structure of MgO and FeO Both MgO and FeO has a rock salt structure (Space Group: Fm 3m) 1, in which both cation s and anion s form a separate FCC sub lattice with lattices interpenetrating with each other. The ionic radii of Mg2+, Fe2+, and O2 are 0.86 0.75 and 1.26 respectively 154. T he ratio of cation/anion radii of MgO and FeO are 0.683 and 0.595 respectively, which fall i nto the range between 0.414 and 0.732. Therefore, each cation or anion has a coordination number of six and is at the vertices of an octahedron. An illustration of atomic arrangement for rock salt structure is shown in Figure 4 1 Figure 4 1 Crystal structure of MgO and FeO. PAGE 94 94 4 .2.2 Crystal Structure of Er2O3 Nb2O3 and Fe2O3 Erbium oxide has a cubic structure, with space group of Ia3 155. The lattice co nstant of Er2O3 is 10.5 155. There are 32 erbium and 48 oxygen in each unit cell (Figure 4 2 ), resulting in a total number of 80 ions. Erbium ions sit on two different Wyckoff positions, 8a and 24d 155 respectively. The details about the Wyckoff positions and fraction coordinates are given in Table 4 1 Er3+ is coordinated with six ox ygen ions and in a distorted octahedron. O xygen has a coordination number of four and occupie s a tetrahedral site. Figure 4 2 Crystal structure of Er2O3. Table 4 1 Wyckoff position and fraction coordinates of Er2O3 155. Er 2 O 3 Wyckoff x y z Er (1) 8a 0 .0000 0 .0000 0 .0000 Er (2) 24d 0.2855 0 .0000 0.25 00 O 48e 0.1293 0.1471 0.9165 The structure of neodymium oxide belongs to the trigonal system (Space Group: P 3m1 156, Figure 4 3). The lattice constants in the [1000] and [0001] direction are 3.8277 and 5.9908 PAGE 95 95 156. The Wyckoff position of Nd3+ is 2d while the oxygen occupy two different Wyckoff site, 2d and 1a respectively. The details of the Wyckoff position and fraction coordinates are given in Table 4 1 Figure 4 3 Crystal structure of Nd2O3. Table 4 2 Wyckoff position and fraction coordinates of Nd2O3 156. Nd2O3 Wyckoff x y z Nd 2d 0.3333 0.6667 0.2462 O 2d 0.3333 0.6667 0.6646 O 1a 0 .0000 0 .0000 0 .0000 Fe2O3 has several phases 157, with corundum Fe2O3 being the most common form 158. It has a trigonal Barvais lattice with space group: R 3c The oxygen ions form hexagonal close packed array with iron sitting on two thirds of the octahedral sites. There are six formu las of units in each unit cell (Figure 4 Fe2O3 is a metastable cubic face centered phase, which converts to the alpha phase at ~500 oC. Fe2O3 is also metastable and cubic and converts to the PAGE 96 96 alpha phase at even higher temperature. To be used as Fe2O3 is chosen because it is the most common and thermodynamically stable phase. Figure 4 4 Crystal Fe2O3. 4 .3 Computational Details The system size for doped LiNbO3 study is 2 x 2 x 1 unit cells (each of 30 atoms), which for perfect LiNbO3 contains 120 atoms and 720 electrons; periodic boundary conditions are applied in all three dimensions. The i ntegration over the Brillouin zone uses a 4 x 4 x 2 Monkhorst Pack 117 k point mesh. Mg has an electronic configuration of [Ne ]. 3s2, in which the two s electrons are considered as valence electrons. This electronic configuration is much simpler than that of the Fe, which is a transition metal with d orbitals For iron, the electronic structure is [ Ar ]. 3d6. 4s2. Both 3d and 4s electrons are treated as valence electrons. Er and Nd are rare earth elements. The electronic structure of Er is [ Xe ] 4f126s2. Due to the localized nature of 4f electrons 159, 160, the frozen core approximation 161, 162 is used. For Er, PAGE 97 97 eleven of twelve f electrons are treated as the frozen core in VASP with only a single f electron treated as a valence electron 163. In addition, the 6s2 and 5p6 are also tre ated as valence electrons, resulting in nine Er valence electrons in total. The electronic configuration of Nd is [ Xe ] 4f36s2, similar to Er; frozen core approximation has also been implemented in VASP. The localized nature of 4f and 3 d electrons results in strong correlation effects. In general, mean field theories, such as the local density approximation (LDA) and GGA, cannot accurately describe the strong correlations of d and f electrons 163165. Therefore, in addressing this problem, either a higher level theory, or a correction to LDA or GGA, such as DFT+U approach 78, 79 is required. The effect of the +U term on defect energetics need s to be assessed. To evaluate the effects of +U, se veral fundamental properties, which include lattice constants, band gap, and magnetic moments, have been calculated using GGA and GGA+U for the reference binary oxides To perform GGA+U calculations, reasonable values for the U and J parameters are needed. In this study Dudarevs approach 79 has been used, in which only the difference between U and J is meaningful. There are no parameters of + U f or Nd2O3 in the literature Only one parameter of +U, in which U J = 6.5 eV, has been applied to Nd, in that case TiO2 using LDA+U. 166. For Fe2O3, there are several studies focusing on parameterization of U and J value s U J = 4 eV has been chosen based on the analysis of heat of formation data 167, which is also probably the most widely used value in the literature. However, this value is used with LSDA+U. As for GGA +U, both U J = 4 eV and U J = 6 eV have been used in the previous st udies 168, 169. Since this section is trying to evaluate the effects of +U, both values have been employed to calculate the fundamental properties of FeO Fe2O3. PAGE 98 98 T he fundamental properties of FeO and Fe2O3 are calculated using GGA and GGA+U and compared with experimental data For FeO, it is observed tha t GGA underestimates the lattice constant by ~ 0.21 while GGA+U using both values overestimate s it by 0.08 GGA predict s its band structure to be metallic, which is a well know n failure of this method 62. On the other hand, GGA+U predict s FeO to be semiconductor with a band gap of 1.7 eV and 2.02 eV respecti vely using U J = 4 eV and U J = 6 eV Although the calculated band gap using GGA+U is lower than the experimental value (2.4 eV 170) it correctly predict s the material to be semiconductor and with an typical of GGA In addition, GGA also severely underestimates the magne tic moments while GGA+U improves it significantly. Table 4 3 Fundamental properties comparison between GGA, GGA+U and experiments for FeO. FeO GGA GGA+(U=4) GGA+(U=6) Experiment Crystal Structure Fm 3m Fm 3m Fm 3m Fm 3m 1 Lattice Constants () 4.088 4.382 4.382 4.3 34 171 Magnetic Structure AFM AFM AFM AFM 172 Magnetic moments 1.009 3.755 3.8 4. 2 173 Band Gap (eV) Metallic 1.7 2.02 2.4 170 Fe2O3, GGA+U has a similar effect. GGA predicts a lower value for the lattice parameter than the experimental results. The lattice constant predicted using GGA+U is in excellent agreement with experimental value. The magnetic moment calculated using GGA is too low while GGA+U improves it to 80% of the experimental value. Table 4 4 Fundamental properties comparison between GGA, GGA+U and experiments for Fe2O3. Fe 2 O 3 GGA GGA+U(4 eV ) GGA+ U( 6 eV ) Experiment Crystal Structure R 3c R 3c R 3c R 3c 174 Lattice Constants () a 4.698 5.02 5.018 5.0355 175 Lattice Constants () c 13.376 13.716 13.685 13.7471 175 Magnetic Structure AFM AFM AFM AFM 174 Magnetic moments 0.734 4.071 4.222 4.9 176 It is observed that GGA+U using both values generate similar results and the difference between two sets of parameters is relatively small. Since U J = 6 eV generally provides slightly PAGE 99 99 better lattice constants, band gap, and magnetic moments. Therefore, this value is used for the further defect formation energies calculations. The comparison of fundamental properties between GGA, GGA+U and experiments for Nd2O3 is given in Table 4 5. Unlike iron oxides, no significant difference in lattice constants and band gap has been observed between GGA and GGA+U. The lattice constants predicted by GGA and GGA+U are very close to the experimental value, while both methods u nderestimate the band gap by ~ 1.6 eV. Therefore, there is no need to use GGA+U for Nd2O3. Table 4 5 Fundamental properties comparison between GGA, GGA+U and experiments for Nd2O3. Nd 2 O 3 GGA GGA+U ( 6.5 eV ) Experiment Crystal Structure P 3m1 P 3m1 P 3m1 156 Lattice Constants () a 3.839 3.934 3.8827 156 Lattice Constants () c 5.998 6.107 6.077 156 Band Gap (eV) 4.12 4.24 5.8 N o +U parameters have been determined for Er in either LiNbO3 or Er2O3. However, there are several sets of parameters available for Er in GaN and in ErxGa1xN 177, 178. For Er in GaN using the density functional based tight banding (DFT TB) method 177, U J = 10.3 eV was used. However, it is known that the DFT TB method yields higher es timates for U and J than does regular DFT 177. For ErxGa1 xN using LSDA+U 178, U=8.6 eV and J=0.75 eV have been used. Since there are no accepted values for U and J for LiNbO3, we have explored a range in order to map out the possible magnitude of the effects of electron localization: (i) U J = 4 eV; (ii) U J = 7.85eV 178, and (iii) U J = 10.3 eV 177. Although none of these parameters may actually be the best physical representation of electron localization for Er in LiNbO3, they do span the range of reasonable values, and should thus enable us to assess the importance of the localization of electrons on the defect energetics. Because the GGA+U has to be used for both the reference PAGE 100 100 state Er2O3 and for Er in LiNbO3, its effects on Er2O3 are discussed first, followed by its effects of Er in LiNbO3. For Er2O3, the band gap calculated using pure GGA is 4.37 eV, which is about 13% less than the experimental value of ~ 5 eV 179. This agreement can be considered to be good, because it is well known that for many systems GGA fails to predict the band gap of a material with partially filled d or f orbital 163 165. There is no significant change in the band gap or electronic density of states fo r any of the parameter sets for U J, which indicates that the Er is reasonably described by GGA (Figure 4 5). Figure 4 5. Electronic density of States of Er2O3 calculated using GGA and GGA+U. The parameters U J = 4 eV, U J = 7.85 eV, and U J=10 .3 eV are considered. For Er in LiNbO3, the representative case of ErLi ..+2 VLi is chosen to illustrate the effects of +U, in which the association effect is also considered. For U J = 4 eV U J = 7.85 eV, and U J = 10.3 eV, the DFE increases by 0.067eV/defect, 0.102 eV /defect and 0.072 eV/defect, PAGE 101 101 respectively. Thus, the stability order predicted from the GGA calculations is correct and the energy values can be considered reliable to within ~0.1 eV/defect. We thus conclude that electron local ization as captured by the +U term has no significant effect on either the electronic structure or the defect energetics for Er systems. In summary, for MgO, no +U is needed as Mg does not have any f or d electrons. For iron oxides, it is necessary to use the GGA+U method to predict reasonable band gap and magnetic moments. Based on the analysis of the lattice constants, band gap and magnetic moments of FeO Fe2O3, U J = 6 eV is chosen to carry out the formation energies calculations. For Er2O3 and Nd2O3, no significant effects of the +U on the lattice constants and band gap of the reference binary oxides or on the formation energy of the impurities in LiNbO3 has been observed. Therefore, GGA will be continued to perform the formation energy calculations. 4 .4 Formation Energies and Site Preference While defect structures are generally written in terms of fully charged constituent defects, partially charged and charge neutral defects are also possib le 118. For example, many experiments on TiO2 have concluded that as the oxygen partial pressure decreases from moderately reducing conditions to extremely reducing conditions, there is a tr ansition from a fully charged oxygen vacancy to a singly charged oxygen vacancy or triply charged titanium interstitial 133. Thus, b efore characterizing the defect clusters, it is important to analyze the individual defects with various charge states. The following sub section s will discuss the results obtained for individual defects of various charges, followed by the defect clusters For the defect clusters, t he constraint of charge neutrality leads to the simplification of not having to consider the effects of the location of the Fermi energy 136. PAGE 102 102 4 .4.1 Single Defects S ingle isolated point defects in the 2x2x2 supercell are investigated which corresponds to the s ituation of the charge compensating defects being far from both the impurities and each other. Assuming the impurities are + 3 charged, depending on the position of ions two different types of point defects could be formed, which are MLi and MNb respectively. To explore how the Fermi level may affect the point defects in the system, an exemplar case of Er in LiNbO3 is discussed in detail in the following. T he DFEs of ErLi and ErNb defects with various charge states in LiNbO3 have been calculated using Equation 3 22. T he influence of the Fermi energy on the stability of each individual defect is considered The reference zero of the Fermi energy is assigned to be the valence band maximum (VBM) of the perfect structure; the highest Fermi energy thu s corresponds to the conduction band minimum (CBM), which for GGA is 3.5 eV, close to the experimental value of 3.78 eV 115, 116. A Fermi energy at the center of the gap corresponds to the pure system, while p type and n type correspond to Fermi energies closer to the valence and the conduction bands, respectively 180. As the Fermi level could shift in the presence of impurities and defects, two extreme cases corresponds to f=0 and f=3.5. The dependence of DFEs on the Fermi lev el is shown in Figure 4 6 using Nb2O5 (the BC line in Figure 3 4) as the reference state. At each value of the Fermi level Figure 4 6 includes only the charge state with the lowest DFE of each individua l defect. The slopes of the lines in Figure 4 6 represent the charge states of the defects .As the Fermi energy increases the thermodynamically stable charge of the cations decreases. ErLi .. has a lower energy than its singly charged counterpart ErLi .. fo r most of the range of Fermi level, while ErNb has the lowest DFE throughout the entire range. The partial density of states of Li, Nb, O, and Er for ErLi .. and ErNb PAGE 103 103 have been also generated and are given in Figure 4 7 No significant effect of Er on the PDOS of lithium, niobium and oxygen is observed. From Figure 4 7 it can be seen that there is a greater overlap in the PDOS of the Er and O ions for ErNb than that for ErLi .., indicating stronger covalent bonding. In the unit cell of LiNbO3, six closed packed planes are stacked 10. The cations lie in the oxygen octahedral in repeati ng sequence: Li, Nb, vacant site, Li Nb, vacant site 10. Within a single plane of cations, there is an ordered arrangement of Li ions, Nb ions, and vacant sites. For ErLi .. the DFT calculations show that the Er is shifted along <0001> towards the vacant site by ~0.17 The observation by Gog et al. using xray standing wave (XSW) spectroscopy 39, 40 showed the Er to be in a position near a Li site but displaced normal to the basal plane by 0.46 39. This difference between the experimental measurement and calculated val ue for the displacement may have several origins. First, the experiments were performed using congruent LiNbO3, while the simulations were carried out for stoichiometric material. Second, the experiments measured the near surface region where Er is incorpo rated through diffusion, while the calculations are for the bulk. Both the composition and surface effects could cause a change in the position of the Er ions. Third, the difference may be also due to the limitation of the current calculation method, which simulates the defects in a finite size supercell and assumes that the defect concentration is in the dilute limit. For partial charged defects, ErLi and ErLi x, the displacements are 0.14 and 0.08 respectively. By comparison, the DFT calculations of ErNb yield a displacement of only ~0.02 towards the vacant site. The displacements for ErNb and ErNb x are less than 0.005 PAGE 104 104 Figure 4 6. Defect formation energies of the lowest energy charge state of ErLi and ErNb as a function of Fermi energy; the lowest energy charge state is given in each case. The Fermi energy ranges from EF=0 (left) at the VBM to EF=3.5eV at the conduction band minimum (right). (a), using Nb2O5 as reference state. (b) using Li2O as reference s tate. PAGE 105 105 Figure 4 7 (a) Partial density of state (PDOS) for ErLi .. ; (b) PDOS for ErNb The absolute values of the number of the DOS in the two figures have been rescaled to make the comparison easier. PAGE 106 106 These results are readily understood in terms of the electrostatic and elastic interactions within the system. Thus, when Er substitutes on a Li site, the repulsive electrostatic interaction between Nb and Er is greater than that between Nb and Li and thus forces the Er towards the vacancy. The ionic radii of Li+, Nb5+ and Er3+ are 0.76 0.64 and 0.89 r espectively 154. Therefore, the larger size of the Er ion also tends to push it closer to the vacancy; thus, both interactions push the Er towards the vacant site. By contrast, when Er substitutes on a Nb site, the smaller Er charge leads to weaker repulsion from the Li, thus tending to move the Er ion away from the vacant site. However, the larger ionic radii will again push the Er towards the vacant site. The electrostatic and elastic interactions thus tend to counteract each other, resulting in only a small net displacement. 4 .4.2 Defect Clusters There are t wo different charge states of the impurities considered in current paper, which are M+2 (Mg2+ and Fe2+) and M+3 (Er3+, Nd3+, and Fe3+). Depending on the charge states, different defect reactions are needed in order to maintain the overall charge balance of the system resulting in different ratio of the charge compensation species in the system Based on the previous studies and literatures, the following defect clusters have been considered in this study For M3+: (i) MLi .. + 2 VLi (ii) MLi .. + MNb (iii) 2 MNb + NbLi For M2+: (i) MLi + VLi (ii) 3 MLi + MNb (iii) 4 MNb + 3 NbLi As in the cases of the single point defect analysis, an exemplar case of Er is pres ented below in full detail. Similar analyses for Mg, Nd, and Fe are given in the later part of this section, with comparison between these impurities. PAGE 107 107 For Er related defects, t he current paper focuses on the following four extrinsic defect clusters that have been discussed in the literat ure: 2 ErLi ..+ 4 VLi = 2( ErLi ..+ 2 VLi ), ErLi ..+ ErNb and 2 ErLi ..+ NbLi 10, 153. T wo completely separated ErLi ..+ 2 VLi cluster are considered because this approach allows all of the defect clusters to be analyzed independently of the chemical potential of the Er. The energies of neutral defect complexes are first determined from the formation energ ies of the individual defects As mentioned above, because these defects are charge neutral, their energies are independent of the location of the Fermi level. However, the DFEs of the clusters do depend on the choice of reference state. From Figure 4 8 it is can be seen that the formation energies of the ErLi ..+ ErNb and 2 ErLi ..+ NbLi clusters are the independent of the reference state. This is because within LiNbO3, the sum of chemical potentials of Li and Nb is constant for any given value of the chemical potential for oxygen. By contrast, the DFE of the 2 ErLi ..+ 4 VLi and 2 ErNb + 4 VLi +2 NbLi clusters do depend on the chemical potential of oxygen. Thus, the ir formation energies are directly influenced by the reference state of chemical potential. Table 4 6 Defect formation energies for defect reactions under different reference states. The chemical potential of Er is assumed to be the same as the value in Er2O3. These value are calculated based on the single defect energy without considering association effects. Defect Reaction DFE(eV) Fully Charged Li 2 O Nb 2 O 5 2 Er Li .. + 4 V Li 8.369 4.607 Er Li .. + Er Nb 1. 9 39 1. 9 39 2 Er Li .. + Nb Li 4.274 4.274 2Er Nb + 4 V Li +2 Nb Li 13.039 0.064 The calculations show that the formation energy for 2 ErLi ..+4 VLi is negative for the Nb2O5 reference state This result is similar to the previously studied intrinsic NbLi .+ 4 VLi cluster Under Nb2O5rich conditions the 2 ErNb + 4VLi + 2 NbLi cluster has a very small positive PAGE 108 108 formation energy, and might thus be expected to also be observed. The calculations also predict that for the Li2O reference state, ErLi ..+ ErNb has the lowest DFE. However, the DFE is positive, which indicates energy is required to create the defects under conditions such as VTE. Based on the thermodynamic calculation, the concentration of such defects at room temperature should be very low. F igure 4 8. Formation energies of Er related defect clusters under Li2O rich conditions (blue) and Nb2O5 rich conditions (red ). To this point, the calculations have assumed that the individual defects that form a defect cluster are infinitely far away from each other; that is, there are no interactions between them. However, due to the charge and elastic interactions induced by the differences in ionic radii, defect association can be expected to change the DFEs. Here, the case of ErLi ..+2 VLi is analyzed While there are an almost innumerable number of possible different compensating structures, it is only necessary to analyze the limiting case of association, in which some of the compensating PAGE 109 109 vacancies are in the first neighbor shell while all of the oth ers are far from the Er defect. Cases in which vacancies lie in the second or third shell or beyond can be assumed to have lower association energies. T here are six Li sites in the first nearest neighbor (FNN) shell around the Er occupied Li site, as indic ated in Figure 4 9 Figure 4 9. Possible lithium vacancy positions in the first nearest neighbor (FNN) around Er sitting on lithium site. The oxygen sub lattice is not shown. The arrangement of the Li vacancies can conveniently be divided into three ca tegories: A0: No lithium vacancies in the FNN shell; A1: One lithium vacancy in the FNN shell; A2: Two lithium vacancies in the FNN shell; In the A0 structure there is no association of the vacancies; this is the case analyzed above and found to have negative formation energy for the Nb2O5 reference state. For A1, there are two distinct configurations for one lithium vacancy sitting around the Er site. The lithium PAGE 110 110 vacancy could lie above the Er ion (4 or 5 or 6 in Figure 4 9 ) or below the Er ion (1 or 2 or 3 in Figure 4 9 ). For the A2 case, there are C6 2 = 15 possible configurations of vacancies. However, due to the three fold rotational symmetry around the ErLi .. only four of these configurations are crystallographically distinct. The four different a rrangements are 1 2, 1 4, 4 5, and 2 5 in Figure 4 9 As illustrated in Figure 4 10, the energy of the A1 structure is lower than that of the A0 structure, while the energy of the A2 structures are lower yet by 0.068 eV/defect (i.e., 0.204 eV for the total cluster). These results indicate an attractive interaction between ErLi .. and VLi which is consistent with the conjecture of Dierolf et al 143 that ErLi .. is c ompensated by a VLi close to ErLi ... Figure 4 10. Defect formation energy for ErLi .. + 2 VLi shows the effects of association of the lithium vacancies around ErLi .. site. PAGE 111 111 Since this association energy is relatively small (~2.5 kBT), it might be supposed that the structures of a defect cluster could fluctuate dynamically among various symmetry equivalent configurations However, our DFT calculations show that the migration barrier for the diffusion of a lithium vacancy between the first nearest ne ighbor in the lithium sublattice is 1.63 eV which is consistent with the experimental value of 1.55 eV 139 and 1.62eV 140 for the ionic conductivity. These values of activation energies indicate that the transitions among configurations are relatively difficult. However, lithium vacancy can still move a considerable distance during long time annealing. For example, experimentally st oichiometric LiNbO3 has been annealed at 250 oC for 5 hours 143. I0 is equal to the thermal energy kBT. T he jump distance between first nearest neighbors is 3.771 Under such condition, the diffusivity of lithium vacancy is calculated to be 7.07 x 10 18, which compares reasonably well with previous estimates of 5.97 x 10 16, (Mehta et al 139) and 4.46 x 10 17 (Halstead et al 140). The average diffusion distance can then be calculated in the usual manner : Dt x 6 [4 1 ] where x is the average diffusion distance, D is the diffusivity and t is the diffusion time. Using the above value of diffusivity, the average diffusion distance is 5.7 nm (this work), 80.3 nm ( Mehta et al. 139), and 21.9 nm ( Halstead et al. 140). The average diffusion distance in the order of tens of nanometers indicates the lithium vacancy can still travel and change the configurations of the defect complex. The calculation is in agreement with the experimental observation that the relative peak emission intensities have been modified by the annealing of the sample 143. PAGE 112 112 Similar analyses have been applied to Mg, Nd, and Fe impurities. A summary of the formation energies of these impurities under Nb2O5 reference state with possible different charge states is given in Figure 4 11. Figure 4 1 1 Defect formation energies of various possible defect reactions for impurities under Nb2O5 reference state. It is noted that under Nb2O5 reference state, the impurities sitting on the lithium site compensated by the lithium vacancies have the lowest formation energies compared with other types of defect clusters. Similar to the intrinsic defects situation, the formation energies are negati ve, which indicates the Nb2O5 may not accurately represent the experimental conditions. However, as the reference state moves from Nb2O5 towards Li2O, the formation energies of the defect clusters increases and impurities on lithium site compensated by lit hium vacancies remains dominant over a certain range of the chemical potentials. This is because the relative stability orders of defects will remain as chemical potential change all of them in a similar way. PAGE 113 113 It is found that the formation energy for Mg s itting on both sites is also negative albeit less negative than the defect cluster consisting of Mg on lithium site compensated by lithium vacancies, unlike the other impurities. This indicates the tendency for Mg to sit on both sites is higher than for Fe Er and Nd impurities. As for iron, both Fe2+ and Fe3+ on lithium site compensated by lithium vacancies are possible, as formation energies of both defect reactions are negative. However, Fe3+ is more energetically favorable. This suggests that Fe doped LiNbO3 may have different charge states of iron and the concentration of Fe3+ is higher than Fe2+. The formation energies of impurities gi ven in Figure 4 12 are under Li2O reference state, which corresponds to the experimental condition of stoichiometric LiNbO3. Figure 4 1 2 Defect formation energies of various possible defect reactions for impurities under Li2O reference state. As the figure shows, all the defects have positive formation energies, which is consistent with the stoichiometric composition as the defect concentration is very small based on energy PAGE 114 114 calculations. By comparing different defect clusters, it is found that impurities occupying both sites have the lowest energy for both M2+ and M3+. Similar to the Nb2O5 cases, Fe3+ is more energetically favorable than Fe2+. The detailed discussion between these calculation and experimental evidence is given in section 4.6. 4 .5 Charge Transfer Levels A charged crystallographic defect associated with one or more electrons could absorb visible light, so that usually transparent materials could become colored. This involves a charge transfer process, which generally absorb s or emit s certain wave length of light The charge transfer level, also called the optical level, is similar to the thermodynamic transition level defined previously. However, they differ in the atomic configuration of final charge state. For the thermodynamic transition level, the atomic configuration of final charge state is fully relaxed. For the charge transfer level, the final state is calculated b a sed on the configuration of the initial charge state This is because the electron movement is m uch faster than the atomic relaxation and the transition from one charge state to the other can be considered spontaneous. The concept of calculating charge transfer levels is shown in Figure 4 13. Electrons in the conduction band can recombine with a hol e on the acceptor, which leads to emission of a photon with energy EPL. This photon can be experimentally observed and measured, for instance, using photoluminescence or absorption spectroscopy. During this process, the atomic configuration of the donor re mains unchanged. The energy difference between this configuration and fully relaxed configuration is the relaxation energy, Erel, which is also termed as Franck Condon shift 118. EA is the t hermal ionization energy of the acceptor 118 and Eg is the band gap of the bulk material. PAGE 115 115 Figure 4 1 3 Schematics of calculations of charge transfer levels. EPL is the charge transfer energy, Erel is Franck Condon shift, Eg is the band gap, and EA i s the thermal ionization energy 118. The thermal ionization energy, or the energy required to remove one electron from dopant at a particular charge state in LiNbO3, can be calculated using the following equation: V E E M E M E Ev corr tot tot ALi Li ] [ ] [ [4 2 ] w here Etot[ ErLi ..] is the energy of the supercell containing an dopant at lithium site with +2 charged Etot[ ErLi ] is energy of the supercell containing an dopant at lithium site with +3 charged Ecorr is the energy difference between the highest occupied state at the point and the special k point where the band energy is lower than at the point 118. This term is needed because the electrons should be removed from the top of the valence band at the point while the electron is actually taken out of the highest occupied KohnSham level at the special k point in the supercell calculations. Ev is the valence band maximum of the bulk materials in the PAGE 116 116 calculation, and V is the electrostatic potential difference between the two supercell s The thermal ionization energy of Fe2+ in LiNbO3 is calculated to be 1.836 eV, which is consistent with the previous estimates 181. The comparison between theoretical calculations and experimental PL measurement of Fe doped LiNbO3 is given in Table 4 7 It is seen that the value of EPL measured from PL experimental is ~ 1.4 eV 182. The charge transfer level (EPL) of Fe2+/Fe3+ has been calculated to be is 1.34, which is in excellent agreement with experime ntal value. This confirms the longtime conjecture that there is a charge transfer effect in the Fe doped LiNbO3. However, due to the well known problem of underestimation band gap using DFT, the results have to be treated with caution. A detailed discussi on is given in section. 4.6.2. Table 4 7 Charge transfer level of impurities calculated by GGA with comparison of experiments. All units in eV Fe 2+ /Fe 3+ Mg 1+ / Mg 2+ Nd 2+ / Nd 3+ Er 2+ / Er 3+ E A 1.84 3.42 3.38 2.75 E g 3.5 3.50 3.50 3.50 E PL 1.34 0.06 0.12 0.75 E rel 0.33 0.0 4 0.01 0.00 Experiment 1.4 } 182 /1.6 183 No Associated Peak 184 No Associated Peak 185 0.80 186 Using a similar approach, the charge transfer levels are calculated for Mg, Nd, and Er. The results of the cal culations are given in Table 4 7 It is noted that for both Mg and Nd, the thermal ionization energies are relatively large and close to the calculated band gap, which result in the EPL being very small. T hese small values of EPL correspond to very long wavelength s, beyond those experimentally accesible Consistent with this, no associated peak has been experi mentally reported. The Er charge transfer levels are an interesting case It is found that the calculated value of EPL is 0.75 eV while there is a experimental peak at 0.80 eV, which is equivalent to 1540 nm 186, 187. However, the calculati on and experiments are believed to describe different process. The PAGE 117 117 experimental peak at 1540 nm has been concluded to be the transition from 4I15/2 to the ground state 4I13/2 These two states differ only in the total spin number. Therefore, it can be con sidered as the local rearrangement of the electrons, which does not involve any charge transfer effects. On the other hand, the calculated EPL simulate a process in which electrons in the condution band jump into the charge transfer level. Therefore, in or der for this process to happen, it is necessary for condu c tion band to have electrons. For Fe doped, there are both Fe2+ and Fe3+ in the materials 182, 183, therefore, there are electrons at the charge transfer level, which is relatively close to the condu c tion band and which can be excited easily. For Er charge transfer, there is no Er2+ in the sytem Thus the incident light energy should have an energy equal to or higher than the band gap to exite electron from the valence band to condu c tion to condution band. However, the experimental incident light usually has a wave length longer than 440nm 182, 186, which is less than the band gap. Thus, there are no electrons in the condution band and the charge transfer level can not be observed from experiments. 4.6 Discussions 4.6.1 Er Doped LiNbO3 The principal conclusion s of this work are that (i ) the Er dopant sits on the Li sites, (ii) various slightly different microstructural configurations of Er sites exist 143, and (ii i ) charge compensation is provided by Li vacancies for Li2O rich conditions or (iv ) by Er on Nb sites under Nb2O5rich conditions. We briefly examine the experimental evidence relevant to each of these predictions. (i ) XSW studies 39, 40 found that Er occupies a Li site with C3 point symmetry, but that it is shifted from the ferroelectric Li position by 0.46 alon g the uniaxial direction of the hexagonal structure. ESR 4 identified a single site for Er3+ ions, but could not distinguish whether it is Li site or Nb site. Rutherford backscattering 38 and EXAFS 41 also yielded Er3+ substituting PAGE 118 118 for Li ions, while ion beam channeling 152 indicate that Er may occupy both Li and Nb sites. The DFT calculation s predicted that Er primarily sits on the Li site under both Li2O rich and Nb2O5rich conditions, but with different concentrations based on the difference values of the DFEs. The site is predicted to be shifted along the Z di rection by 0.17 This is consistent with the experimental data measured using XSW and EXAFS which show Er occupies the Li site and moves along Z direction. However, the magnitude of displacement is quite different. Possible reasons for this is discussed in section 4.4. The case of Er site on both Li and Nb site, which is relevant to ion beam channeling, is discussed in (iv) (ii) Optical and magnetic resonance spectroscopy revealed multiple slightly different sites for the Er ion. 143, 188 Baumann et al. using XSW and optical spectroscopy 43 found that there are four energetically distinguishable sites for Er, resulting from different local environments. Gill et al. 45, 46 using total site selective spectroscopy reported six different sites. Furthermore, Dierolf et al. ,143 by combini ng excitation emission as well as absorption, emission, and Raman spectroscopy (CEES), found a large number (>11) of different Er3+ sites. Moreover, the results of Dierolf et al. 143 suggest that the charge compensation can be classified into two categories : a first close compensation which influences the local symmetry of the Er defect and a secondary probably somewhat further away which influences the local electric field along t he z axis. The inclusion of a second more distant charge compensating Li vacancy further increases the number of potential sites as well. These different arrangements are associated with different local and distant charge compensation species A2 arrangements (see Sec. 3.B) are considered as the first category while A1 and A0 arrangements fall into the second category. The experimental structure analysis supports there are many slightly different sites for Er3+. When considering only lithium vacancies aroun d the Er3+ site, seven different configurations could be identified using PAGE 119 119 DFT calculations. If the presence of impurities, niobium antisites, Er Er cluster or even domain walls were to be considered, the number of different Er3+ site would be considerably larger. The optical spectroscopy measurements show that the relative number of differently charge compensated Er ions depends on the stoichiometry of the sample. While in congruent samples the more perturbed A2 type arrangements dominate, in stoichiometric samples the less perturbed arrangements are more prominent. Furthermore, under domain inversion at room temperature during which the charge compensating defects are immobile, drastic changes are observed for some defect configurations and no changes are observed for others. In particular, the majority defect configuration in stoichiometric samples remains unchanged (aside from a small spectral shift). For this reason, it is likely that this defect has distant charge compensation as described by configurati on A0. (iii) A consensus as to the precise nature of the charge compensation mechanism for ErLi .. has not yet been achieved. It has been reported by Gil l et a l. 46 that the site redistribution a ssociated with an increase in Er concentration is very similar to that associated with an increase in Li deficiency, which indicates that ErLi .. is compensated by lithium vacancies 143, 189. This is based on the fact that for congruent LiN b O3, lithium vacancies already exist in the material 136, 190. Therefore, it is much easier for Er3+ to occupy a vacant site. The DFT calculations predicted that under N b2O5rich conditions, which includes ambient temperature and oxygen partial pressure, the dominant defect are ErLi .. compensated by VL i based on the formation energies, which provides theoretical support to this conjecture. It could also explain the observation that increasing the Er concentration tends to increase the Li deficiency in the crystal. (iv ) No direct observation for ErNb has bee n reported, although ion beam channeling indicates that Er may occupy a Nb site 152 and electron paramagnetic resonance (EPR) indicates PAGE 120 120 Er sits on both Li and Nb site 191. This absence may be partially due to the methods to incorporate the Er ions into the crys tal. The Er3+ doping is often achieved during Czochralski growth or diffusiondoping 43, the latter involving the diffusion of Er ions into congruent LiNbO3. It has been reported that the diffusion of Er3+ is through the Li sublattice 142, 192. The diffusion of Nb was also found to be much slower than that of lithium vacancies 142, 193, 194. Therefore, it is probably much easier for a Er ion to move into a Li sublattice and stay there rather than to push Nb ions off their lattice sites during these kinetic processes; this would make the incorporation of Er on a Nb site difficult. Furthermore, both Czoc hralski growth and diffusion doping predetermine the dominant defects to Er on Li site compensated by lithium vacancies. These existing defects can be expected to play a significant role when the sample is converted to the stoichiometric composition using VTE under Li2O rich conditions. Taking nearly stoichiometric Li0.498Nb0.5 + xO3 for example, after VTE process the concentration of lithium vacancy is still 0.4 %. If we assume that the lithium vacancies are primarily compensated by ErLi .., 0.2 % of ErLi .. would be require d. The system would thus be equivalent to 0.2 mol % Er:LiNbO3. Since the concentration of ErLi .. and ErNb is a constant, the concentration ErNb of can be calculated from the mass action law as: ) / exp(' 'kT G Er Erf Nb Li [4 3 ] where, Gf is the defect formation energy, k is the Boltzmann constant, and T is the temperature in Kelvin. At the room temperature, the concentration of ErNb is estimated to be 9.7 x 10 26, which are ~23 orders of magnitude smaller than ErLi ... Increas ing the temperature will increase the concentration of ErNb such that at 800 K, the concentration of ErNb is calculated to be 1.01 x 10 5 However, this is still two orders of magnitude smaller than the concentration of PAGE 121 121 ErLi .. Thus, c ompared to the concentration of ErLi .. t he ErNb concentration is negligible. These calculations indicate that even a smal l lithium vacancy concentration could significantly affect the defect chemistry of the system. Furthermore, it may provide an expla nation why the Er on Nb site has not been observed. 4.6.2 Fe Doped LiNbO3 The key findings about Fe doped LiNbO3 from DFT calculations are the following. (i) Fe prefer s to occupy a lithium site, compensated by lithium vacancies for the Nb2O5 reference state and by Fe on a niobium site for the Li2O reference state; (ii) Fe3+ is more energetically favorable than Fe2+; (iii) There is an charge transfer level provided by iron at the position of 1.4 eV below the conduction band. (i ) In spite of exten sive effect s to determine the occupied site of iron using various techniques, such as optical absorption 195, electron spin resonance (EPR) 196, Mossbauer spectroscopy 36, 37, and Proton Induced X ray Emission (PIXE) 42, no experimental consensus has yet emerged Most of the st udies suggest that iron occupies a lithium site, which is supported by electron nuclear double resonance (ENDOR) 35, PIXE 42, X ray standing waves 197, and X ray absorption fine structure analysis 198. On the contrary several studies using EPR and ESR claim that iron occupies the niobium site 196, 197. In addition there are some studies that support two site model in which iron o ccupies both lithium and niobium site s 199. Our calculations have shown that iron on lithium site is favorable under both Nb2O5 and Li2O rich conditions, but with different charge compensation mechanisms. Under Nb2O5 reference state, in which the intrinsic defects are predicted to be niobium antisites compensated by lithium vacancies and a congruent composition is formed the compensating defects are lithium vacancies. This supports the studies using ENDOR, PIXE, Xray related techniques 35, 42, 197. Under Li2O rich conditions, which PAGE 122 122 correspond to growth condition of stoichiometric LiNbO3, the compensation mechanism is considered as Fe on a niobium site. However, since the dopants are usually introduced into system through a di ffusion process on a congruent sample, it is very difficult to observe iron on a niobium site due to existing intrinsic defects Similar discussions about how the growth method, impurities introduction process, and the intrinsic defects may affect the domi nant defects have been given in section 4.6.1 Overall the DFT calculations indicate that by controlling the experimental condition or the composition of the material, the site occupancy of the iron can be altered (ii) As a transition metal, iron ions have d ifferent valence charges (Fe2+ and Fe3+), which is also the case in Fe doped LiNbO3. In most of the studies, the concentration of iron is very low 182, 183. Furthermore, the concentration of Fe2+ is usually at least one order of magnitude lower than F e3+ concentration 182. Although the concentration ratio of Fe2+/Fe3 can be varied using thermal treatments or other types of technique 182, 200, the naturally occurred concentration ratio of Fe2+ and Fe3+ indicate s that Fe3+ is more stable DFT calculations show that the defect reactions involved with Fe3+ have lower formation energies than the reactions involve with Fe2+, which corroborates experimental observation. (iii) The Fe2+/Fe3+ charge transfer level is then calculated to be1.34 eV below the conduction band using the calculated band gap of 3.5 eV. (Figure 4 14) However, a s indicated previously, the band gap is 0.28 eV lower than the experimental value. The difference also represents the possible error bar in t he charge transfer level. Using the experimental value of the band gap(3.78 eV 115, 116), the charge transfer level will be 1.62 eV. The experimental photoluminescence measurement of peak associ ated with iron is determined to be 1.4 eV at temperatures below 300K and for various iron concentrations 182. With treatment at 520 oC in PAGE 123 123 lithium carbonate, a luminescence ba nd at 770 nm (which is equivalent to 1.7 eV) due to iron center has been observed 183. The difference between the experiments also indicates the difficulty of determining this value. However, both values fall into the range predicted by DFT calculations. Therefore, the DFT calculations confirm the longt ime conjecture that the peak is associated with charge transfer between Fe2+ and Fe3+. Figure 4 1 4 The relative position of Fe2+/Fe3+ charge transfer level with regard to the band gap. PAGE 124 124 CHAPTER 5 POINT DEFECTS IN CERIA BASED ELECTROLYTES 5 .1 Solid Oxide Fuel Cell Electrolytes The first test of a solid oxide electrolyte for solid oxide fuel cells (SOFCs) was carried out in the 1930s Since then, research has focused on finding an appropriate electrolyte. To be used in practical applications, solid el ectrolytes need to satisfy several major requirements, which include a high ionic conduction, negligible electronic conduction, and thermodynamic stability over a wide range of temperature and oxygen partial pressure. In addition, solid electrolytes must h ave a thermal expansion compatibility with electrodes and suitable mechanical properties 50. Among various electrolytes, yttria stabilized zirconia (YSZ) 5154 has been most extensively investigated and is widely used Compared with other solid electrolytes, YSZ exhibits a very low electronic contribution to total conductivity in the oxygen partial pressure range most important for practical applications 50. However, YSZ electrolytes require a very high operating temperature (~1000 oC), which could cause many interface reactions such as the interaction between solid electrolytes and electrodes Thus, research focus on other potential electrolyte materials is of great interest. As an alternative to YSZ, doped ceria 5659 h as been the subject of a great deal attention due to its higher ionic co nductivity than that of YSZ. The electrical conductivity of Ce0.80Gd0.20O1.9 (gadolinia stabilized ceria) is about one order of magnitude higher than that of YSZ in the intermediate temperature (IT) range (500oC700oC) (Figure 5 1). Although Bi2O3based ma terials showed even higher ionic conductivity than YSZ and doped ceria, these materials are unstable at high temperatures and low partial pressures 60, which have severely limited their applications. PAGE 125 125 Figure 5 1. Electrical conductivity of v arious electrolytes as a functional of temperature 201. The higher ionic conductivity of ceria based electrolytes with respect to YSZ is the main advantage for their use in SOFCs. However the doped ceria systems undergo reduction under low oxygen partial pressures, and the magnitude of electrical conductivity and stability are greatly dependent on the type and concentration of dopants 55. Overall, ceria based elec trolyte shows considerable potential for practical application. 5 .2 Structure of Ceria Based Materials Ceria has a fluorite structure (space group: Fm 3m Figure 5 2 ) which c an be considered as an anion cube cage encapsulated in a face centered cubic (FC C) cation lattice (Fig 5 2 ). In the lattice, Ce atoms sit at the corner of the cubic lattice, while the O atoms sit at the tetrahedral sites. This indicates that each Ce atom is coordinated with eight O atoms and each O atom is coordinated with four Ce ato ms. From the radius ratio calculation r Cation / r Anion = 0.8347 (crystal radius for cerium and oxygen are 1.01 and 1.21 respectively 202), which falls in the 0.7321 PAGE 126 126 regime. Therefore, the cubic arrangement of the O atoms around the Ce is consistent with t he model established by Pauling 203. Figure 5 2. Schematics of fluorite structure for CeO2. Another instructive way of describing the crystal structure is by looking at the stacking sequence through different plane s. Fig ure 5 3 shows the Kagome 204, 205 net formed by Ce atoms. As the Ce atoms forming the net are in different planes, a shift vector is needed to relate the projected positions of Ce atoms in different planes. Thus, from the atomic positions, one can identify Ce atoms belonging to different planes in the net. The c eria unit cell contains four ceria molecules (Z = 4) Considering the cation as refere nce, the location of the atoms, site symmetry, multiplicity and Wyckoff letter in an ideal fluorite structure are given in the Table 5 1. In reality, the atomic positions may vary a little from the ideal lattice position. A n eutron diffraction study on cer ia by Yashima et al showed that the PAGE 127 127 O atoms are distributed at 8c and 32f positions 206. The structural details of the neutron diffraction study are summarized in Table 5 2. Figure 5 3. Kag ome net formed by Ce from <1 1 1> view Only the Ce sublattice is shown. Table 5 1 Potential Fluorite structure data from International Tables of Crystallography, Volume A: space Group Symmetry, Fourth Edition 207. Atom Site Site symmetry Coordinates O 8c 43m 1/4, 1/4, 1/4 1/4, 1/4, 3/4 4b m3m 1/2, 1/2, 1/2 Ce 4a m3m 0.0 0.0 0.0 Table 5 2 Atomic position of ceria from neutron diffraction 206. Atom Site U( 2 ) x y z Ce 4a 0.0340 0.0 0.0 0.0 O 1 8c 0.041 0.25 0.25 0.25 O 2 32f 0.041 0.325 0.325 0.325 Doped ceria has a very open structure and it exhibits large tolerance for high levels of atomic disorder. Oxygen vacancies c an be introduced by doping with oxides of metals with low valences. For instance the elements with 2+ or 3+ valence as MO and M2O3: x o O CeO V M MO [5 1] PAGE 128 128 x o O CeO V M O M 3 2' 3 2 [5 2] Many studies suggest that vacancies introduced by doping are bound to dopant cations to form some complex defect associates 55, 208. This association energy is m ainly due to the electrostatic interaction of the defects caused by their charge The association energy is also dependent on the ionic radii and dielectric polarizability of the doped cations which may cause elastic strain in the system 58. Because of the presence of oxygen vacancies in the doped ceria structure and because of the role of defect in determining physical properties of ceramics, the mechanical properties are quite different from that of the perfect structure. This is supported by the study of A t kinson et al 209 using a simulation approach, which showed that ceria based electrolytes will be mechanically unstable above 700oC. This indicates that the mechanical property of doped ceria is inferior to YSZ 52. However, the origins of this are not fully understood. 5 .3 Computational Details Molecular dynamics simulations and lattice statics calculations are carried out in the present study. The interatomic interactions between cations and anions are described by combining both longrange Coulombic interaction and short range repulsive contributions. The long range Coulombic interactions are : N i i j ij j i ij LRr q q r E12 1 [5 3 ] where N is the total number of ions in the system, qi and qj are the charge on ion i and j respectively, rij is the distance between ion i and j The summations over the system are carried out using the Ewald method. All of the potentials used in the present study employ the formal ionic charge for cerium and oxygen (i.e., qCe = +4 and qO = 2). PAGE 129 129 S hort ranged interactions are typically described by one of two types o f potentials. The first one is the Buckingham potential : 6 ) / (/ij r ij Buckr C Ae r Eij [5 4 ] where and C are adjustable parameters, chosen to reproduce pertinent physical properties of the real material. The other one used is the Born Mayer Huggins (BMH) potential: 6 0/ )} /( ) exp{( ) (ij j i j i ij j i j i ij BMHr c c b b r a a b b f r E [5 5 ] Here ai, bi and ci are parameters for individual atom species. Although the BMH form can be simply recast into the Buckingham form, it does have the advantage that, at least in principle, parameters for indi vidual ions can be defined, rather than for ion pairs as in the Buckingham form. For ceria based systems, parameters in both the Buckingham and BMH forms have been developed. Since the quality of the simulations is directly determined by the quality of the potential, a n evaluation of the available potentials is given in Section 5 .4. The calculations used in the evaluation of the potentials are carried out using General Utility Lattice Program (GULP) 210, 211. All simulations a re performed on a 6 x 6 x 6 supercell of cubic non primitive fluorite unit cells, each of which contains four CeO2 formula units, for a total of 2592 ions. Periodic boundary conditions are applied in all three spatial directions. To simulate sub stoichiometry, Ce4+ ions are randomly replaced by Ce3+ (q=+3) ions. Within a single simulation th e arrangement of Ce4+ and C e3+ ions is fixed; this is a limitation of the current approach, since experimentally electron transfer allows the Ce3+ ions to develop an equilibrium arrangement, which may not be random. An app ropriate number of oxygen ions a re then removed to maintain c harge neutrality. Computational annealing of the system by MD at high temperature allows the oxygen vacancies to diffuse to form a structure in which the anion and PAGE 130 130 vacancy arrangements are in equilibrium with respect to the fixed cation sublattice. For sub stoichiometric ceria and doped ceria, since the properties may depend on the arrangements of dopant cations and oxygen vacancies, several random cation structures are generated and the properties of each structure are calculated to capture the range of ef fects. For T > 0 K calculations, the free energy is calculated using the quasi harmonic approximation to lattice dynamics, as implemented in the GULP code. 5 .4 Interatomic Potentials Evaluation F ive different parameterizations of the Buckingham potential and one parameterization of the BMH potentials are currently avail able for the ceria based systems. The parameters for the five Buckingham potentials, denoted as Grimes 212, Vyas 213, Butler 214, Gotte (2004) 215, and Gotte (2007) 216, are given in Table 5 3. The Grimes, Gotte (2004) and Gotte (2007) potentials have also been parameterized for Ce3+, allowing the effects of off stoichiometry to be assessed. In addition, the Gr imes potential has been parameterized for a number of rare earth ions, allowing doping effects to be determined. The parameters for Inabas BMH potential are given in Table 5 4. All of the Buckingham potentials have shell model parameters; unless otherwis e indicated all of the simulation results for these systems are based on shell model calculations. Inabas BMH potential is a rigid ion potential. T o evalua te the quality of each potential, the lattice constant, thermal expansion, chemical expansion, oxyge n migration energ y and dielectric properties are determined. The detailed comparison for these properties is given in the following subsections 5.4.1 Lattice Constant and Thermal Expansion The zero temperature lattice constants calculated using these po tentials are listed in Table 5 5. The values are all essentially identical and in excellent agreement with the experimental data This is to be expected since the potentials were fitted to the lattice parameter and some PAGE 131 131 other fundamental properties. Howev er, it should be noted that the experimental data is 293 K while simulation data are for 0 K. Table 5 3 Parameters of the five Buckingham potentials for ceria based systems. Species A (eV) r () C (eV 6 ) Y e (e) K 1 ( eV 2 ) Ref Grimes O 2 O 2 9547.96 0.2192 32.0 2.04 6 .3 217 218 Ce 4+ O 2 1809.68 0.3547 20.40 0.20 177.84 213 Ce 3+ O 2 2010.18 0.3449 23.11 58 In 3+ O 2 1495.65 0.3327 4.33 219 Y 3+ O 2 1766.40 0.33849 19.43 217 Gd 3+ O 2 1885.75 0.3399 20.34 58 La 3+ O 2 2088.79 0.3460 23.25 58 Vyas O 2 O 2 9547.9 2 0.2192 32.0 2.04 10.3 213 Ce 4+ O 2 2531.5 0.335 20.40 0.20 177.84 Butler O 2 O 2 22764.3 0.149 45.83 6.10 419.9 214 Ce 4+ O 2 1986.8 0.3511 20.40 7.7 291.8 Gotte ( 2004 ) O 2 O 2 9547.92 0.2192 32.0 60.78 3.0 215 Ce 4+ O 2 1809.68 0.3547 2 4 .40 166.021 7.0 Ce 3+ O 2 1809.68 0.3547 2 4 .40 166.021 7.0 Gotte ( 2007 ) O 2 O 2 9533.421 0.234 224.88 1759.8 6.5667 216 Ce 4+ O 2 755.1311 0.429 0.0 43.451 4.6475 Ce 3+ O 2 1140.193 0.386 0.0 2172.5 15.092 Table 5 4 Parameters of the Inaba 220 Born Huggins Mayer potential for CeO2. No shell model is available for this potential. z a () b () c ( J 0.5 (nm) 3 mol 0.5 ) Inaba Ce 4+ 2.700 1.33 0.0454 0.00 O 2 1.350 1.847 0.166 1.294 To determine the thermal expansion, the temperature dependence of the lattice parameter for each potential (normalized to the corresponding T = 0 K values) is shown in Figure 5 4 The quasi harmonic approximation used for these calculations includes quantum effects, which leads to zero thermal expansion at T = 0 K. The coefficient of thermal expansion values given in Table 5 5 are taken from T= 4 00 K to T= 14 00 K over which the lattice parameter is approximately linear in temperature and over which quantum effects are no longer important. None of these potentials is able to reproduce the experimental value quantitatively It is found that the Butler potential overestimate s the experimental thermal expansion whereas the others underestimate the thermal expansion. The thermal expansion calculated by the But l er and PAGE 132 132 Gotte (2007) potentials which are 112% and 82% of the experimental value, are better than other values. The thermal expansion calculated using Gotte (2007), Gotte (2004), Vyas and Grimes potential s predict values that are 63%, 57%, 10% and 11% o f the experimental value, respectively. Figure 5 4. Temperature dependence of the lattice parameter of CeO2 obtained from MD simulation using different potentials compared with experiment data 221. 5 .4.2 Chemical Expansion The chemical expansion measures the lattice expansion produced by off stoichiometry. To simulate this, we introduce d oxygen vacancies into the pure ceria system and replace d the requisite number of Ce4+ ions with Ce3+ ions. The effect of this sub stoichiometry on the lattice parameter is four fold (i) T he larger i onic radius of the Ce3+ compared to the Ce4 + ions ( 1.143 vs. 0.97 ) tends to increase the lattice parameter (ii) T he lower charge of the Ce3+ reduces the strength of the attractive Coulombic interactions, thereby also tending to increase the lattice PAGE 133 133 parameter (iii) the presence of oxygen vacancies decreases the attractive Coulombic interactions, also tending to an increased lattice parameter (iv) T he oxygen vacancies also have the effect of lowering the atomic density thereby tend ing to decrease the lattice parameter. We note that three of these four effects can be expected to lead to an increase in the lattice parameter with increasing sub stoichiometry. Figure 5 5. The effect of composition on the lattice constants of CeO2x obt ained from MD simulations at 1073K. The values are compared with experimental results 221, 222. As mentioned above in section 5.4 only three of the potentials have parameters for Ce3+: Grimes, Gotte (2004) and Gotte (2007). Figure 5 5 compares the simulated lattice parameters at 1 0 73 K using the Grimes Gotte (2 004) and Gotte (2007) potential s with the experimental value for CeO2 x (x = 0.0 0. 15) 222, 223. The simulation results at each temperature show that the lattice parameter increases with increasing vacancy concentration consistent with the qualitative argument above The slope of the chemical expansion line with respect to x in CeO2 x is given in PAGE 134 134 Table 5 5 It is found that the both Grimes and Gotte (2007) potential give a chemical expansion close to the experimental results which indicates that both potentials successfully capture the effects of sub stoichiometry. Table 5 5 Properties of CeO2 from experiment and from various interatomic potentials. Exp. Empirical Potential Grimes 58 Vyas 213 Butler 214 Got te (2004) 215 Gotte (2007) 216 Inaba 220 Lattice Constant () 5.411 5.411 5.414 5.411 5.411 5.411 5.414 Coefficient of Thermal Expansion (X 106 K1) 11.6 1.27 1.16 13.1 6.65 7.31 9.58 Coefficient of Chemical Expansion @ 1173K (/Oxygen Vacancy Concentration) 0.535 0.41 0.505 0.43 Static Relative Permittivity 18.6 20.0 18.6 12.71 19.53 16.2 18.6 High Frequency Dielectric 4 4.54 3. 98 4 5.59 4 Oxygen Migration Energy (eV) 0.70.1 0.3 (shell) 0.63 (rigid ion) 0.54 0.74 0.41(shell) 0.58(rigid ion) 0.58 0.54 5 .4.3 Dielectric Properties and Oxygen Migration Energy As shown in Table 5 5, all the potenti als reasonably reproduc e both the static and highfrequency dielectric constants The oxygen migration energy predicted by Grimes and Gotte (2004) are significant ly smaller than the experimental value. Interestingly, using the Grimes potential in the rigid ion approximation, rather than with a shell model, the predicted oxygen migration energy is in excellent agreement with the experimental results. The Gotte (2004) PAGE 135 135 potential also predicts higher migration energy when the shell model is not used. Physically, this is quite reasonable since the polarizability of the ions in the shell model allows charge s of ions to polarize in a manner that lowers the energy barriers. The oxygen migration energy calculated using Vyas, Butler, Gotte (2007) and Inaba are also in reasonably good agreement with experimental data. 5 .4.4 Elastic Constants Table 5 6 compares the zero temperature elastic constants as determined directly fro m GULP for the various potentials, with the experimental results and with DFT. The e xperimental values for the elastic prop erties are also given in Table 5 6 The values of C11, C12 and C44 we re obtained from Raman scattering, which probes the fundamental vibrations of the system 224. The value of the bulk modulus determined from these values is in quite good agreeme nt with independent determinations from pressure volume curves 225. Likewi11se, the Youngs modulus obtained f rom analysis of the Raman scattering results is in good agreement with the value from nanoindentation 226. For the simulations, the values of C11, C12 and C44 were determined in GULP, from which all of the moduli and the anisotropic factor Z were determined using standard analytic elastic relationships. The values predicted by the Grimes, Vyas, Butler, Inaba and Gotte (2004) potentials are all larger than the corresponding experimental values. The values obtained for the Gotte (2007) are in essentially perfect agreement with the experimental values. Although not discussed in their paper, it can reasonably be assumed that the experimental elastic properties were part of the fitting set for the potential. The bulk modulus B, tetragonal shear modulus G, and Youngs modulus, Y, show corresponding levels of agreement with experiment. PAGE 136 136 Table 5 6 Comparison of elastic properties, including the bulk modulus, B, Young s modulus, Y, and shear modulus, G= (C11C12), of CeO2 from experiment, DFT calculations and as predicted by empirical potentials. Those experimental and simulated values that were independently determined are given in bold. Others were determined using standard relationships among elastic properties Exp. Buckingham BMH DFT All units in GPa Grimes 58 Vyas 213 Butler 214 Gotte (2004) 215 Gotte (2007) 216 Inaba 220 C 11 403.0 224 554.2 573.0 504.4 554.0 403.0 451.0 399.5 227 C 12 105.0 224 124.6 14 7.7 143.1 125.0 105.0 103.6 127.5 227 C 44 60.0 224 123.6 14 6.8 16.1 71.7 60.0 99.1 63.5 227 B 204 2364 225 23010 228 267.9 289.4 263.6 268.0 204 .3 207.0 218 G 149.0 215.0 212.6 180.7 214.5 149.0 173.7 136 Y 508.0 512.5 439.7 511.2 362.4 412.3 338 G Po ly 96.0 160.2 173.4 81.8 128.2 96.0 124.3 92.5 Y Po ly 248.95 255 226 400.7 433.6 222.4 331.7 249.0 310.6 243.1 Z 0.40 0.58 0.69 0.09 0.19 0.40 0.57 0.46 Also included in this table are the results of previous electronic structure calculations at the level density functional theory (DFT) using the local density approximation (LDA) with the projector augmented wave (PAW) method 62. Both the atomistic simulations and the electronic structure calculations were performed on ideal crystals. We therefore used the Voigt Reuss Hess method 229 to determine polycrystalline average elastic constants from the single crystal values. Because C11 = C33 and C12 = C13 = C23 PAGE 137 137 for a cubic system, the polycrystalline average of the bulk modulus is identical to the single crystal value. The H ess average is shown in Table 5 6 The polycrystalline averaged shear modulus li es within the bounds 229: 5 3 3 4 544 12 11 12 11 44 44 12 11C C C GG C C C C C C Gv R [5 6 ] The elastic anisotropy of the material can be characterize d by the directionality of the Youngs modulus, E, can also be determined using: 2 sin sin 4 1 cos sin 1 22 4 2 2 44 12 11 12 11 12 11 1 Z C C C C C C C E [5 7 ] where between [010] axis and the projection of the chosen direction on the (100) plane. Here Z is the Zener anisotropy ratio: 12 11 442 C C C Z [5 8 ] For Z=1, the system is isotropic, and the Youngs modulus is independent of direction. For Ceria, the Yongs modulus along different crystallographic directions is shown in Table 5 7. Since all of these metrics of the elastic properties were determined analytically from the elastic constants, they show corresponding degrees of agreement with experimental values. Table 5 7 Anisotropy in the Youngs modulus of CeO2 using different methods. Numbers in brackets indicate fractional value relative to that in the <001> direction. Directions Grimes 58 Gotte(2007) 216 DFT 227 Exp. <0 0 1> 504.9 (1.0) 359.6 (1.0) 355.9 (1.0) 359.6 (1.0) <1 1 0> 351.6 (0.70) 189.8 (0.53) 196.7 (0.55) 189.8 (0.53) <1 1 1> 319.3 (0.63) 163.9 (0.46) 171.2 (0.48) 163.9 (0.46) 5 .4.5 Potential Summary None of the potentials is able to reproduce all of the experimental properties evaluated. Overall, the Gotte (2007) has the highest materials fidelity. In particular, the elastic and dielectric PAGE 138 138 properties are in essentially perfect agreement with experiment, and the predicted oxygen migration energy is in good agreement with experiment. The chemical expansion is also in reasonable agreement with the experimental data. However, the predicted thermal expansion is smaller than the experimental value, i ndicating that the potential is significantly less anharmonic than the experimental system. Interestingly, the Grimes potential predicts the chemical expansion to within 20%. It also yields dielectric properties and the oxygen migration energy in good agr eement with experiment. However, it does a very poor job in predicting the coefficient of thermal expansion, predicting a value that is only 11% of the experimental value. In terms of the elastic properties, the Grimes potential is clearly not as good as t he Gotte (2007) or Inaba potentials, but has a similar level of fidelity as the other potentials. To its advantage, the Grimes potential is the only potentials providing parameters for some rare earth ions, which allows the study of doped ceria system. In the next section, we evaluate the effects of off stoichiometry on the elastic properties. Because both the Grimes and Gotte (2007) potentials do predict some of the salient properties with reasonable accuracy, we will compare the results predicted for eac h. In Section 5 .5 we examine the effects of rare earth doping on the elastic properties, for which the Grimes potential is the only one with potentials parameters. 5 .5 Mechanical Softening Effects A key finding in the experimental measurements of Wang et al. 226 for the Youngs modulus of ceria was that there is significant elastic softening with decreasing oxygen partial pressure (PO2), i.e., with increasing sub stoichiometry 226. The open square symbols in Fig ure 5 6 denote the room temperature Youngs modulus for different levels of substoichiometry determined from their data The experiments were actually performed at fixed PO2, at which the la ttice parameter was determined. The corresponding degree of off stoichiometry can be PAGE 139 139 calculated from the method developed by Bevan 230 and Panlener 231. Figure 5 6 also shows the corresponding results predicted from the Grimes and Gotte (2007) potentials at 300 K. In Figure 5 6, the Youngs moduli are plotted as a function of lattice constants rather than of the composition. There are dif ferences between the predicted chemical expansions and the experimental value; it thus turns out that the direct calculated Youngs moduli agree with the experimental values quite well for the Gotte (2007) potential The Grimes potential agrees less well, mainly due the significant over estimate in the elastic constant s of the pure system. Figure 5 6. Composition dependence on the room temperature Youngs modulus of CeO2 x obtained from MD simulation using Gotte (2007) potential, Grimes potential and from experiment 226. In contrast to the Raman scattering experiments the DFT calculations and the MD simulations, the nanoindentatio n experiments of Wang et al. 226 did not show significant elastic anisotropy. It is possible that microstructural factor in the indented samples le d to effectively PAGE 140 140 isotropic elastic behavior, rather tha n the anisotropic behavior seen in the spectroscopy experiments 224, which probe the intrinsic atomic level elastic response. We have also used atomic simulations to determine the effect of temperature on the elastic properties of CeO2 using the Gotte (2007) potential, which has the highest fidelity for elastic properties. Figure 5 7 shows the tempera ture dependence of Youngs modulus (Y) in the <001> direction the shear modulus (G), the bulk modulus (B), and the anisotropy factor (Z) The trends are in accord with conventional trends for the softening of the elastic properties with temperature. The elastic anisotropy factor also decreases with increasing temperature. Figure 5 7. The effects of temperature on the normalized bulk modulus (B) shear modulus (G), Youngs modulus (Y), and anisotropic factor (Z) for CeO2, as determined by MD simulation using Gotte (2007) potential We have determined the effects of doping of ceria doped with In3+, Y3+, Gd3+, La3+ ions, for which the potential parameters compatible with the Grimes potential were provided by Minervini et al. 58. The ionic radii o f these dopant ions are given in Table 5 8. Even though t he PAGE 141 141 Grimes potential cannot accurately reproduce the thermal expansion, dielectric constant, and oxyg en migration energy accurately; it does accurately reproduce the most two important properties for addressing the effects of doping are the chemical expansion, for which it does a good job, and the elastic constants, for which it captures the trends, albeit being overall too stiff. Figure 5 8. The effect of dopa nt ion size on the lattice parameter of Ce0.8M0.2O1.9 (where M = In3+, Y3+, Gd3+, Ce3+, La3+) obtained from MD simulation, compared with the experimental results 232. Table 5 8 Ionic Radii for relevant ions from Shannon Table 154. Dopant Ions Ionic Radius (VIII) Ce 4+ 0.97 In 3+ 0.92 Y 3+ 1.019 Gd 3+ 1.053 Ce 3+ 1.143 La 3+ 1.16 PAGE 142 142 We compared the lattice constants of various doped ceria systems with the experimental values. In a manner analogous to that for CeO2x, the doped structures are generated by randomly replacing the Ce4+ ions by dopant cations and compensating the charge by creating appropriate number of oxygen vacancies in the system. Since the experimental results for dopants are already available at 298 K and the effect is more pronounced for a larger dopant concentration, room temperature Ce0.8M0.2O1.9 was simulated, where M = In3+, Y3+, Gd3+, La3+ ions. Figure 5 8 compares the simulated lattice parameters using the Grimes potential with the experimental values 232. The simulated values are within 0.2 of the experimental values over the full range of ionic radii 154, which can be considered to be in a adequate agreement, As is evident from Figure 5 8, the simulations show a consistently weaker dependence of the lattice parameter on the ionic radius of the dopant ion. This is consistent with the calculated lower chemical expansion in CeO2x (see Table 5 5). For fuel cell electrolyte applications, it is also important to understand the effect of dopants on the Youngs modulus of doped ceria system. From the point of view of the atomic level simulation the approach is essentially identical to that for sub stoichiometric CeO2 since the Ce3+ ions and M3+ ions act in the same manner. Experimentally, the Youngs modulus of Ce0.8Gd0.2O1.9 at room temperature is determined to be 147.2 16.1 GPa 233, whereas the Grimes potential predict twice this value. This indicates that Grimes potential overestimates the Youngs modulus of pure and reduced ceria system by simil ar amounts. However, the relative trend of Youngs modulus dependence on ionic radii is considered should be more reliable than the absolute value s As Figure 5 9 shows the system becomes elastically softer with increasing ionic radius of the dopant. Simi lar elastic softening effects have been observed experimentally in doped ceria systems 226. PAGE 143 143 Figure 5 9. Dopant size dependence on the Youngs modulus of Ce0.8M0.2O1.9 (where M = In3+, Y3+, Gd3+, Ce3+, La3+) obtained using Grimes Potential. 5 .6 Oxygen Diffusion T he diffusivity of doped ceria system is calculated based on the positions of oxygen ions at different time in the MD simulation. The mean square displacement (MSD) of cations and anions in the s ystem can be extracted from the position of oxygen ions in the usual manner as : 20 r t r MSD [5 9 ] where r (t ) r (0 ) is the distance traveled by an oxygen ion over time t To calculate the diffusivity, t his quantity needs to be averaged over all the oxygen ions in the system for a long period. The n, the diffusivity ( D ) can be obtained from 20 6 1 lim r t r t Dt [5 10] PAGE 144 144 After determin ing the diffusivity at different temperatures, the activation energy are calculated from the Arrhenius equation : ) / exp(0RT E D Da [5 11] where D0 is the pre factor and Ea is the activation energy. The conductivity can be then c alculated from diffusivity by Ne rnst Einstein equation 234 and compared with experiment data. For doped ceria systems, the diffusivities at different temperatures are given in Fig ure 5 10. From Equation 5 11, the activation energy of each system is calculated and compared with DFT results and experimental value (Figure 5 11). It can be seen that DFT predict Sm3+ and Pm3+ have the lowest activation energies, unlike experimental resu lts. The values predicted by DFT are also lower than experiments. It is noted that Gd doped ceria has the lowest activation energy among the trivalent doped systems that have been studied in the present paper, which is consistent with the experimental obse rvation. The underlying reason that Gddoped ceria is energetically favorable is attributed to the small mismatch of ionic radii between the host cations and dopant ions. Figure 5 10. Diffusivity dependence on temperature for various dopedceria systems PAGE 145 145 Figure 5 11. Activation energies of oxygen diffusion of various ceria based electrolytes. The results are compared with experimental values and DFT results. Since activation energy consists of the migration energy and association energy, to separate these two factors, the migration energies along different diffusion paths: <100>, <110>, and <111> have been calculated using the nudged elastic band (NEB) method (Figure 5 12). It is the observed that t he migration energy along <100> path is smaller than the other two, which indicated <100> is the favorable path. Furthermore, the migration energy barrier is calculated to 0.62 eV which is in excellent with experimental value of 0.5 0.6 eV 49. The association energy, which is due to the ela stic and electrostatic interactions between dopant and host cat ions, can be ca lculated by taking the difference between activation energy and migration energy. PAGE 146 146 Figure 5 12. Migration barrier along [001], [110], and [111] directions obtained from NEB calculations. 5 .7 Summary The materials fidelity of six different empirical potentials have been evaluated, based on lattice constants, thermal expansion, chemical expansion, dielectric properties, oxygen migration energy and mechanical properties. No s ingle potential can reproduce all fundamental properties. Therefore, choosing the appropriate potential(s) for specific applications is critical. Overall, the Gotte (2007) has the highest overall fidelity, while the Grimes can be applied to the widest range of systems. The simulations, using the best available potential(s) for each physical property, show that both substoichiometry and aliovalent doping of ceria lead s to large decreases in the elastic constants of the material. These decreases arise from the significantly reduced strength of ionic interactions. If this elastic softening were to be so large as to threaten the mechanical integrity of PAGE 147 147 the system, then it would pose a challenging problem for the design of ceria based electrolytes. It appears that the only way to mitigate this softening would be to strengthen the ionic interactions in the system. Of course, reducing the amount of doping or the substoichiometry would do this; however it would also significantly reduce the ionic conductivity of the system, which is the key performance metric of the material. It would thus be of significant value to identify a strategy for increasing the ionic conductivity that does not result in a significant softening in the elastic properties. PAGE 148 148 CHAPTER 6 CASCADE SIMULATIONS IN TITANIUM 6 .1 Background A fast reactor is a specific type of nuclear reactor which has no moderator and relies on fast neutrons alone to sustain the fission chain reaction 235, 236. Compared with thermal reactors, fast reactors produce more neutrons per fission, resulting in more neutrons than required to sustain the chain reaction. These extra neutrons can be used to produce extra fuels or to tra nsmute longhalftime waste to less troublesome isotopes 237. In fast reactors, the key component is the nuclear fuel which can be enriched uranium or plutonium In addition to nuclear fuel, there are many components in fast reactors, such as coolant, contr ol rods, cladding materials Cladding is t he outer layer of the fuel rods The cu rrent standard mater ial for cladding is zircal oy 238 239, which is corrosion resistant and is able to maintain the mechanical integrity du ring the operation. Zirca loy mainly consists of zirconium; other elements can be tin, niobium, iron, chromium, and nickel 240. Similar to zirconium, titanium also has a HCP structure, with even better corrosion resistance and mechanical properties. Therefore, titanium has been considered as a potential cladding m aterial for the next generation of fast reactors However, b efore any practical applications, it is of great importance to understand how radiation effects may affect its physical properties Radiation effect s in nuclear materials are of great complexity, as they are coupled with various microstructural phenomenon and processes. Therefore it is quite challenging to dis en tangle each process and to understand the corresponding mechanism. For instance, high energy radiation c an knock the atoms off their equi librium sites with extremely high kinet ic energy, generating an extremely high local temperature and pressure. Accordingly thermal and PAGE 149 149 pressure gradient s can be form ed. In addition the electrons can be stimulated from the ground state to an excited state initiating charge transfer process es or chemical reactions. All of these processes are entangled with each other and the complexity is overwhelming. Cascade simulations offer a powerful means of charactering the details of some of these sophisticated process and decoupling some of the phenomena 241, 242. The typical spatial and time resolut ions of cascade simulations match very well with the dimension s and lifetime of experimental c ascades by neutrons or heavy ions Moreover cascade simulations can provide information about production, evolution and subsequent diffusion of the defects informatio n which prove s to be of great importance in determining the microstructure evolution and mechanical properties of the materials. Based on previous studies of c ascade simul ations in a number of materials some consistent trends have been identified 243: The total number of stable surviving point defects is found to have a power law d ependence on the cascade energy; i t is also noted that the efficiency of generating stable defects decreases with cas cade energy until subcascade becomes prominent. In addition, the temperature dependence is rather weak. However, p revious cascade simulations of metallic materials mainly focus ed on cubic structure, such as iron 241, 244. Few studies are available for HCP structure 245, which has the potential to show anisotropy under the radiation. Therefore, it is of great importa nce to car ry out cascade simulations in H C P metals. In addition improvements, particularly with the accuracy of the in teratomic potential, are also needed for cascade simulation s In the present study, titanium have been carried out using high fidelity modified embedded atom method (MEAM) 83, 246 potentials combined with the potential developed by Ziegler, Biersack and Littmark (ZBL) 247 for short range interactions The PAGE 150 150 simul ations are performed with energies of primary knock on atom (PKA) ranging from 0.1 k eV up to 30 k eV at 100K. To obtain statistical ly meaningful results, for each value of PKA energy, cascade simulat ions using several crystallographic directions and positions of PKA have been carried out The results are then compared with Norgett RobinsonTorrens ( NRT ) model and previous simulation studies using a Finnis Sinclair potential 245. The effects of PKA orientation and grain boundaries are also discussed. 6 .2 Structure of Titanium Titanium has several phases (Figure 6 1 ). Among these titanium is the stable phase at room temperature and atmosphere; ature bodycentered cubic phase; is a high pressure phase with hexagonal structure In addition phase phase is martensitic and strongly hyster etic, as indicated in Figure 6 1 Figure 6 1. Phase diagram of titanium 248. titanium has a hexagonal closepacked (HCP) structure, with lattice constants of 2.95 and 4.68 in a and c directions. T he three planes with the highest densities are basal plane PAGE 151 151 (0001), prismatic plane (10 10), and pyramidal plane (10 11) respectively. The closed packed directions are <11 20>. Along the [0001] direction, the structure forms alternating layer with ABAB sequence. There are two atoms in each unit cell, with Wyckoff positions of (1/3, 2/3, 1/4) and (2/3, 1/3, 3/4) respectively. The ratio of c/a is 1.588, which is slightly lower than the ideal ratio 1.633, resulting in six neighbors at 2.90 and six at 2.95 titanium transform s to body centered cubic (BCC) titanium at high temperature which has a lattice constant of 3.32 at 900 oC 248. The most densely packed pla nes and directions f or BCC are (110) and <111> respectively Figure 6 248 titanium has a hexagonal structure, with space group P6/mmm 248. Each unit cell contains three atoms. the c axis. The first layer consist s of atoms in the (0, 0, 0) Wyckoff position. The second layer is formed from atoms at (2/3, 1/3, 1/2) and (1/3, 2/3, 1/2). The atoms in the first layers have a coordination number of PAGE 152 152 fourteen whereas the atoms in the next layer have eleven n eighbors. The density of the higher titanium, as it is a high pressure phase. As the The collapse of bcc planes in the [111] direction can generate the str ucture. Figure 6 3 248. 6 .3 Computational Details The cascade simulations are carried out using an in house MD code and two MEAM potentials, one of which is especially accurate in describing the transformation s among the phases of titanium. The details of the MEAM potential are given in section 6. 3.1. Since the PAGE 153 153 interatomic distances in cascade simulations are sometimes very short, the MEAM potentials are used in combination with the ZBL potential, whose details are prese nted in section 6.3.2. Periodic boundary conditions are applied in three dimensions with a Parrinello Rahman 91 barostat The Berendsen thermostat which is more reasonable in both physics and statistical mechanics than velocity rescaling, is employed to maintain the temperature of the syste m. The motions of the atoms are integrated via fifth order predictor corrector algorithm 64. The time step is chosen to be 0.2 fs to balance accuracy with the computational efficiency A t ime step less than that does not have a significant effect on the number of surviving Fre nkel defects (Figure 6 4 ). Prior to the cascade simulations, the systems are equilibrated at 100K for 2 0 ps. Typical cascade simulations times are 2 0 ps, which is considered to be long enough for the system to research a stable state, based on the analysis of the time evolution of the defects. Figure 6 4. The evolution of number of Frenkel defect with time in cascade simulations at 100K with PKA = 1 keV using different time steps. PAGE 154 154 The number of point defects (vacancies and interstitials) and defect clu sters are analyzed using two different approaches: common neighbor analysis and lattice matching analysis respectively (section 6.4). To achieve reasonable statistics, several cascade simulations are simulated for each case. The number of simulation s performed for each PKA energy and corresponding system size are given in Table 6 1. The equivalent neutron energies to PKA energies are also listed, assuming elastic collisions between two species. Table 6 1 The system size and number of simulations pre formed for each PKA energy. The neutron energy that is equivalent to average PKA energy is also listed. Neutron Energy (MeV) Average PKA Energy (KeV) Atoms in Simulations No. of Simulations Performed 0.0025 0.1 64000 6 0.0125 0.5 64000 6 0.025 1 64000 12 0.075 3 64000 12 0.125 5 216000 12 0.25 10 216000 18 0.75 30 216000 3 6 .3.1 MEAM Potential In the present study, two types of MEAM potentials have been used. One is the traditional MEAM, named MEAM_I hereafter. The total energy of the system is given as 246: i i j ij ij i i iR F E) (2 1 [6 1 ] w here Fi is the embedding function, i is the background electron density at the site i ij is the pair function between atom i and atom j Rij is the distance between atom i and j The embedding fu nction is given as 82: 0 0ln i i c iAE F [6 2 ] PAGE 155 155 w here A is an adjustable parameter, Ec is the cohesive energy, 0 is the background electron density for a reference structure, which is usually the equilibrated structure of titanium i is the background electron density at given site, which consists of several partial electron densities for different angular contributions, corresponding to s, p, d, and f orbital respectively. The partial electron density term are expressed as 82: 2 ) 0 ( 2 0 i j ij a jRi [6 3 ] 2 ) 1 ( 2 1 i j ij a j ij ijR R Ri [6 4 ] 2 ) 2 ( 2 ) 2 ( 2 2 23 1 i j ij a j i j ij a j ij ij ijR R R R Ri [6 5 ] 2 ) 3 ( 2 ) 3 ( 3 2 35 3 i j ij a j ij ij i j ij a j ij ij ij ijR R R R R R R Ri [6 6 ] w here j is the atomic electron density from atom j at distance Rij, R ij is the component ( =x, y, z) of the distance vector between atom i and j There are several ways of combining the partial densities to give the total background electron densities 249. The following formula is widely used 246: ) (0 Gi i [6 7 ] e G 1 2 ) ( [6 8 ] 3 1 2 0 h i h i h it [6 9 ] w here ti h are adjustable parameters. The atomic electron density from atom j is calculated by PAGE 156 156 ) 1 / ( ) () () ( e ij hr R ij h a je R [6 10] w here (h) are the decay le ngths and adjustable parameters; re is the near neighbor distance in the equilibrium reference structure. Although the embedding function is given in the specific form above, the pair interaction function given in Equation 6 1 is not uniquely defined. Instead, the total energy per atom of the equilibrated reference structure is given from a universal functional proposed by Rose et al. 250. The functional is expressed as *a c ue da a E R E ) 1 (3 [6 11] w here d is an adjustable parameter and ) 1 / (* er R a [6 1 2 ] 2 19 cE B [6 13] where B is the bulk modulus and is the equilibrium volume Once the total energy per atom and embedding function are computed for the reference structure, the pair interaction as a function of the nearest neighbor distance can be computed using Equation 6 1 It is noted that the formulas above only consider nearest neighbor interactions. The second and more distant neighbor interaction are described using a many bodyscreening function 251 253. For given atom i j and k a unique plane can be defined. An ellipse passing through atom i and j can be also defined using the following equation: 2 2 2) ( 2 1 1ijR y C x [6 14] PAGE 157 157 The central idea of this screeni ng function is to define two limiting values, which are Cmax and Cmin respectively. If the atom k is outside the ellipse defined by Cmax, it is assumed that atom k does not have any effect on the interactions between atom i and j If atom k is inside the ellipse define by Cmin, it is assumed that the atom k completely screens the interaction between atom i and j Between Cmax and Cmin, the screening interaction gradually changes. The screening function is defined as ikj j i k ijS S, [6 15] min max minC C C C f Sc ikj [6 16] w here 0 0 ) ( 1 0 ) 1 ( 1 )( 1 1 ) (2 4 x x f x xx f x x fc c c [6 17] Another type of MEAM potential is in the function form of cubic spline 83, thereafter named MEAM_II, which removes the constraints of fixed angular feature and allows additional flexibility of the potential. The total energy in MEAM_II is c alculated by 84 ij i i ijU r E ) ( ) ( [6 18] w here ij) is the pair interaction between atom i and j is the self energy function depending on the electron density of atom i The electron density of atom i is calculated using j jk jik ik ij ij ig r f r f r )] [cos( ) ( ) ( ) ( [6 19] PAGE 158 158 w here jik is the angle between atom j i and k centered on atom i All five functions ij) U ij) f(rij), and g [cos (jik)] are represented by cubic spline functions, which are fitted to the fundamental properties of the system. MEAM_II can be related to the Stillinger Weber 254 potential [ U(x)=x and =0 ] and EAM potential [ f=0 or g=0 ]. The parameters of MEAM_II are given in Table 6 2. Table 6 2 The parameters of spline function for MEAM_II 83. t t_min t_max N phi r [A] 1.742693 5.5 13 rho r [A] 2.055802 4.41 11 f r [A] 2.055802 4.41 10 U rho tot 55.1423 23.9383 4 g cos(theta) 1 0.928437 8 i phi(r_i) rho(r_i) f(r_i) U(rho_i) g(x_i) 0 3.7443 1.7475 0.1485 0.2975 0.0765 1 0.9108 5.8678 1.6845 0.1545 0.1416 2 0.3880 8.3376 2.0113 0.0510 0.7579 3 0.0188 5.8399 1.1444 0.5734 0.6301 4 0.2481 3.1141 0.2862 0.0905 5 0.2645 1.7257 0.3459 0.3574 6 0.2272 0.4429 0.6257 0.6529 7 0.1293 0.1467 0.6120 6.0091 8 0.0597 0.2096 0.3112 9 0.0311 0.1442 0.0000 10 0.0138 0.0000 11 0.0032 12 0.0000 i phi'(r) rho'(r) f'(r) U'(rho) g'(x) 0 20 1 2.773275 0.007769 8.336423 N 0 0 0 0.105198 60.4025 Several properties are calculated using these two MEAM potential and compared with other types of approaches, which includes EAM 255 and FS potential 256, DFT 83, and experimental results. The calculated properties include lattice constants, elastic constants and modulus, defect formation and migration energies, surface energies, stacking fault energy, thermal expansion coefficient, specific heat, melting point and enthalpy of melting. The comparison s are given in Table 6 3 and Table 6 4 In general, both MEAM potentials have a PAGE 159 159 titanium. However, further tests show that MEAM_I predicts a very low titanium, which results in incorrect structures during the MD tensile test. On the other hand, MEAM_II provides a reasonable description of the tra nsformation among different phases. Table 6 3 The fundamental properties calculated using MEAM_I and MEAM_II, compared with EAM, FS, DFT and experimental data. Ti EAM 255 F S 256 MEAM_I 246 MEAM_II 83 DFT 83 Exp. E c 4.831 257 a 2.922 2.945 2.931 2.948 2.951 258 c 4.772 4.687 4.678 5.171 4.679 258 c/a 1 0.975 0.974 1.596 1.583 0.971 258 B 1.099 1.097 1.097 259 C 11 1.954 1.8 1.701 1.74 1.72 1.761 259 C 12 0.737 0.873 0.804 0.95 0.82 0.869 259 C 44 0.481 0.514 0.421 0.58 0.45 0.508 259 C 33 2.067 2.17 1.871 1.88 1.90 1.905 259 C 13 0.609 0.766 0.748 0.72 0.75 0.683 259 E hcp 0.024 0.07 260 E hcp 0.048 0.06 260 Table 6 4 The fundamen tal properties calculated using MEAM_I and MEAM_II, compared with EAM, FS, DFT and experimental data. Ti EAM 255 F S 256 MEAM_I 246 MEAM_II 83 DFT 83 Exp. E v 1.48 1.43 1.79 2.24 2.03 >1.50 261 E m 0.82 1.09 0.96 0.87 E dv 2.73 3.87 4.00 2.58 C 3.82 3.070.2 3.9 2.21 2.87 O 3.45 3.070.2 4.53 2.64 2.58 T 3.39 3.070.2 unstable unstable unstable B asal (0001) 993 2144 1474 1939 2100 261 Prism(1 1 00) 1061 2145 1554 2451 1920 261 P rism(11 2 0) 1187 2352 1682 1875 6 .3.2 ZBL Potential The interatomic distances in the cascade simulations may be very short compared with traditional MD, due to the extremely high kinetic energy of the atoms. Therefore, to reproduce a realistic potential surface energy under such extreme conditions, a short rang e interaction has PAGE 160 160 been developed by Ziegle r Biersa c k and Littmark (ZBL) 247, in which the interatomic interactions are described by a universal form over a large number of the materials. The ZBL potential is based on the Thomas Ferm i and Linhard model and can be considered as a screened electrostatic potential of the nucleon nucleon interaction 262. The form of the ZBL potential is: ijij j ir r e Z Z r V 2 [6 20] w here Zi and Zj are the atomic number of the atoms, e is the charge of electron, rij is the interatomic distance, and (rij) is the screen function. The screen ing function can be expressed as 4 1) / exp(i u ij i i ija r b A r [6 21] Bohr j i uaZ Z a23 0 23 08854 0 [6 22] w here Ai and bi are predefined paramet ers, aBohr is the Bohr radius (0.53917 7 ). Table 6 5 The predefined parameter s of the ZBL potential 247. i A i b i 1 0.1818 3.2 000 2 0.5099 0.9423 3 0.2802 0.4028 4 0.02817 0.2016 6 .4 Defect Analysis Method In the current studies, both common neighbor analysis (CNA) and lattice matching analysis (LMA) have been employed. In essence, these two methods characterize similar aspects of the bonding and local environments of the vacancies and inters titials. The numbers of surviving stable defects using both methods are found to be consistent with each other. However, they differ in several aspects, which are discussed in details in the following PAGE 161 161 6 .4.1 Common Neighbor Analysis The CN A 263 method identifies a given atoms nearest neighbors within a certain cutoff radius. titanium, a cutoff value of 3.37 has been used to identify the first nearest neighbors. This is much longer than the bond length even with consideration of the thermal vibrations but shorter t han the second nearest neighbor distance In general, CNA examines the common neighbors of specific atom and provide s the following information : (1) the total number of the atoms in the system, (2) the total number of the bonds of each atom, (3) the number of the atoms in each subset with a different coordination number. Furthermore, CNA has been extended to distinguish the stacking sequence of HCP and FCC 263. In both structures, the numbers of first nearest neighbors and sec ond nearest neighbors are the same. Therefore, third nearest neighbor information is needed to achieve this. It is noted this extended analysis has been solely applied to the atoms with twelve neighbor s. If a n atom with tw elve coordination number is neithe r HPC nor FCC, it is indicated as O thers A typical analysis results is given in Table 6 6 Unfortunately, CNA cannot directly provide either the accurate positions for or the number of vacancies and interstitials. For instance, if a titanium vacancy is created in the system, CNA characterizes a certain number of the atoms having a less neighbor than they should, providing the local environment of vacancy position, from which the position of vacancy may be indirectly extracted. However, at finite temperature, the thermal vibration and diffusion process will make the detection of the accurate posit ion of a vacancy very difficult. Furthermore, if many vacancies are present, there is no systematic way to determine how many vacancies are in the system, as the number depends on the spatial distribution of the vacancies and interstitials. Overall, the nu mber of the defects can be indirectly estimated using CNA, the accurate number and positions are yet very difficult to obtain. PAGE 162 162 Table 6 6 Coordination number distributions provided by CNA with distinguishing HPC and FCC stacking sequence. *Others means o ther types of stacking sequence, not FCC and HCP A) Value of cut off: 3.370 B) # of Coordination # of atoms 10 : 12 11 : 110 12 : 63835 13 : 43 14 : 0 15 0 Total : 64000 C) Types of structure # of atoms among 12 coord ination HCP : 63779 FCC : 2 Others : 54 Total : 63835 6 .4.2 Lattice Matching Analysis LMA is an approach that compares the positions of atoms with reference lattice sites. The reference lattice is usually chosen to be the structure after long time equilibration at certain temperature. The schematic of LMA is illustrated in Figure 6 5 For a given lattice site of reference structure, LMA searches atoms within a certain radius of that site typically 40% of the bond length If there is an atom within the radius, the site is considered as be ing occupied Otherwise, the site is regarded as a v acancy. If an atom does not belong to any site, then it is considered as an interstitial. Therefore, LMA provides the positions of vacancies and interstitial. Moreover, the number of defect s is easily calculated within LMA. Nevertheless, the quality of LMA depends on the reference lattic e, unlike CNA, which solely relies on the positions of atoms. In some cases, the reference lattice is difficult to define, such as polycrystalline system s with grain boundaries, in which CNA is more straightforward. The refore, LMA and CNA are complement ary to each other. Due to this, both methods are employed in this study. PAGE 163 163 Figure 6 5 The schematics of lattice matching analysis (LMA). Both interstitials and vacancies are shown. 6 5 Single Crystal Results and Discus sions The threshold energy, Ed, is probably the most important single parameter that can be used to provide estimates about the production of survival Frenkel defects in cascade simulations 245. If the PKA energy is less than the threshold energy, atoms undergo large amplitude vibrations, without forming stable Frenkel pairs. To obtain its value, cascade simulations are carried out along various crystallographic directions. The orientation dependence of Ed titanium using MEAM_I potential is shown i n Figure 6 6 compared with a previou s study using the FS potential. PAGE 164 164 Figure 6 6 .Threshold energies at various crystallographic orientations using MEAM_I, compared with previous results using FS potential. It is found that there is a significant variation in the threshold energy, ranging from a minimum value of 14 eV in [0122] direction to the maximum number of 52 eV in [ 1211] direction using MEAM_I In addition, several orientatio ns have threshold energies less 20 eV. It is found that there is a head on collision between atoms when PKA follows these low energy directions. In the case of high index PKA directions, such collision does not happen. Since the displacements in the cascade simulations are easily affected by these directions, this can be used to estimate the upper limit of the surviving defects. d dam NRTE E N 2 / 8 0 [6 2 3 ] where Edam is the damage energy and Ed is the average threshold energy. It is noted that the damage energy is only a portion of the PKA, which is lost during the elastic collision process. The other part energy lost is due to the electronic excitation and ionization. However, MD PAGE 165 165 simulations do not consider such electronic effects. Therefore, in MD simulations, the damage energy is approximated equivalent to the PKA energy. Overall, MEAM_I and FS potential s predict similar trend, although certain differences exist in different orientations. The average value of the threshold energy is estimated to be approximately 25 eV with MEAM_I which is smaller than value predicted using the FS potential (30 eV). Therefore, from the NRT relationship (Equation 6 23), the number of stable survival Frenkel defects predicted by MEAM_I is estimated to be higher than that of FS potential. Once the threshold energy is obtained, the cascade simulations are carried out in single crystal titanium. During such simulations, a PKA with high kinetic energy is introduced into the system after long time equilibration, creating a large number of defects within a short period of time Initially, high index direction s of PKA, such as [13 45], [1124], and [1121], are chosen to avoid channeling effects. However, to explore anisotropic effects, some low index direction simulations are also carried out. For a particular PKA energy at a specific orientation, the position of the PKA is also varied to achieve average effects. For all cascade simulations in a single cr ystal system, the defect production and evolution process are very similar. A detailed description of cascade simulation when EPKA = 1 keV at 100K is given below. The snapshot of cascade simulation with PKA = 1keV at 100 K is shown in Figure 6 7. It is ob served that initially the number of defects increases very rapidly until a maximum number of Frenkel defect reached. This is the ballistic phase of the cascades usually last 0.3 ~ 0.4 ps. A highly disordered region is formed around the PKA within such phas e. After ~ 0.4 ps, recombination of vacancies and interstitial begins, which results in a decrease in the number of defects. This subsequent process is considered as the relaxation phase. Soon after, the number of surviving Frenkel defects PAGE 166 166 reaches a stead y state, accompanied by lattice restoration, leaving a small amount of defects and defect clusters inside the system. Figure 6 7 Snapshot of cascade simulation in single crystal system at different simulation time. Only the atoms with coordination rather than 12 are shown. To achieve statistically meaningful results, several cascades simulations with different PKA directions and positions have been carried out for a given value of PKA energy. The numbers of Frenkel defects as a function of simulation ti me for PKA energy of 1 keV, 3 keV, 5 keV and 10 keV of PKA are shown in Figure 6 8. It can be seen that there is a peak for each curve, which defines the transition from the ballistic phase to relaxation phase. The height of the peak, which is the maximum number of Frenkel defects created during the simulations, increases PAGE 167 167 with increases in the PKA energy. The time to reach a steady state also increases as a function of PKA energy. For instance, for PKA energy with 1keV, the numbers of surviving Frenkel def ects reach a steady state after 4 ps. For 10 keV simulations, it requires approximately 6 ps to reach steady state. From Figure 6 8, it can be also seen that all simulations with the same PKA energy nearly overlap with each other, indicating reasonable sta tistics has been achieved. The average value for a given PKA energy has been taken to establish the relationship between PKA energy and surviving Frenkel defects. Figure 6 8 .The number of survival Frenkel defects as a function of simulation time with va rious directions, positions, and energies of primary knockoff atom. In principal high index direction s for the PKA should be used to avoid channeling effects. To investigate the effects of orientations and lattice symmetry, simulations with low index d irections have also been carried out. For each value of PKA energy, PKA with different positions a long [13 4 5], [11 2 4], [00 0 1], [11 2 1] and [112 0] are simulated. The average PAGE 168 168 numbers of survival Frenkel defects with various PKA energies are shown in Figure 6 9 No significant variation is observed at any energy Therefore, the effects of anisotropy of the HCP structure are rather weak, which indicate s that the cascade feature is insensitive to anisotropy of the crystal structure. Figure 6 9 .The e ffect of PKA directions on survival Frenkel defects. The dependence of the number of surviving Frenkel defects on the P KA energy is shown in Figure 6 10. It is seen that the number of Frenkel pairs using MEAM_I and MEAM_II is higher than that of FE potenti al at given PKA energy. This is expected, as the threshold energy is lower for FS potential. Since the number of Frenkel defects and the number predicted by NRT model have a defin ite relationship, NNRT prediction is also shown in Figure 610. It is found that the relative ratio of NF with respect to NNRT using MEAM potentials is higher than FS potential PAGE 169 169 Figure 6 10.The comparison of number of survival Frenkel defects as a function of PKA energies between MEAM_I, MEAM_II, and FS potentials. The numbers predicted by NRT model with difference threshold energies are also indicated. Detailed analysis leads to the conclusion that t he large difference in surviving Frenkel defects is considered to be caused by the potentials. Although these three potential predict similar value s of cohesive energy, lattice parameters, c/a ratio, elastic constants, vacancy and interstitial energy, and surface energies for equilibrium structure of Ti, they are different in several aspects. First, the potentials were fitted t o different structure s and data. The FS potential w as fitted to an imaginary FCC lattice which does not exist in the phase diagram of titanium, with the nearest neighbor distance equal to the lattice parameter, a, of the hexagonal lattice and with a negat ive stacking fault energy, assuming the thirdnearest neighbor interaction is just a small perturbation compared with first and second nearest neighbor interactions. This can be done PAGE 170 170 because the first and second nearest neighbor environment of FCC and id eal HCP structure are identical; the negative stacking fault energy assures the preference for the HCP structure. By contrast the MEAM_I parameters were determined by fitting to the fundamental properties of HCP structure of titanium, which includes ela stic constants, surface energy, stacking fault energy, and vacancy formation energy. In addition, the stability of the HCP structure at finite temperatures was verified by considering a thermal expansion coefficient for MEAM_ I. For MEAM_ II, the parameters were fitted to the physical properties of several phase of titanium, FCC and simple hexagonal phases, to improve the accuracy and extend the applicability. In addition, a DFT datab ase for test of the potential was employed during the fitting process of MEAM_II. It included As a result the phase transformation is reasonably descr ibed by MEAM_II. Second, the short range interactions of these potentials are very different. The potential energy (Eu) as a function of nearest neighbor distance is given in Figure 6 1 1 It is noted that the Eu (in Equation 6 11) is function for uniform expansion or contraction in the reference structure. From Figure 6 1 1 it can be seen that MEAM_ I is least repulsive at short range. Consistent with this, MEAM_ I predicts the highest number of stable Frenkel pair after cascade simulations. MEAM_ II is more repulsive than MEAM_ I, especially in the range of 1 1.5 However, it is less repulsive than FS when the interatomic distance is shorter than ~ 0.7 or larger than ~1.6 The FS potential is the most repulsive one when the first nearest neighbor dist ance is larger than ~1.6 In agreement with this, FS predicts the least number of stable Frenkel defects in the cascade simulations. PAGE 171 171 Figure 6 1 1 The potential energy as a functional of nearest neighbor distance using FS, MEAM_2006 and MEAM_2008 potentials. In addition, different cutoff s have been used when connecting these potentials with ZBL potential. For the FS potential, the ZBL potential is used for spacing less than 1.1 and FS potential is used when spacing greater than 2.09 An exp onential potential with four parameters is employed to join these two potential. For MEAM_2008, the potential is used when interatomic distance is larger than 1.7427 and ZBL is used when distance is shorter than 1.5 From Figure 6 11, it is also obse rved that MEAM_2008 and FS are very similar when interatomic distance is longer than 2.2 Therefore, if this is used as a cutoff for both MEMA_2008 and FS, they are expected to predict similar results. In summary, t he differences in the cutoff and connec ting functions lead to significant changes in the overall potential energy functions, which have a significant effect on the radiation damage and defect distributions PAGE 172 172 6.6 Effects of Grain Boundaries Grain boundaries can be considered as sources and sinks of defects. Therefore, they can have a significant influence on the defect generation, evolution, and distribution which will further determine the microstructure s and mechanical properties of materials. There are experimental indicati ons that metals with grain size s in the nano scale may be more radiation resistant than that with coarse grains 264. However, the underlying reasons are not fully understood. Therefore, it is of great importance to perform cascade simulations with the presence of grain boundaries. Polycrystalline titanium is generated with four grain s for cascade simulations, with a 2D hexagonal texture. Each hexagonal is filled with a single crystal titanium orientated with [0001] along the columnar direction. The in plane orientations are chosen in such a way that t he grain boundaries between grains are high energy tilt grain boundaries, ensuring the microstructure is stable. The grain size of polycrystalline system is approximately 20 nm, with 252838 atoms in the system. Figure 6 1 2 Polycrystalline structure of titanium in the cascade simulations. Courtesy of D. H. Kim. PAGE 173 173 Similar to the single crystal study, various direction and position of PKA have been simulated for a given value PKA energy. It is found that the effects of grain boundaries can be classified in to two cases. First, when the PKA energy is low (~ 1 keV), if the PKA is located at the center region of one grain, the effects of grain boundaries are insignificant. This case is very similar to that in the single crystal, as the amorphous region created by the PKA does not interact with the grain boundaries. Consistent with this, the number of surviving Frenkel defects is similar to that in the single crystal. Second, there are significant effects of the grain boundaries. This could happen when the PKA energy is larger than 3 keV or the position of PKA is close to the grain boundaries. For example, when PKA = 10 keV, the defect evolution process is shown in Figure 6 1 3 It can be seen that the defects created by PKA quickly reach the grain boundaries an d interact with them. Furthermore, it is found that grain boundaries also have a significant effect on the relaxation p hase, limiting the ability of atoms to return to their equilibrium positions. Therefore, the number of surviving Frenkel defects in polyc rystalline system for such case is much higher than that in the single crystal. In addition, a portion of the defects created by PKA become part of the grain boundaries leading to a grain boundary growth. The cascade morphologies are found to be strongl y influenced by the structure of grain boundaries. In summary cascades simulations have been carried in polycrystalline titanium system s including various PKA energies, directions, and positions. It is found that when the PKA energy is small and PKA i s far from the grain boundaries, there is no strong interaction between them, generating similar number of surviving Frenkel defects as in the single crystal. However, when PKA energy is relatively large or PKA is close to the grain boundaries, numerous points defects are created along the interfaces, leading to a grain boundary growth. The numbers of surviving PAGE 174 1 74 Frenkel defects in such cases are much higher than that in single crystal cascade simulations. It is also found that grain boundary structure s in th e polycrystalline system s after cascade simulations are highly disordered and the morphologies strongly depend on these grain boundary structures. Figure 6 13. Snapshot of cascade simulation of polycrystalline system when PKA = 10 keV at 100K. Only the atoms with coordination rather than 12 are shown. PAGE 175 175 CHAPTER 7 SUMMARY AND FUTURE W ORK In this dissertation, the structure and stability of defects in structurally complex materials have been determined and predicted using atomic level simulations. The production, evolution, and dynamic behavior of defects have been explored A great effort has been focused on establish ing the relationship between fundamental properties of bulk materials and effects of point defects. 7.1 Lithium Niobate A large experimental body of literature on lithium niobate, a technologically important ferroelectric, suggests that non stoichiometric defects dominate its physical behavior, from macroscale switching to nanoscale wall structure. The exact structure and energetics of such proposed intrinsic def ects and defect clusters remained unverified by either first principles calculations or experiments. Here, densityfunctional theory (DFT) was used to determine the dominant intrinsic defects in LiNbO3 under various conditions. In particular, in an Nb2O5rich environment, a cluster consisting of a niobium antisite compensated by four lithium vacancies was predicted to be the most stable defect structure, thereby verifying a conjecture in the lite rature Under Li2O rich conditions the lithium Frenkel defect was predicted to be the most stable with positive defect formation energy (DFE). This was proposed as the underlying reason that the vapor transport equilibration (VTE) method can grow stoichi ometric LiNbO3. The effects of temperature and oxygen partial pressure were also explored by combining the DFT results with thermodynamic calculations The relative stabilities of various defect cluster arrangements of lithium vacancies around a niobium antisite in LiNbO3 were determined. Their effects on the ferroelectricity of the system were also discussed. It was found that at room temperature the non uniaxial dipole moments associated with the defect clusters could locally PAGE 176 176 affect the properties of the system. The diffusion mechanism was predicted to be through first nearest neighbor jumps on the Li sublattice. The diffusivity of lithium vacancy was found to be extremely low at r oom temperature, which indicated that the defect complexes should be rather stable. These predictions provide d a picture of a very rich defect structure in lithium niobate, which has important effects on its physical behavior at the macroscale. The incorporation site preference and corr esponding charge compensation mechanisms of several impurities in LiNbO3 were predicted based on DFT calculations of possible defect clusters at congruent and stoichiometric growth conditions respectively. In general, it was found that impurities occupy lithium sites compensated by lithium vacancies under Nb2O5 reference state. In addition impurities on both lithium can niobium site were predicted to be dominant under Li2O reference state. However, due to the effects of intrinsic defects and method of in corporation impurities, the concentration of impurities on niobium site was predicted to be negligible for stoichiometric sample. Interestingly, it was found that Mg had a strong tendency to occupy both lithium and niobium sites even under Nb2O5 rich condi tions. Increasing the concentration of Mg would lead to an increase in the value of its chemical potential, which would further result in a switch of dominant defects from Mg on lithium site compensated by intrinsic lithium vacancies to Mg on both sites. The prediction is consistent with the experimental observation that Mg is mainly on lithium site s at low dopant concentration and on both sites at high concentrations. It was also found that Fe3+ is more energetically favorable than Fe2+ under both referen ce states. Therefore, without special treatment, the Fe3+ was predicted to be the majority species. In agreement with this, the experimental concentration of Fe2+ is usually at least one order of magnitude lower than the Fe3+ concentration The thermal ionization energies and charge transfer levels of impurities were calculated, especially for transition metal PAGE 177 177 iron. The charge transfer level of Fe2+/Fe3+ was calculated t o be 1.34 eV, which corroborated the lo ng time conjecture that there wa s a charge transfer effect in the Fe doped LiNbO3. The relative position of this charge transfer level with regard to the band structure was illustrated, with consideration of the effects of intrinsic defects NbLi 4 +/NbLi 5 +. The charge transfer level s of Mg Nd, and Er and were also calculated and discussed. The suggested future work is to investigate the effects of defect s at interface s such as polarization domain wall s It is known that domain wall s are atomic sharp from simulations, usually a few unit c ells wide However, the experimental measurements of such width are of the order or dozens nanometers. Studies of point defects and defect clusters would contribute to resolving this discrepancy. These point defects could introduce both short range electr ostatic interaction and long range elastic interaction into the system, which could locally change the domain wall structure and may cause a broadening effects. Therefore, it is of great significance to simulate the defects at domain wall to determine ho w they may affect the domain wall width. Another direction would be to explore the effects of external field s in particular electrical field In experiments, the cohesive field of congruent materials is almost a n order of magnitude larger than that of stoichiometric lithium niobate. The hypothesis is that such difference s are caused by point defects. DFT calculations can provide theoretical evidence to corroborate or refute such a conjecture. It has been shown tha t the defect clusters have dipole moments For switching the polarization, first nearest neighbor jumping is requisite for defect clusters. Due to the high energy barrier, this jumping process is rather difficult at room temperature. Nevertheless o nly sh ort distance movements are needed f or polarization inversion of bulk atoms which is relatively easy. Therefore, by simulating defects under electrical field, the PAGE 178 178 polarization reversal process can be elucidated. The corresponding structural changes of defect clusters and their effects on the other fundamental properties on the materials can be determined. 7.2 Ceria based Materials We critically assess ed the materials fidelity of six interatomic potentials for ceria, based on predicted lattice constants, t hermal expansion, chemical expansion, dielectric properties, oxygen migration energy and mechanical properties. While, no potential can reproduce all fundamental properties, the Gotte (2007) and Grimes potentials display the combination of highest fidelity widest range of applicability. The simulations show ed that sub stoichiometry lead to a significant softening of the elastic constant s which was consistent with experimental results Similar results were observed for doped ceria systems. The simulations show ed that both sub stoichiometry and aliovalent doping of ceria lead to large decreases in the elastic constants of the material using the best available potential(s) for each physical property These decreases arise from the significantly reduced stren gth of ionic interactions. If this elastic softening were to be so large as to threaten the mechanical integrity of the system, then it would pose a challenging problem for the design of ceria based electrolytes. It would thus be of significant value to i dentify a strategy for increasing the ionic conductivity that does not result in a significant softening in the elastic properties. A co doping strategy is a possible solution to achieve both high ionic conductivity and mechanical strength. U sing multiple trivalent dopants could suppress the ordering effects of oxygen vacancy and increase in the configurational entropy of the system, which result s in enhance ment of the ionic conductivity. In addition, due to the mismatch between host cations and dopant ions strain s can be created, which may strengthen the materials. PAGE 179 179 7.3 Titanium Cascade simulations were titanium, including various energies, positions and orientations of primary knock on atoms to obtain statistical ly meaningful results. The dependence of Frenkel defects on the time was established and two different phase s during the cascade simulation were observed : the ballistic and relaxation phase. The effects of PKA energies on survival Frenkel defects were evaluated and compared with previous study using FS potentials. It was found that both MEAM potential in this study predict ed higher value of the survival Frenkel defects. This is because both MEAM potentials predict lower average threshold energ ies whic h result s in a higher number of defects in the NRT model. The orientation effects of the PKA were also evaluated and no significant effect were found. The effects of grain boundaries were evaluated For the practical application of titanium as a cladding m aterials, the fundamental understanding of the mechanical properties, especially the plastic deformation behavior of the initial stage during the radiation, is required. There are some experimental indications that nano crystalline materials have a higher tolerance of radiation damage. Therefore it is crucial to perform mechanical test on the polycrystalline system after cascade simulations. Some simulations might p redict microstructure evolution with atomic resolution. Fu rthermore, these simulations might also provide useful insights about the deformation mechanisms such as dislocation and twinning. PAGE 180 180 LIST OF REFERENCES 1 M. W. Barsoum, Fundamentals of Ceramics (Taylor & Francis, Bodmin,Cornwall,UK, 2002). 2 W. D. Callister, Materials Science and Engineering: An Introduction (John Wiley & Sons Inc., 2002). 3 J. H. Spaeth and H. Overhof, Point Defects in Semiconductors and Insulators (Springer, New York, USA, 2003). 4 J. Garcia Sole, L. E. Bausa, and D. Jaque, An Introduction to the Optical Spectroscopy of Inoganic Solids (John Wiley & Sons Ltd., West Sussex, England, 2005). 5 G. H. Dieke, Spectra and Energy Levels of Rare Earth Ions in Crystals (Interscience Publishers, New York, USA, 1968). 6 Y. Tanabe and S. Sugano, J. Phys. Soc. Jpn. 9 753 (1954). 7 Y. Tanabe and S. Sugano, J. Phys. Soc. Jpn. 9 766 (1954). 8 H. Holleck, Journal of Vacuum Science & Technology a Vacuum Surfaces and Films 4 2661 (1986). 9 X. Jiang, M. Wang, K. Schmidt, E. Dunlop, J. Haupt, and W. Gissler, J. Appl. Phys. 69, 3053 (1991). 10 V. Gopalan, V. Dierolf, and D. A. Scrymgeour, Annu. Rev. Mater. Res. 37, 449 (2007). 11 K. Duncan, in Department of Materials Science and Engineering (University of Florida, Gainesville, 2001), Vol. Ph.D. 12 K. K. Wong, Pro perties of Lithium Niobate (INSPEC, Stevenage, Herts, UK, 2002). 13 S. Erdei and V. T. Gabrieljan, Cryst. Res. Technol. 24, 987 (1989). 14 H. L. Wang, Y. Hang, J. Xu, L. H. Zhang, S. N. Zhu, and Y. Y. Zhu, Mater. Lett. 58, 3119 (2004). 15 S. C. Abrahams an d P. Marsh, Acta Crystallogr., Sect. B: Struct. Sci. 42, 61 (1986). 16 K. Kitamura, J. K. Yamamoto, N. Iyi, S. Kimura, and T. Hayashi, J. Cryst. Growth 116, 327 (1992). 17 K. Kitamura, Y. Furukawa, K. Niwa, V. Gopalan, and T. E. Mitchell, Appl. Phys. Lett. 73, 3073 (1998). 18 F. S. Chen, Journal of Applied Physics 40, 3389 (1969). PAGE 181 181 19 R. J. Holmes and D. M. Smyth, J. Appl. Phys. 55, 3531 (1984). 20 R. L. Holman, Processing of Crystalline Ceramics (Plenum, New Yourk, 1978). 21 V. Gopalan and M. C. Gupta, Appl. Phys. Lett. 68, 888 (1996). 22 V. Gopalan, T. E. Mitchell, Y. Furukawa, and K. Kitamura, Appl. Phys. Lett. 72, 1981 (1998). 23 A. Prokhorov and I.Kuzminov, Physics and Chemistry of Crystalline Lithium Niobate (Hilger, Bris tol, New York, 1990). 24 O. F. Schirmer, O. Thiemann, and M. Wohlecke, J. Phys. Chem. Solids 52, 185 (1991). 25 R. Mouras, M. D. Fontana, P. Bourson, and A. V. Postnikov, J. Phys.: Condens. Matter 12, 5053 (2000). 26 D. A. Bryan, R. Gerson, and H. E. Tomas chke, Appl. Phys. Lett. 44, 847 (1984). 27 D. L. Staebler and W. Phillips, Appl. Opt. 13, 788 (1974). 28 E. Lallier, Appl. Opt. 31, 5276 (1992). 29 E. Lallier, J. P. Pocholle, M. Papuchon, C. Grezesbesset, E. Pelletier, M. Demicheli, M. J. Li, Q. He, and D B. Ostrowsky, Electron. Lett. 25, 1491 (1989). 30 L. F. Johnson and A. A. Ballman, J. Appl. Phys. 40, 297 (1969). 31 S. J. Field, D. C. Hanna, D. P. Shepherd, A. C. Tropper, P. J. Chandler, P. D. Townsend, and L. Zhang, Opt. Lett. 16, 481 (1991). 32 V. D ierolf, T. Morgus, C. Sandmann, E. Cantelar, F. Cusso, P. Capek, J. Spirkova, K. Polgar, W. Sohler, and A. Ostendorf, Radiat. Eff. Defects Solids 158, 263 (2003). 33 V. Dierolf, C. Sandmann, S. Kim, V. Gopalan, and K. Polgar, J. Appl. Phys. 93, 2295 (2003). 34 V. Dierolf, C. Sandmann, V. Gopalan, S. Kim, and K. Polgar, Radiat. Eff. Defects Solids 158, 247 (2003). 35 H. Sothe and J. M. Spaeth, J. Phys.: Condens. Matter 4 9901 (1992). 36 W. Keune, S. K. Date, I. Dezsi, and U. Gonser, J. Appl. Phys. 46, 3914 (1975). 37 W. Keune, S. K. Date, U. Gonser, and H. Bunzel, Ferroelectrics 13, 443 (1976). 38 L. Kovacs, L. Rebouta, J. C. Soares, M. F. Dasilva, M. Hageali, J. P. Stoquert, P. Siffert, C. Zaldo, Z. Szaller, and K. Polgar, Mater. Sci. Eng., B 9 505 ( 1991). PAGE 182 182 39 T. Gog, M. Griebenow, and G. Materlik, Phys. Lett. A 181, 417 (1993). 40 T. Gog, T. Harasimowicz, B. N. Dev, and G. Materlik, Europhys. Lett. 25, 253 (1994). 41 C. Prieto and C. Zaldo, Solid State Commun. 83, 819 (1992). 42 L. Rebouta, M. F. Dasi lva, J. C. Soares, M. Hageali, J. P. Stoquert, P. Siffert, J. A. Sanzgarcia, E. Dieguez, and F. Agullolopez, Europhys. Lett. 14, 557 (1991). 43 I. Baumann, R. Brinkmann, M. Dinand, W. Sohler, L. Beckers, C. Buchal, M. Fleuster, H. Holzbrecher, H. Paulus, K H. Muller, T. Gog, G. Materlik, O. Witte, H. Stolz, and W. vonderOsten, Appl. Phys. A: Mater. Sci. Process. 64, 33 (1997). 44 T. Nolte, T. Pawlik, and J. M. Spaeth, Solid State Commun. 104, 535 (1997). 45 D. M. Gill, J. C. Wright, and L. Mccaughan, Appl. Phys. Lett. 64, 2483 (1994). 46 D. M. Gill, L. McCaughan, and J. C. Wright, Phys. Rev. B 53, 2334 (1996). 47 O. Witte, H. Stolz, and W. vonderOsten, J. Phys. D 29, 561 (1996). 48 V. V. Kharton, F. M. Figuei redo, L. Navarro, E. N. Naumovich, A. V. Kovalevsky, A. A. Yaremchenko, A. P. Viskup, A. Carneiro, F. M. B. Marques, and J. R. Frade, J. Mater. Sci. 36, 1105 (2001). 49 M. Mogensen, N. M. Sammes, and G. A. Tompsett, Solid State Ionics 129, 63 (2000). 50 V. V. Kharton, F. M. B. Marques, and A. Atkinson, Solid State Ionics 174, 135 (2004). 51 T. H. Etsell and S. N. Flengas, Chemical Reviews 70, 339 (1970). 52 V. V. Kharton, E. N. Naumovich, and A. A. Vecher, Journal of Solid State Electrochemistry 3 61 (1999). 53 S. P. S. Badwal, Solid State Ionics 52, 23 (1992). 54 R. M. Dell and A. Hooper, Solid Electrolytes (Academic Press, New York, 1978). 55 H. Inaba, R. Sagawa, H. Hayshi, and K. Kawammura, Solid State Ionics 122, 95 (1998). 56 T. Inoue, T. Setoguchi, K. Eguchi, and H. Arai, Solid State Ionics 35, 285 (1989). 57 A. Trovarelli, Catalysis Reviews Science and Engineering 38, 439 (1996). 58 L. Minervini, M. O. Zacate, and R. W. Grimes, Solid State Ionics 116, 339 (1999). 59 E. A. Kummerle and G. Heger, Journa l of Solid State Chemistry 147, 485 (1999). PAGE 183 183 60 B. C. H. Steele, High Conductivity Solid Ionic Conductors (World Scientific, Singapore, 1989). 61 S. A. Fabritsiev and A. S. Pokrovsky, Fusion Eng. Des. 65, 545 (2003). 62 R. M. Martin, Electronic Structure: Basic Theory and Practical Methods (Cambridge University Press, New York, USA, 2004). 63 D. S. Sholl and J. A. Steckel, Density Functional Theory: A Practical Introduction (Wiley, Hoboken,New Jersey, 2009). 64 M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Oxford University Press, Oxford, 1987). 65 C. Kittel, Introduction to Solid State Physics (John Wiley & Sons Inc., Hoboken, NY, 2005). 66 M. Born and R. Oppenheimer, Annalen Der Physik 84, 0457 (1927). 67 P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964). 68 W. Kohn and L. J. Sham, Phys. Rev. 140, 1133 (1965). 69 J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992). 70 J. P. Perdew, K. Burke, and M. Ernzerhof, Chemical Applications of DensityFunctional Theory 629, 453 (1996). 71 G. B. Bachelet, D. R. Hamann, and M. Schluter, Phys. Rev. B 26, 4199 (1982). 72 M. Veithen and P. Ghosez, Phys. Rev. B 65 (2002). 73 K. Parlinski and Z. Q. Li, Phys. Rev. B 61, 272 (2000). 74 Q. K. Li, B. Wang, C. H. Woo, H. Wang, and R. Wang, J. Phys. Chem. Solids 68, 1336 (2007). 75 P. E. Blochl, Phys. Rev. B 50, 17953 (1994). 76 G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999). 77 N. A. W. Holzwarth, G. E. Matthews, R. B. Dunning, A. R. Tackett, and Y. Zeng, Phys. Rev B 55, 2005 (1997). 78 A. I. Liechtenstein, V. I. Anisimov, and J. Zaanen, Phys. Rev. B 52, R5467 (1995). 79 S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, and A. P. Sutton, Phys. Rev. B 57, 1505 (1998). PAGE 184 184 80 D. Wolf, P. Keblinski, S. R. Phil lpot, and J. Eggebrecht, J. Chem. Phys. 110, 8254 (1999). 81 S. M. Foiles, M. I. Baskes, and M. S. Daw, Phys. Rev. B 33, 7983 (1986). 82 M. I. Baskes, Phys. Rev. B 46, 2727 (1992). 83 R. G. Hennig, T. J. Lenosky, D. R. Trinkle, S. P. Rudin, and J. W. Wilki ns, Phys. Rev. B 78 (2008). 84 T. J. Lenosky, B. Sadigh, E. Alonso, V. V. Bulatov, T. D. de la Rubia, J. Kim, A. F. Voter, and J. D. Kress, Modelling and Simulation in Materials Science and Engineering 8 825 (2000). 85 M. P. Allen and D. J. Tildesley, Com puter Simulation of Liquids (Oxford Science, Oxford, 2004). 86 C. Gear, Numerical Initial Value Problems in Ordinary Differential Equation (Prentice Hall, Englewood Cliffs, NJ, 1971). 87 Y. H. Hu and S. B. Sinnott, J. Comput. Phys. 200, 251 (2004). 88 K. H uang, Statistical Mechanics (John Wiley & Sons Inc., New York, 1987). 89 H. J. C. Berendsen, J. P. M. Postma, W. F. Vangunsteren, A. Dinola, and J. R. Haak, J. Chem. Phys. 81, 3684 (1984). 90 H. C. Andersen, J. Chem. Phys. 72, 2384 (1980). 91 M. Parrinello and A. Rahman, J. Appl. Phys. 52, 7182 (1981). 92 G. D. Boyd, K. Nassau, R. C. Miller, W. L. Bond, and A. Savage, Appl. Phys. Lett. 5 234 (1964). 93 L. Tian, V. Gopalan, and L. Galambos, Appl. Phys. Lett. 85, 4445 (2004). 94 P. F. Bordui, R. G. Norwood, D. H. Jundt, and M. M. Fejer, J. Appl. Phys. 71, 875 (1992). 95 C. Baumer, C. David, A. Tunyagi, K. Betzler, H. Hesse, E. Kratzig, and M. Wohlecke, J. Appl. Phys. 93, 3102 (2003). 96 Y. H. Han and D. M. Smyth, Am. Ceram. Soc. Bull. 62, 852 (1983). 97 H. Donnerberg, S. M. Tomlinson, C. R. A. Catlow, and O. F. Schirmer, Phys. Rev. B 44, 4877 (1991). 98 N. Iyi, K. Kitamura, F. Izumi, J. K. Yamamoto, T. Hayashi, H. Asano, and S. Kimura, J. Solid State Chem. 101, 340 (1992). PAGE 185 185 99 N. Zotov, H. Boysen, F. Frey, T. Metzger, and E. Born, J. Phys. Chem. Solids 55, 145 (1994). 100 G. S. Zhdanov, S. A. Ivanov, E. V. Kolontsova, and A. E. Korneev, Ferroelectrics 21, 463 (1978). 101 S. A. Ivanov, A. E. Korneyev, E. V. Kolontsova, and Y. N. Venevtsev, Kristallog rafiya 23, 1071 (1978). 102 E. M. Ivanova, N. A. Sergeev, and A. V. Yatsenko, Crystallogr. Rep. 43, 303 (1998). 103 A. V. Yatsenko, E. N. Ivanova, and N. A. Sergeev, Physica B 240, 254 (1997). 104 T. W. D. Farley, W. Hayes, S. Hull, R. Ward, M. T. Hutching s, and M. Alba, Solid State Ionics 28, 189 (1988). 105 J. B. Goodenough and A. Hamnett, Physics of NonTetrahedrally Bonded Binary Compounds III (Springer Verlag, Berlin, 1984). 106 I. P. Zibrov, V. P. Filonenko, P. E. Werner, B. O. Marinder, and M. Sundbe rg, J. Solid State Chem. 141, 205 (1998). 107 F. Holtzberg, A. Reisman, M. Berry, and M. Berkenblit, J. Am. Chem. Soc. 78, 1536 (1956). 108 G. Kresse and J. Furthmuller, Comput. Mater. Sci. 6 15 (1996). 109 G. Kresse and J. Furthmuller, Phys. Rev. B 54, 1 1169 (1996). 110 P. Blaha, K. Schwarz, P. Sorantin, and S. B. Trickey, Comput. Phys. Commun. 59, 399 (1990). 111 D. Vanderbilt, Phys. Rev. B 41, 7892 (1990). 112 P. Pulay, Chem. Phys. Lett. 73, 393 (1980). 113 H. Boysen and F. Altorfer, Acta Crystallogr., Sect. B: Struct. Sci. 50, 405 (1994). 114 W. Y. Ching, Z. Q. Gu, and Y. N. Xu, Phys. Rev. B 50, 1992 (1994). 115 D. Redfield and W. J. Burke, J. Appl. Phys. 45, 4566 (1974). 116 A. Dhar and A. Mansingh, J. Appl Phys. 68, 5804 (1990). 117 H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976). 118 C. G. Van de Walle and J. Neugebauer, J. Appl. Phys. 95, 3851 (2004). PAGE 186 186 119 C. G. Vandewalle, P. J. H. Denteneer, Y. Baryam, and S. T. Pantelides, Phys. Rev. B 39, 10791 (1989). 120 S. G. Louie, M. Schluter, J. R. Chelikowsky, and M. L. Cohen, Phys. Rev. B 13, 1654 (1976). 121 W. E. Pickett, K. M. Ho, and M. L. Cohen, Phys. Rev. B 19, 1734 (1979). 122 A. F. Kohan, G. Ceder, D. Morgan, and C. G. Van de Walle, Phys. Rev. B 61, 15019 (2000). 123 G. Makov and M. C. Payne, Phys. Rev. B 51, 4014 (1995). 124 U. Gerstmann, E. Rauls, T. Frauenheim, and H. Overhof, Phys. Rev. B 67 (2003). 125 D. Segev and S. H. W ei, Phys. Rev. Lett. 91 (2003). 126 K. Reuter and M. Scheffler, Phys. Rev. B 65, 035406 (2002). 127 P. Erhart and K. Albe, J. Appl. Phys. 102 (2007). 128 T. Tanaka, K. Matsunaga, Y. Ikuhara, and T. Yamamoto, Phys. Rev. B 68 (2003). 129 R. DeHoff, Thermodynamics in Materials Science (Taylor & Francis Group, Boca Taton, FL, 2006). 130 J. Malcolm W. Chase, NIST JANAF thermochemical tables (American Chemical Society, Washington, D.C., 1998). 131 I. Barin, O. Knacke, and O. Kubaschewski, Thermochemical propertie s of inorganic substances: supplement (Springer Verlag, New York, 1977). 132 O. Knacke, O. Kabaschewski, and K. Hesselmann, Thermochemical properties of inorganic substances (Spring, Berlin, 1991). 133 J. He, R. K. Behera, M. W. Finnis, X. Li, E. C. Dickey S. R. Phillpot, and S. B. Sinnott, Acta Mater. 55, 4325 (2007). 134 P. W. Anderson, Phys. Rev. Lett. 34, 953 (1975). 135 S. Kim, V. Gopalan, K. Kitamura, and Y. Furukawa, J. Appl. Phys. 90, 2949 (2001). 136 H. X. Xu, D. Lee, J. He, S. B. Sinnott, V. Gopalan, V. Dierolf, and S. R. Phillpot, Phys. Rev. B 78, 174103 (2008). 137 J. H. Yao, J. J. Xu, G. Y. Zhang, X. J. Chen, B. Li, and G. G. Li, Chin. Phys. Lett. 17, 513 (2000). 138 R. L. Byer, J. Nonlinear Opt. P hys. Mater 6 549 (1997). PAGE 187 187 139 A. Mehta, E. K. Chang, and D. M. Smyth, J. Mater. Res. 6 851 (1991). 140 T. K. Halstead, J. Chem. Phys. 53, 3427 (1970). 141 N. Schmidt, K. Betzler, M. Grabs, S. Kapphan, and F. Klose, J. Appl. Phys. 65, 1253 (1989). 142 D. P Birnie, J. Mater. Sci. 28, 302 (1993). 143 V. Dierolf and C. Sandmann, J. Lumin. 125, 67 (2007). 144 J. Zou, S. M. Zhou, J. Xu, L. H. Zhang, Z. L. Xie, P. Han, and R. Zhang, J. Appl. Phys. 98 (2005). 145 H. Suche, R. Wessel, S. Westenhofer, W. Sohler, S. Bosso, C. Carmannini, and R. Corsini, Opt. Lett. 20, 596 (1995). 146 J. G. Sole, L. Bausa, D. Jaque, E. Montoya, H. Murrieta, and F. Jaque, Spectrochim. Acta, Part A 54, 1571 (1998). 147 V. Dierolf, C. Sandmann, A. B. Kutsenko, and T. Troster, Radiat. Eff Defects Solids 155, 253 (2001). 148 S. Balsamo, S. Maio, I. Montrosset, H. Suche, and W. Sohler, Opt. Quantum Electron. 31, 29 (1999). 149 I. Baumann, S. Bosso, R. Brinkmann, R. Corsini, M. Dinand, A. Greiner, K. Schafer, J. Sochtig, W. Sohler, H. Suche, and R. Wessel, IEEE J. Sel. Top. Quantum Electron. 2 355 (1996). 150 J. Amin, J. A. Aust, D. L. Veasey, and N. A. Sanford, Electron. Lett. 34, 456 (1998). 151 V. Voinot, R. Ferriere, and J. P. Goedgebuer, Electron. Lett. 34, 549 (1998). 152 L. Rebouta, M. F. daSilva, J. C. Soares, D. Serrano, E. Dieguez, F. AgulloLopez, and J. Tornero, Appl. Phys. Lett. 70, 1070 (1997). 153 R. M. Araujo, K. L engyel, R. A. Jackson, L. Kovacs, and M. E. G. Valerio, J. Phys.: Condens. Matter 19, 046211 (2007). 154 R. D. Shannon, Acta Crystallographica Section A 32, 751 (1976). 155 R. M. Moon, W. C. Koehler, H. R. Child, and Raubenhe.Lj, Phys. Rev. 1 722 (1968). 156 M. Faucher, J. Pannetier, Y. Charreire, and P. Caro, Acta Crystallogr., Sect. B: Struct. Sci. 38, 344 (1982). 157 J. S. Olsen, C. S. G. Cousins, L. Gerward, H. Jhans, and B. J. Sheldon, Phys. Scr. 43, 327 (1991). PAGE 188 188 158 C. T. Prewitt, R. D. Shannon, D. B. Rogers, and A. W. Sleight, Inorg. Chem. 8 1985 (1969). 159 A. M. Boring and J. L. Smith, Challenges in Plutonium Science I 91 (2000). 160 K. T. Moore and G. van der Laan, Rev. Mod. Phys. 81, 235 (2009). 161 U. Vonbarth and C. D. Gelatt, Physical Review B 21, 2222 (1980). 162 A. Kiejna, G. Kresse, J. Rogal, A. De Sarkar, K. Reuter, and M. Scheffler, Phys. Rev. B 73 (2006). 163 J. Hafner, J. Comput. Chem. 29, 2044 (2008). 164 T. Oguchi, K. Terakura, and A. R. Williams, Ph ys. Rev. B 28, 6443 (1983). 165 M. Cococcioni and S. de Gironcoli, Phys. Rev. B 71 (2005). 166 Y. Wang and D. J. Doren, Solid State Commun. 136, 186 (2005). 167 L. Wang, T. Maxisch, and G. Ceder, Phys. Rev. B 73, 195107 (2006). 168 R. Grau Crespo, F. Cora, A. A. Sokol, N. H. de Leeuw, and C. R. A. Catlow, Phys. Rev. B 73 (2006). 169 A. Rohrbach, J. Hafner, and G. Kresse, Phys. Rev. B 70 (2004). 170 R. Zimmermann, P. Steiner, R. Claessen, F. Reinert, S. Hufner, P. Blaha, and P. Dufek, J. Phys.: Condens. Matt er 11, 1657 (1999). 171 C. A. Mccammon and L. G. Liu, Physics and Chemistry of Minerals 10, 106 (1984). 172 J. Kubler and A. R. Williams, J. Magn. Magn. Mater. 547 603 (1986). 173 P. D. Battle and A. K. Cheetham, Journal of Physics C: Solid State Physics 12, 337 (1979). 174 C. N. R. R. a. B. Raveau, Transition Metal Oxixdes (Wiley VCH, New York, 1995). 175 E. N. Maslen, V. A. Streltsov, N. R. Streltsova, and N. Ishizawa, Acta Crystallogr., Sect. B: Struct. Sci. 50, 435 (1994). 176 G. Rollmann, A. Rohrbach, P. Entel, and J. Hafner, Phys. Rev. B 69 (2004). 177 B. Hourahine, S. Sanna, B. Aradi, C. Kohler, and T. Frauenheim, Physica B 376, 512 (2006). 178 A. Lazreg, Z. Dridi, F. Benkabou, and B. Bouhafs, Physica B 403, 2702 (2008). PAGE 189 189 179 M. Losurdo, M. M. Giangregorio, P. Capezzuto, G. Bruno, R. G. Toro, G. Malandrino, I. L. Fragala, L. Armelao, D. Barreca, E. Tondello, A. A. Suvorova, D. Yang, and E. A. Irene, Adv. Funct. Mater. 17, 3607 (2007). 180 J. P. McKelvey, Solid State Physics For Engineering and Materials Science (Krieger Publishing Company, Malabar, Florida, 1993). 181 M. Carrascosa and F. Agullolopez, Journal of the Optical Society of America B Optical Physics 7 2317 (1990). 182 A. Harhira, L. Guilbert, P. Bourson, and H. Rinne rt, Applied Physics B Lasers and Optics 92, 555 (2008). 183 F. Jermann and J. Otten, Journal of the Optical Society of America B Optical Physics 10, 2085 (1993). 184 I. W. Kim, S. S. Yi, V. F. Pichugin, V. Y. Yakovlev, and M. S. Dmitriev, J. Cryst. Growth 253, 319 (2003). 185 J. E. Alfonso, M. J. Martin, and C. Zaldo, Appl. Phys. Lett. 71, 2904 (1997). 186 L. S. Qiang, H. X. Zhang, and C. Q. Xu, Mater. Chem. Phys. 77, 6 (2003). 187 D. L. Zhang, D. C. Wang, and E. Y. B. Pun, J. Appl. Phys. 97 (2005). 188 D. M. B. P. Milori, I. J. Moraes, A. C. Hernandes, R. R. Desouza, M. S. Li, M. C. Terrile, and G. E. Barberis, Phys. Rev. B 52, 3802 (1995). 189 K. Polgar and A. P. Skvortsov, Opt. Spektrosk. 58, 229 (1985). 190 J. Blumel, E. Born, and T. Metzger, J. Phys. Ch em. Solids 55, 589 (1994). 191 S. Y. Wu and W. C. Zheng, Phys. Rev. B 65 (2002). 192 C. Buchal and S. Mohr, J. Mater. Res. 6 134 (1991). 193 V. B. Ptashnik, T. Y. Dunaeva, and I. V. Myasnikov, Inorg. Mater. 21, 1814 (1985). 194 V. I. Lapshin and A. P. Rum yantsev, Inorg. Mater. 12, 1797 (1976). 195 G. I. Malovichko, V. G. Grachev, and O. F. Schirmer, Solid State Commun. 89, 195 (1994). 196 H. Wang, X. Y. Kuang, D. Dong, Y. Xiong, and K. W. Zhou, Physica B 367, 53 (2005). 197 T. Gog, P. Schotters, J. Falta, G. Materlik, and M. Grodzicki, J. Phys.: Condens. Matter 7 6971 (1995). 198 D. S. Yang, N. Sung, and T. H. Yeom, J. Phys. Soc. Jpn. 78, 114605 (2009). PAGE 190 190 199 G. I. Malovichko, V. G. Grachev, O. F. Schirmer, and B. Faust, J. Phys.: Condens. Matter 5 3971 (19 93). 200 G. E. Peterson, A. M. Glass, and T. J. Negran, Appl. Phys. Lett. 19, 130 (1971). 201 S. Omar, in Department of Materials Science and Engineering (University of Florida, Gainesville, 2008), Vol. Ph.D. 202 N. E. Brese and M. Okeeffe, Acta Crystallographica Section B Structural Science 47, 192 (1991). 203 L.Pauling, The Nature of the Chemical Bond (Cornell University Press, New York, 1960). 204 I. Syozi, Progress of Theoretical Physics 6 306 (1951). 205 P. G Watson, Physica 75, 627 (1974). 206 M. Yashima, S. Kobayashi, and T. Yasui, Solid State Ionics 177, 211 (2006). 207 I. U. o. Crystallography, The International Tables of Crystallography (Kluwer Academic Publishers, 1996). 208 H. Hayashi, R. Sagawa, H. Inaba, and K. Kawamura, Solid State Ionics 131, 281 (2000). 209 A. Atkinson and A. Selcuk, Solid State Ionics 134, 59 (2000). 210 J. D. Gale, Journal of the Chemical Society Faraday Transactions 93, 629 (1997). 211 J. D. Gale and A. L. Rohl, Molecular Simulation 29, 291 (2003). 212 R. W. Grimes, G. Busker, M. A. McCoy, A. Chroneos, J. A. Kilner, and S. P. Chen, Berichte Der Bunsen Gesellschaft Physical Chemistry Chemical Physics 101, 1204 (1997). 213 S. Vyas, R. W. Grimes, D. H. Ga y, and A. L. Rohl, Journal of the Chemical SocietyFaraday Transactions 94, 427 (1998). 214 V. Butler, C. R. A. Catlow, B. E. F. Fender, and J. H. Harding, Solid State Ionics 8 109 (1983). 215 A. Gotte, K. Hermansson, and M. Baudin, Surface Science 552, 2 73 (2004). 216 A. Gotte, D. Spangberg, K. Hermansson, and M. Baudin, Solid State Ionics 178, 1421 (2007). 217 R. W. Grimes, G. Busker, M. A. McCoy, A. Chroneos, J. A. Kilner, and S. P. Chen, Physical Chemistry Chemical Physics 101, 1204 (1997). 218 R. W. G rimes, D. J. Binks, and A. B. Lidiard, Philosophical Magazine A 72, 651 (1995). PAGE 191 191 219 M. A. McCoy, R. W. Grimes, and W. E. Lee, Philosophical Magazine A 76, 1187 (1997). 220 H. Inaba, R. Sagawa, H. Hayashi, and K. Kawamura, Solid State Ionics 122, 95 (1999). 221 J. R. Sims and R. N. Blumenthal, High Temperature Science, 99 (1976). 222 H. W. Chiang, R. N. Blumenthal, and R. A. Fournelle, Solid State Ionics 66, 85 (1993). 223 J. R. Sims and R. N. Blumenthal, High Temperature Science 8 99 (1976). 224 A. Nakajim a, A. Yoshihara, and M. Ishigame, Physical Review B 50, 13297 (1994). 225 L. Gerward, J. S. Olsen, L. Petit, G. Vaitheeswaran, V. Kanchana, and A. Svane, Journal of Alloys and Compounds 400, 56 (2005). 226 Y. L. Wang, K. Duncan, E. D. Wachsman, and F. Ebrahimi, Solid State Ionics 178, 53 (2007). 227 V. Kanchana, G. Vaitheeswaran, A. Svane, and A. Delin, Journal of Physics Condensed Matter 18, 9615 (2006). 228 S. J. Duclos, Y. K. Vohra, A. L. Ruoff, A. Jayara man, and G. P. Espinosa, Physical Review B 38, 7755 (1988). 229 D. J. Green, An introduction to the mechanical properties of ceramics (Cambridge University Press, 2004). 230 D. J. M. Bevan and J. Kordis, Journal of Inorganic & Nuclear Chemistry 26, 1509 (1964). 231 R. J. Panlener, R. N. Blumenthal, and J. E. Garnier, Journal of Physics and Chemistry of Solids 36, 1213 (1975). 232 S. Sameshima, H. Ono, K. Higashi, K. Sonoda, Y. Hirata, and Y. Ikuma, Journal of the Ceramic Society of Japan 108, 1060 (2000). 2 33 N. Sammes, G. Tompsett, Y. J. Zhang, A. Cartner, and R. Torrens, Denki Kagaku 64, 674 (1996). 234 M. W. Barsoum, Fundamentals of Ceramics (McGraw Hill, 1997). 235 M. Katsuragawa, H. Kashihara, and M. Akebi, J. Nucl. Mater. 204, 14 (1993). 236 R. N. Hill, D. C. Wade, J. R. Liaw, and E. K. Fujita, Nucl. Sci. Eng. 121, 17 (1995). 237 T. Inoue, M. Sakata, H. Miyashiro, T. Matsumura, A. Sasahara, and N. Yoshiki, Nucl. Technol. 93, 206 (1991). 238 B. Cox and Y. M. Wong, J. Nucl. Mater. 270, 134 (1999). PAGE 192 192 239 S. I. Hong, K. W. Lee, and K. T. Kim, J. Nucl. Mater. 303, 169 (2002). 240 E. Vitikainen and P. Nenonen, J. Nucl. Mater. 78, 362 (1978). 241 R. E. Stoller, G. R. Odette, and B. D. Wirth, J. Nucl. Mater. 251, 49 (1997). 242 J. Gan, G. S. Was, and R. E. Stoller, J. Nucl. Mater. 299, 53 (2001). 243 R. E. Stoller, J. Nucl. Mater. 276, 22 (2000). 244 R. E. Stoller and A. F. Calder, J. Nucl. Mater. 283, 746 (2000). 245 S. J. Wooding, D. J. Bacon, and W. J. Phythian, Philosophical Magazine a Physics of Cond ensed Matter Structure Defects and Mechanical Properties 72, 1261 (1995). 246 Y. M. Kim, B. J. Lee, and M. I. Baskes, Phys. Rev. B 74 (2006). 247 J. F. Ziegler, J. P. Biersak, and U. Littmark, The Stooping and Range of Ions in Matter Pergomon, New Youk, 1985). 248 D. R. Trinkle, (The Ohio State University, 2003), Vol. Ph.D. 249 M. I. Baskes, Mater. Chem. Phys. 50, 152 (1997). 250 J. H. Rose, J. R. Smith, F. Guinea, and J. Ferrante, Phys. Rev. B 29, 2963 (1984). 251 B. J. Lee and M. I. Baskes, Phys. Rev. B 62, 8564 (2000). 252 B. J. Lee, M. I. Baskes, H. Kim, and Y. K. Cho, Phys. Rev. B 6418 (2001). 253 B. J. Lee, J. H. Shim, and M. I. Baskes, Phys. Rev. B 68 (2003). 254 F. H. Stillinger and T. A. Weber, Phys. Rev. B 31, 5262 (1985). 255 R. A. Johnson, Phi losophical Magazine a Physics of Condensed Matter Structure Defects and Mechanical Properties 63, 865 (1991). 256 G. J. Ackland, Philosophical Magazine a Physics of Condensed Matter Structure Defects and Mechanical Properties 66, 917 (1992). 257 E. A. Bran des, Smithells Metals Reference Book (Butterworths, London, 1983). 258 C. S. Barrett and T. B. Massalski, Structure of Metals (McGraw Hill, New York, 1966). 259 G. Simmons and H. Wang, Single Crystal Elastic Constants and Calculated Aggregate Properties: A Handbook (MIT Press, Cambridge, 1971). 260 A. T. Dinsdale, Calphad Computer Coupling of Phase Diagrams and Thermochemistry 15, 317 (1991). PAGE 193 193 261 F. R. de Boer, W. C. M. Mattens, A. R. Miedema, and A. K. Niessen, Cohesion in Metals: Transition Metal Alloys (Horth Holland, Amsterdam, 1988). 262 R. Smith, Atomic & Ion Collisions in Solids and at Surfaces (Cambridge University Press, Cambridge, UK, 1997). 263 D. Faken and H. Jonsson, Comput. Mater. Sci. 2 279 (1994). 264 A. Meldrum, L. A. Boatner, and R. C. Ewi ng, Nuclear Instruments & Methods in Physics Research Section B Beam Interactions with Materials and Atoms 207, 28 (2003). PAGE 194 194 BIOGRAPHICAL SKETCH Haixuan Xu was born in Anshan, Liaoning Province, China. He received his Bachelor of Engineering degree in Metallic Materials Engineering from Dalian University of Technology, Dalian, in 2005. He was also a member of Institute of University Students Innova tion of the university. Then he came to the United States of America for higher education. He joined Computational Materials Science Focus Group (CMSFG) in the Department of Materials and Engineering at University of Florida, with Prof. Simon R. Phillpot as his advisor. He received Master of Science from University of Florida at 2007. He was awarded the Doctoral Degree in Materials Science and Engineering in May 2010. 