UFDC Home  myUFDC Home  Help 



Full Text  
PAGE 1 1 EFFECT OF SURFACES, DOMAIN WALL S AND GRAIN BOUNDARIES ON FERRO ELECTRICITY IN LEAD TITANATE USING ATOMIC SCALE SIMULATIONS By RAKESH KUMAR BEHERA 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 2009 PAGE 2 2 2009 Rakesh Kumar Behera PAGE 3 3 To my parents PAGE 4 4 ACKNOWLEDGMENTS I would like to express my since re gratitude to my advisor Prof. Simon R. Phillpot for his invaluable guidance and encouragement extended throughout the study. His tenacious supervision, helpful suggestion, patience and time deserve a special attention. His advice on any topic mixed with extraordinary humor made my time at UF extremely special. I admire him as a great researcher, teacher, mentor and motivator I would like to express my appreciation to my committee members Prof. Susan B. Sinnott, Prof. David Norton, Prof. Juan C. Nino and Prof. Aravind Asthagiri for their continuous support and advice. I would like to convey my special thanks to my masters advisor Prof. Dorel Moldovan at LSU, who recommended me to work w ith this wonderful group at UF and Dr. Venkatraman Gopalan of PSU for his guidance in research. I was fortunate enough to work with a great team of motivated students and Post docs in the Computational Materials Science Focus Group. The curricular and extracurricular activities in the group made me experience a wide variet y of research topics and simultaneously enjoy my life to the fullest. I would like to specifically thank Prof. Marina Yao; Dr s Tao Liang, Boris Ni, Alex Chernatinsky, Wen Dung Hsu, SeongJun Heo, Jun He, Taku Watanabe, Sharon Pregler Dilpuneet Aidhy, Pank a j Nerikar, Saurav Pathak, Jennifer Wohlwend, Patrick Chiu, Sankara Tatiparti, Shobit Omar, Abhijit Par manick; Bryce Devine, Peter Barry, Donghwa Lee, Priyank Shukla, Haixuan Xu, Jasmine Crenshaw, Fernanda Foertter, Chanwoo Lee, D onghyun Kim, Mahesh Tanni ru, Krishna Ganeshan, Ray Shan and Travis Kemper. I am also thankful to my friends at LSU and India for their const ant support and encouragement Finally, I would like to express my gratitude towards my parents, brothers, sisters and their families for th eir unconditional love and s upport They have been a constant source of PAGE 5 5 i nspiration during my stay at UF This work was supported by NSF ITR grant number DMR 0426870. PAGE 6 6 TABLE OF CONTENTS page ACKNOWLEDGMENTS .................................................................................................................... 4 LIST OF TABLES ................................................................................................................................ 9 LIST OF FIGURES ............................................................................................................................ 10 ABSTRACT ........................................................................................................................................ 16 CHAPTER 1. GENERAL INTRODUCTION .................................................................................................. 18 1.1 What is Ferroelectricity? ...................................................................................................... 18 1.2 Effect of Electric Field in Ionic Crystals ............................................................................. 19 1.3 Origin of Spontaneous Polarization ..................................................................................... 21 1.4 Basics of Ferroelect ricity in Perovskites ............................................................................. 27 1.5 Theoretical Approach to Describe Ferroelectricity ............................................................. 30 1.6 Applications of Ferroelectric Materials ............................................................................... 31 1.7 Motivation for This Work .................................................................................................... 32 1.7.1 Ferroelectric Memories .............................................................................................. 32 1.7.2 Possible Future Applications of Ferroelectric Materials ......................................... 34 1.8 Objective of This Work ........................................................................................................ 35 2 An OVERVIEW of PEROVSKITE STRU CTURED MATERIALS ..................................... 36 2.1 Introduction to Perovskites ................................................................................................... 36 2.2 Perovskite Oxides ................................................................................................................. 37 2.2 The Use of Perovskite Oxides .............................................................................................. 39 2.2.1 Structure Field Maps .................................................................................................. 39 2.2.2 Ferroelectricity ............................................................................................................ 40 2.3 PbTiO3 ................................................................................................................................... 41 2.3.1 Crystallography of PbTiO3 ........................................................................................ 41 2.3.2 Structural Details for Cubic PT ................................................................................. 42 2.3.3 Structural Details for Tetragonal PT ......................................................................... 43 3 SIMULATION METHODS ....................................................................................................... 46 3.1 Periodic Boundary Conditions ............................................................................................. 47 3.2 Simulation Methods Used in This Study ............................................................................. 48 3.2.1 First Principles Calculations ...................................................................................... 48 3.2.2 Molecular Dynamics (MD) Simulation .................................................................... 57 3.3 Model Validation (Bulk PbTiO3) ......................................................................................... 72 PAGE 7 7 4 EFFECT OF SURFACES ON FERROELECTRICITY .......................................................... 77 4.1 Introduction ........................................................................................................................... 77 4.2 (001) Surface Terminations in PbTiO3 ................................................................................ 78 4.3 Critical Issues ........................................................................................................................ 81 4.3.1 Use of Fixed Charge Simulations for Surfaces ........................................................ 82 4.3.2 Short circuit vs. Open circuit Electrical Boundary Condition ................................ 83 4.3.3 Dipole Effect in Symmetric and Asymmetric Slabs ................................................ 84 4.4. Characterization of Surfaces ............................................................................................... 86 4.4.1 Atomic Displacements ............................................................................................... 86 4.4.2 Surface Rumpling and Interlayer Distance ............................................................... 87 4.4.3 Surface Polarization ................................................................................................... 89 4.4.4 Surface Energy ........................................................................................................... 94 4.5. Di scussion on Stability of Free Standing Films ................................................................. 98 4.5.1 Normal Polarization (InOut Polarization) ............................................................... 98 4.5.2 Tail Tail Polarization (Out Out Polarization) .......................................................... 99 4.5.3 Parallel Polarization ................................................................................................. 100 4.6. Conclusions ........................................................................................................................ 101 4.7. Future Work ....................................................................................................................... 102 5 EFFECT OF (100) 180 DOMAIN WALLS IN LEAD TITANATE ................................... 103 5.1 Introduction ......................................................................................................................... 103 5.2 Types of Domain Walls in PbTiO3 .................................................................................... 103 5.2.1 Pbcentered and Ti centered 180 Domain Walls .................................................. 105 5.2.2 Domain Wall Set up ................................................................................................. 107 5.3 Simulation Methodology .................................................................................................... 107 5.4 Motivation for the Domain Wall Simulations ................................................................... 108 5.5 Structural Optimization Schemes ....................................................................................... 110 5.6 Model Validation (Characterization of Domain Walls with R Z Scheme) ..................... 111 5.6.1 Domain Wall Energies ............................................................................................. 111 5.6.2 Atomic Displacements ............................................................................................. 112 5.6.3 Ferroel ectric Distortion ............................................................................................ 113 5.6.4 Summary of USPP Calculations .............................................................................. 115 5.7 Characterization of Domain Walls with PAW Potential (R Z Scheme) ......................... 116 5.7.1 Domain Wall Energies ............................................................................................. 116 5.7.2 Polarization Calculation ........................................................................................... 117 5.7.3 Polarization Fitting to Calculate Domain Width .................................................... 118 5.8 Characterization of Domain Walls with PAW Potential (R XYZ Scheme) ................... 119 5.8.1 Domain Wall Energies ............................................................................................. 119 5.8.2 Polarization Calculation ........................................................................................... 120 5.8.3 Atomic Displacements ............................................................................................. 125 5.9 Characterization of Domain Walls with PAW Potential (A XYZ Scheme) ................... 127 5.10 Discussion ......................................................................................................................... 128 5.10.1 Domain Walls Energy Comparison between (100) and (110) Domain Walls ... 129 5.10.2 Polarization Variation Comparison between (100) and (110) Domain Walls ... 1 30 5.11 Critical Assessment of the Inter atomic Potentials ......................................................... 130 PAGE 8 8 5.11.1 Domain Wall Energies ........................................................................................... 131 5.11.2 Polarization Analysis ............................................................................................. 132 5.12 Conclusions ....................................................................................................................... 133 5.13 Future Work ...................................................................................................................... 133 6 STRUCTURE AND PROPERTIES OF 5 (310) TILT GRAIN BOUNDARIES IN STRONTIUM TITANATE AND LEAD TITANATE .......................................................... 134 6.1 Effect of Grain Boundaries on Ferroelectricity ................................................................. 134 6.2 Coincident Site Lattice (CSL) ............................................................................................ 136 6.3 Bicrystal Simulation Setup ................................................................................................. 138 6.3 .1 Geometry of the (310) Planes in SrTiO3 ................................................................. 139 3 ...................................................... 140 6.4 Characterization of Bulk SrTiO3 ........................................................................................ 141 6.5 Characterization of Cubic SrTiO3 Grain Boundary .......................................................... 145 6.5.1 R Z Relaxation Scheme ........................................................................................... 146 6.5.2 A XYZ Relaxation Scheme ..................................................................................... 147 6.5.3 Structural Analysis ................................................................................................... 149 6.6 Characterization of AFD SrTiO3 Grain Boundary ............................................................ 152 6.6.1 Energetics of G Bs ..................................................................................................... 152 6.6.2 Polarization and Structure Analysis ........................................................................ 152 3 ............................................................................... 156 6.7.1 R Z Relaxation Scheme ........................................................................................... 157 6.7.2 A XYZ Relaxation Scheme ..................................................................................... 160 6.8 Conclusions ......................................................................................................................... 162 6.9 Future Work ........................................................................................................................ 162 7 PARALLELIZATION OF THE SHELL MODEL MOLECULAR DYNAMICS CODE .. 163 7.1 Introduction ......................................................................................................................... 163 7.2. Parallel Computing ............................................................................................................ 165 7.2.1. Link Cell Algorithm and Spatial Decomposition ................................................. 165 7.2.2 Dynamic Memory Allocation .................................................................................. 167 7.3. Energy Comparison ........................................................................................................... 168 7.4. Perfo rmance of Parallel Code ........................................................................................... 169 8 GENERAL CONCLUSIONS .................................................................................................. 174 8.1. Surfaces .............................................................................................................................. 174 8.2. Domain Walls ..................................................................................................................... 175 8.3. Grain Boundaries ............................................................................................................... 177 LIST OF REFERENCES ................................................................................................................. 180 BIOGRAPHICAL SKETCH ........................................................................................................... 193 PAGE 9 9 LIST OF TABLES Table page 3 1 Potential parameters used for describing interatomic interactions ..................................... 62 3 2 Summary of the eight thermodynamical ensembles and their corresponding fixed and dependent variables.. .............................................................................................................. 69 3 3 Quantitative analysis of bulk PT predicted from both interatomic potentials. ................... 75 4 1 List of simulation details for PT thin films ........................................................................... 77 4 2 Surface characte rization for PbO terminations using DFT (in % of c lattice parameter) ............................................................................................................................... 85 4 3 Surface characterizations for TiO2terminations using DFT (in % of c lattice parameter) ............................................................................................................................... 86 5 1 Details of DFT simulation approaches for PT domain wall structures ............................. 103 5 2 Calculated Born effective charges of the optimized Cubic and Tetragonal P T structures as calculated with a PAW potential, and compared with literature values. .... 117 7 1 Comparison of energy/atom for core and shell of each atom type in PT between the existing serial c ode and modified parallel code ................................................................. 168 7 2 Comparison of energy/atom of each atom type in PT by turning off the shell contributions between the existing serial code and modified parallel code ..................... 169 PAGE 10 10 LIST OF FIGURES Figure page 1 1 Schematic of different sources of electric dipoles in an ionic crystal ................................. 20 1 2 Schematic of frequency dispersion behavior of a typical ionic crystal ............................. 20 1 3 Schematic of t he energy contribution due to dipole dipole coupling, elastic energy, and the total energy to explain the origin of ferroelectricity ............................................... 23 1 4 Schematic representation of some of the possible eigen lattice vibration mod es in a perovskite material ................................................................................................................ 25 1 5 Phonon disp ersion curves calculated for cubic, and tetragonal PbTiO3 ............................. 26 1 6 Schematic of the unit cell of BaTiO3 .................................................................................... 27 1 7 Phase transition and spontaneous polarization variation with temperature in BaTiO3. ..... 28 1 8 Schematic of surface charge associated with spontaneous polarization, and formation of 180 domains to minimi ze the electrostatic energy ...................................... 29 1 9 Hysteresis behavior of the polarization, in relation to the applied electric field f or a ferroelectric crystal ................................................................................................................. 30 1 10 Schematic of the temperature dependence of spontaneous polarization, permittivity and inverse permittivity in a ferroelectric material .............................................................. 32 1 11 The development of FeRAM technology with respect to the des ign rule and memory density, and t he memory capacity roadmap for a standard FeRAM ................................... 33 2 1 Schematic repr esentation of an ideal ABX3 perovskite structure ....................................... 36 2 2 Summarized list of various applications of perovskite materials emphasizing different functionality. ........................................................................................................... 37 2 3 Schematic representation of ABO3 perovskite structure. .................................................... 38 2 4 Lattice energy contour maps for predicting the stability of different compounds with A+3B+3O3 2 structures ............................................................................................................ 40 2 5 Schematic of the tetragonal perovskite crystal structure of PT. .......................................... 41 2 6 Schematics showing three different s ymmetry ope rations represented by cubic PT ......... 42 2 7 Schematic two dimensional representation of P4mm symmet ry operation. ...................... 43 2 8 Sc hematic representation of the break of high symmetry in tetragonal PT. ....................... 44 PAGE 11 11 3 1 Sc hematic of a 2D periodic system ...................................................................................... 47 3 2 2D Schem atic of the charge neutralization schem e used in the direct sum method .......... 59 3 3 Illustation of local minima and maxima of a function in an interval .................................. 65 3 4 Schematic showing the golden search approach. ................................................................. 66 3 5 Schematic of steepest descent method. ................................................................................. 67 3 6 Schem atic of a 2D simulation system showing the piston setup to maintain pressure. ..... 70 3 7 Lattice parameter as a function of temperature for PT. ....................................................... 73 3 8 Polarization as a function of temperature for PT obtained with both the potentials. ......... 74 4 1 Schematic of all four possible free st anding ferroelectric (001) PT films .......................... 79 4 2 <010> view of P bO terminated PT with vacuum, change in the interlayer d istance and C) surface rumpling ........................................................................................................ 81 4 3 Comparison of the charge variation in each ion in surface and bulk with Bader charge analysis ....................................................................................................................... 82 4 4 Cell bycell atomic displacement of each individual atom type from their original relaxed polarized condition with Potential 2 ........................................................................ 87 4 5 Interlayer distance and surface rumpling calculated for ferroelectric PT with surface PbO In termination and PbO Out termination. .................................................................... 88 4 6 Interlayer distance and surface rumpling calculated for fe rroelectric PT with surface, TiO2In termination, and TiO2Out termination ................................................................... 88 4 7 Schematic of a Ti centered unit cell, and Pb centered unit cell used for polarization calculation. .............................................................................................................................. 89 4 8 Cell bycell c direction polarization obtained for a 6 x 6 x 32.5 tetragonal PT system with Potential 1 and Potential 2. ............................................................................................ 91 4 9 Polarization induced due to surface rumpling for tetragonal PT with no ferroelectricity for P bO termination and TiO2termination ................................................. 93 4 10 Comparison of average Inand Out polarization surface energies for PbO and TiO2 termination in tetragonal PT. ................................................................................................. 95 4 11 Cell bycell excess energy per unit ce ll compared to the bulk PT energy for asymmetric slabs calculated with Potential 2. ...................................................................... 96 PAGE 12 12 4 12 Individual In and Out polarization surface energies predicted from both the empirical potentials. ............................................................................................................... 97 4 13 Schematic of asymmetric Tail Tail (Out Out ) polarization condition, and asymmetric parallel polarization condition in free standing ferroelectric films ..................................... 99 4 14 Total surface energy of PbO TiO2 terminated PT thin film for different termination polarization combinations with empirical potentials and DFT simulations ..................... 100 5 1 Schematic representation of the generation of 90 and 180 domain walls during the cubic to tetragonal phase transition below Curie temperature. ......................................... 104 5 2 Schematic representation of Pb centered, and Ti centered (100) 180 domain wall. The middle plane in each schematic indicates the domain wall. ...................................... 105 5 3 Schematic of (010) oriented Pb centered 180 domain wall showing the positions of the Pb atoms ......................................................................................................................... 106 5 4 Schematic representation of Ising type, Bloch type, and Nel type walls ........................ 109 5 5 Comparison of domain wall energies of Pband Ti centered domain walls using USPP with R Z scheme ....................................................................................................... 112 5 6 Relative z displacement of atoms with respect to the Pb atom in the Pb centered domain wall as a fu nction of the distance from the domain wall ...................................... 113 5 7 Schematic of tetragonal ferroelectric PT unit cell highlighting the Ti bonding .............. 114 5 8 Ferroelectric lattice distortion of TiO2planes in the z direction for a Pb centered 180 domain wall ................................................................................................................. 115 5 9 Comparison of domain wall energies of Pband Ti centered domain walls using PAW potential with R Z scheme ........................................................................................ 116 5 10 Cell bycell z polarization variation across the Pb centered domain wall for differen t system sizes with PAW pseudo potential ........................................................................... 118 5 11 Comparison of domain wall energies of Pband Ti centered domain walls using PAW potential with R Z and R XYZ schemes. ................................................................. 119 5 12 Cell bycell z polarization variation across the Pb centered domain wall for different system sizes with PAW pseudo potential within R XYZ scheme. ................................... 120 5 13 Cell bycell y polarization (Pn) variation acros s the Pb centered domain wall for different system sizes with PAW pseudo potential within R XYZ scheme .................... 121 5 14 Cell bycell normalized polarization fitting performed for both Pz and Pn for Pb centered domain wall using R XYZ scheme ...................................................................... 124 PAGE 13 13 5 15 Cell bycell normalized polarization fitting performed for both Pz and Pn for Ti centered domain wall using R XYZ scheme ...................................................................... 125 5 16 Schematic of the atomic displacement of each type of ion around the Pb centered domain wall .......................................................................................................................... 126 5 17 Change in tetragonality with increase i n system size for 180 Pbcentered domain wall systems using A XYZ scheme .................................................................................... 127 5 18 Change in Pn using R XYZ and A XYZ are the the same for Pb c entered domain wall ........................................................................................................................................ 128 5 19 Comparison of domain wall energies of (100) and (110) domain structures using R XYZ scheme ......................................................................................................................... 129 5 20 Cell bycell z polarization variation across the (100) Pbcentered and (1 10) OO centered domain wall with PAW pseudopotential within R XYZ scheme. .................... 130 5 21 Comparison of domain wall energies of (100) domain walls using Potential 1 and Pot ential 2 ............................................................................................................................. 131 5 22 Cell bycell polarization normal to the domain wall (Pn) predicted with Potential 1 and Potential 2 ...................................................................................................................... 132 6 1 Sc hematic of the relative rotati on of the grains to generate twist and symmetrical tilt bounda ries ............................................................................................................................. 135 6 2 Computed Egb as a function of rotation angle of symmetrical < 110> tilt boundaries in Al .......................................................................................................................................... 136 6 3 Schematic of a 36. tilt boundary and twist boundary ....... 137 6 4 Schematic of a bicrystal simulation set up ......................................................................... 138 6 5 Schemat ic of the two types of (310) surface terminations ................................................ 139 6 6 Schematic of the four possible twinned structures for (310) twist boundary in SrTiO3 .. 140 6 7 Prediction of cubic lattice parameter of SrTiO3 using various interatomic interactions 141 6 8 Comparison of energy grain due to the formation of the AFD phase in SrTiO3 for Akhtar and Tinte potentials with first principle calculations ............................................ 143 6 9 Compar ison of a lattice parameter and c/a ratio for the AFD phase in SrTiO3 ............... 143 6 10 Comparison of angle of rotation of the octahedral during AFD in SrTiO3. ..................... 144 6 11 AFD order parameter comparison between experiment, firs t principles calculation s and interatomic potential predictions in SrTiO3. ................................................................ 144 PAGE 14 14 6 12 Egb for different twinned STGBs with respect to system size for Grimes, Sekiguchi, and Thomas potentials in the R Z scheme. ......................................................................... 146 6 13 Relaxed Egb predicted from four different potential descriptions for a fully dense (310) STGB in SrTiO3 ......................................................................................................... 147 6 14 Relaxed GB str uctures with Sekiguchi potential for cation twin, anion twin, reverse cation twin, and reverse anion twin. ................................................................................... 149 6 15 Cell bycell polarization variation for anion twin GB using Sekiguchi potent ial............ 151 6 16 Cell bycell polarization variation for the reverse anion twin GB using Sekiguchi potential ............................................................................................................................... 151 6 17 Relaxed Egb p redicted for AFD SrTiO3 for a fully dense (310) STGB ............................ 152 6 18 Cell bycell polarization variation for the reverse anion twin GB using Akhtar potentia l ............................................................................................................................... 153 6 19 Structure of a relaxed reversed anion twin in (310) STGB of SrTiO3 with Akhtar potential showing octahedral rotation ................................................................................. 154 6 20 Cell bycell in plane polarizatio n variation for the reverse anion twin GB using Tinte potential ................................................................................................................................ 155 6 21 Structure of a relaxed reversed anion twin in (310) STGB of SrTiO3 with Tinte potential showing octahedral rotation ................................................................................. 156 6 22 Relaxed Egb for (310) STGBs in PT wi th R Z scheme using Potential 2 ......................... 158 6 23 Cell bycell in plane polarization variation for the Up Up reverse anion twin in PT using Potential 2 ................................................................................................................... 158 6 24 Cell bycell in plane polarization variation for the Up Dn reverse anion twin in PT using Potential 2 ................................................................................................................... 159 6 25 Comparison of Egb with R Z and A XYZ relaxation schemes for (310) STGBs in PT with Up Up polarization using Potential 2. ........................................................................ 160 6 26 Cell bycell polarization variation for the reverse anion twin using Potential 2 with A XYZ optimization scheme .............................................................................................. 161 7 1 The effect of the total number of atoms in the simulation system on the total time fo r 1000 MD runs for a bulk PT system using Potential 2 using serial code ......................... 164 7 2 Schematic repre sentation of 1 D spatial decomposition and 2D spatial decomposition .. 166 7 3 2D spatial decomposition of the 6 x 24 x 60 system ......................................................... 170 PAGE 15 15 7 4 Total simulation CPU time for running 1000 MD steps with number of processors for a 2D spat ial decomposition ............................................................................................ 171 7 5 Total simulation CPU time for running 1000 MD steps with number of processors for a 2D spatial decomposition with 2160 atoms in each processor ................................. 172 7 6 2D spatial decomposition of the 6 x 24 x 48 system with 16 processors ......................... 173 PAGE 16 16 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 EFFECT OF SURFACES, DOMAIN WALL S AND GRAIN BOUNDARIES ON FERRO ELECTRICITY IN LEAD TITANATE USING ATOMIC SCALE SIMULATIONS By Rakesh Kumar Behera December 2009 Chair: Simon R. Phillpot Major: Materials Science and Engineering Ultrathin ferroelectric (FE) films have potential application in non volatile random access memories and microelectromechanical systems. Miniaturization of such devices demands the establishment of the thickness property relationship, a detailed understanding of the surface termination, and the effects of domain walls and grain boundaries Density functional theory (DFT) and molecular dynamics (MD) methods are use d to investigate these effects i n lead titanate. The surface effects are characterized for the (100) terminations with open circuit electrical boundary condition. The resulting surface energies indicate that the surface terminations with polarization out of the surface are more stable than the cases where polarization occur into the surface. Analysis of atomic relaxation, surface rumpling and interlayer distance provide insight into the surface effects for each termination including polarization. The effect of domain walls are carried out for (100) and (110) 180 degree domains The (100) Pbcentered and (110) OO centered domain walls are predicted to be the favorable structure with comparable energies. The examination of polarization variation across the domain wall results in the presence of in plane polarization, which develops a small polarization rotation PAGE 17 17 analogous to that observed experimentally in ferromagnetic materials. The in plane polarization is perpendicular to the domain wall and points away from the wall. The effect of grain b oundaries on polarization are calculated for sigma 5 coincident site lattice (310) tilt boundaries in strontium titanate and lead titanate T he reverse twin ned structures are found to be energetically more favorable than the regular twins. The presence of the boundary resulted in a local polarization around the boundary (normal and parallel to the boundary), which vanishes away from the boundary. This study also focuses on the challenge to simulate complex ferroelectric structures. As a first step towards achieving the goal, the existing parallel molecular dynamics code was modified to include the core shell model, thereby allowing complex ferroelectric systems to be simulated. PAGE 18 18 CHAPTER 1 GENERAL INTRODUCTION Ferroelectricity was discovered in 1920 by J. V alasek [1, 2] Ever since its discovery, there has been a continuous effort to understand the phenomenon by researchers all over the world. These investigations have resulted in the development of new classes of devices. For example, expitaxially grown ultrathin films of ferroelectric materials, such as BaTiO3 (BT) and PbTiO3 (PT), have attract ed considerable attention due to their numerous potential applications [3 5] including non volatile memory components [6] The continuous miniaturization of such electronic devices demands an atomic level understanding of ferroelectric properties In thin films, surfaces are expected to affect the ferroelectric properties, thereby affecting the device performance. By contrast memory application demands polarization switching, which can be achieved by the movement of domain walls. Therefore, a detailed understanding of the domain wall structure, kinetics and dynamics is necessary to tailor the ferroelectric materials. Moreover the movement of these domain walls can be affected by the presence of grain boundaries and other defects in the system. This work focuses on providing an atomic level understanding of the surface effect in ferroelectric thin films, the structure and energetics for 180 domain wall and the effect of grain boundaries on the ferroelectric properties. 1. 1 What is Ferroelectricity? Ferroelectricity is the phe nomenon in which a material show s a spontaneous polarization that can be reversed by the application of a physically realizable external electric field [7, 8] This phenomenon is observed in materials which belong to the noncentrosymmetric point groups. In the point group classification, 21 out of 32 point groups do not possess c entrosymmetry. Crystals belonging to 20 of these point groups ( all except the 4 32 point group) can generate electricity in response to applied mechanical stress. This characteristic is known as PAGE 19 19 piezoelectricity. Out of these 20 point groups only 10 groups develop spontaneous polarization as the temperature of the crystal is changed. This response of the crystal is known as pyroelectricity Ferroelectrics are those pyroelectric cystals whose spontaneous polarization can be reversed by an external electric field. Therefore, ferroelectric materials also show piezoele ctric and pyroelectri c effects. A comparison between the ferroelectric and ferromagnetic material is inevitable. Ferroelectric materials possess a spontaneous electrical polarization that is switchable with an external electric field; where as ferromagnetic materials possess a spontaneous magnetization that is switchable with an external magnetic field. In the presence of a domain wall (separation of spin up with spin down states), the magnetization vector rotates either parallel to the domain wall (Bloch wall) or normal to the domain wall (Nel wall). The width of the magnetic domain walls are of the order of microns, as this is the length scale required for the spins to overcome the weak spontaneous magnetization that couples with the lattice strain. However, the ferroelec tric domain walls are traditionally represented as very sharp and as having narrow width. This is because the coupling between the ferroelectric polarization and lattice strain (spontaneous electrostriction) imposes a significantly larger energy cost for r otation of the polarization parallel or perpendicular to the domain wall. Thus, ferroelectric domain walls are predominantly considered as Ising type. 1.2 Effect of Electric Field in Ionic Crystals In ionic crystals (dielectric materials) the application of an electric field aligns the cations and anions based on the electrostatic interaction. This rearrangem ent results in electric dipoles, where p olarization of the crystal is calculated by the sum of these electric dipoles per unit volume. Figure 1.1 shows the three primary contributions to the electric dipole in ionic crystals. PAGE 20 20 Figure 1 1. Schematic of different sources of electric dipoles in an ionic crystal by the application of external electric field [9]. Figure 1 2. Schematic of frequency dispersion behavior of a typical ionic crystal with multiple rel axation and resonance modes [10] PAGE 21 21 The electronic, ioni c and dipole contributions to the overall polarization of the ionic crystal depend on the frequency of the applied electric field because the relaxation or resonance time for each of these responses is different. The electronic polarization can follow electric fields with THz PHz (1012 1015 cycles/second) frequencies, whereas the ionic polarization can respond to a GHz THz (109 1012 cycles/second) frequency electric field. Figure 1.2 illustrates the effect of the electric field on die lectric constant. T he resonance response for electronic, ionic and dipole reorientation (represented as space charge relaxation in Figure 1 2) with frequency is shown. 1. 3 Origin of Spontaneous Polarization At the macroscopic level, ferroelectricity can be understood in term s of a competition between electrostatics and elasticity. The origin of spontaneous polarization can be described by the total dipole and elastic energy contributions due to the displacement of atoms in the crystal. The dipole energy contribution for the displacement of one type of atom in the crystal is described as [9]: loc dipole dipoleE N w N W (1 1) Where, N is the total number of atoms per unit volume in the crystal, is the dipole moment of the cry s tal, and Eloc is the local electric field which is the driving force for this ion displacement. The and Eloc are defined as P E and Ploc 0 03 3 (1 2 ) Where, is the ionic polarizability of the displaced ions, is the Lorentz factor 0 is the permittivity of vacuu m (8.854 x 1012 F/m) and P is the spontaneous polarization For an isotropic cubic system, the Lorentz fac t or takes the value = 1 [11] Therefore, combining Equation 1 2 with Equation 1 1, the dipole energy can be written as [9]: PAGE 22 22 2 2 0 29 P N Wdipole (1 3) The displacement of atoms due to this dipole interaction also generates an elastic energy. Assuming a displacement of u the increase in the elastic energy of the crystal per unit volume can be written as [9]: 4 22 / 2 / u k u k N Welastic (1 4) Where, k and k are the force constants. The displacement u and the cha rge of the atoms ( q ) can be related to the spontaneous polarization as Nqu P (1 5) Therefore, the electrostatic energy can be written in terms of spontaneous polarization by combi ning Equation s 1 4 and 15. 4 4 3 2 24 2 P q N k P Nq k Welastic (1 6) Thus, the total energy is given by adding the dipole ( Equation 1 3) and elastic ( Equation 1 6) contributions. 4 4 3 2 2 2 2 0 24 2 9 P q N k P Nq k P N W W Welastic dipole total (1 7) Figure 1 3 shows the individual and total energy contribution towards the origin of fe rroelectricity. From Equation 1 7 and Figure 1 3 it is clear that if the coefficient of the harmonic term ( P2 terms) of the dipole dipole coupling is equal to or less than the elastic energy contributions, then the system will not show any spontaneous polarization. However there will be a spontaneous polarization when the dipole dipole harmonic term dominates the elastic energy harmonic term The Lorentz factor () and the ionic polarizability () play an important PAGE 23 23 role in stabilizing the spontaneous polarization. Since the Lorentz factor is much higher (102) for perovskite like crystal structures (like BaTiO3, PbTiO3) than that found in othe r crystal structures, perovskites are widely investigated for ferroelectric properties. This research focused on PbTiO3, which is also a perovskite. Figure 1 3. Schematic of the energy contribution due to A) dipole dipole coupling, B) elasti c energy, and C ) the total energy to explain the origin of ferroelectricity [9]. 1.3.1 Phonon Modes Atoms in a crystal can vibrate, and these vibrations are not random but rather a superposition of vibrations of the atoms and their neighbors. A combined effect of the vibrations of atoms in a crystal corresponds to waves of s pecific wavelength and amplitude. Ph onons are the quantized packets of energy generated by the vibration of atoms in a solid [12, 13] A B C PAGE 24 24 There are two types of modes of vibrations of atoms in a crystal: longitudinal and transverse. If the displacement of atoms from their equilibrium lattice position is along the direction of the wave propagation, it is known as the longitudinal mode. On the other hand, if the atoms move perpendicular to the direction of w ave propagation, it is termed a transverse mode. Along with th e vibration modes, the relative movement of atoms can generate different waves. Considering a diatomic system, an acoustic branch is developed when both the atoms move in phase, and an optical branch is developed when the atoms move out of phase. For a system with only one atom per unit cell, the phonon dispersion curves are represented by only the acoustic branch. In general of a system with N atoms per unit cell, there will be 3 acous tic branches (1 longitudinal, 2 transverse) and 3N 3 optica l branches (N 1 longitudinal, 2N 2 transverse) Phonon dispersion provides a relationship between t he phonon frequency and the wave vector. The concept of s oft phonon mode can be used to understand phase transitions in ferroelectric materials. Considering only one type of cation movement for the sake of discussion (only B type for an ABO3 structure), Figure 1.4 represents a schematic of some of the lattice vibrations, or eigenmodes that this lattice might display [9]. Figure 1.4(a) shows the initial cubic structure with all the atoms sitting at their respective centrosymmetric positions. Figure 1.4(b) shows only the elongation of the lattice along the z direction, while th e lattice positions are maintained in their respective symmetrical positions. Figure 1.4(c) structure exhibits coherently shifted B cations, while Figure 1.4(d) structure shows an antipolarized shift of the B cations. If any of these particular states lowe rs the crystal energy with respect to the undistorted lattice, the ions will adjust correspondingly to minimize the system energy, thereby stabilizing one of the corresponding distorted crystal structures. Assuming the lattice changes from Figure 1.4(a) to Figure 1.4(b), only oxygen octahedra are PAGE 25 25 z A B C D distorted without generating any dipole moment in the system. This effect can be observed as an acoustic mode in the phonon dispe rsion curve. T ransformation of Figure 1.4(a) to either Figures 1.4(c) or (d) generate s a dipole moment in the system due to the relative displacement of the cations and anions. These are observed as optical modes in the phonon dispersion curves. Figure 1.4(c) corresponds to a ferroelectric state, while Figure 1.4(d) represents an antiferroelectric state. Figure 1 4. Schematic representation of some of the possible eigen lattice vibration modes in a perovskite material (ABO3 structure), A) cubic structure, B ) symm etrically elongated structure, C ) coherent ly shifted central cation and D ) antipolarized shift of central cations [9]. Considering an ABO3 structure, the A cations are not shown for clarity and the B cations are shifted for the discussion. Therefore, if the high temperature cubic phase transforms to a ferroelectric state as the temperature decreases, the vibration mode frequency of the cub ic phase will decrease PAGE 26 26 B A (representing a soft phonon mode), and at the phase transition temperature have a negative frequency squared (i.e., imaginary frequency), corresponding to a vibrational mode with an amplitude that grows with time; that is the system is unstable [9]. In particular, a long wavelength (zone center) transve rse optical mode softens for a ferroelectric instability. It is important to note that this description is applicable only to displacive models [8]. Figure 1 5. Phonon dis persion curves calculated for A) cubic, and B ) tetragonal PbTiO3 using shell model (solid lines) and ab initio (solid circles) calculations [14] The vectors in the figure on the right represent the path of the Brillouin zone in the reciprocal space. PAGE 27 27 Calculations of phonon frequencies to describe the f erroelectric phase transition have been widely reported in the literature. The characteristic feature of this phenomenon in the phonon dispersion curves is represented with a negative phonon frequency (imaginary frequency), which indicates that the eigen s tate is unstable. Therefore, the phonon frequency calculations for a cubic (paraelectric) to tetragonal (ferroelectric) phase transition should show negative phonon frequencies for the cubic phase which is unstable. Figure 1.5 represents the phonon dispers ion curves for cubic PbTiO3 showing negative phonon frequencies, whereas there are no negative frequencies observed for tetragonal PbTiO3 [14] The instabilities at the and R points represent the ferroelectric (FE) and antiferrodistortive (AFD) modes respectively. Tinte et al. [15] have shown the FE and AFD instabilities in BaTiO3 and SrTiO3 respectively using atomic level simulation methods. Further discussion on AFD SrTiO3 is presented in Se ction 6.4.1. 1 .4 Basics of Ferroelectricity in Perovskites BaTiO3 (BT) was the first ceramic material in which ferroelectric behavior was observed [16] BT belongs to the perovskite family, having a gen eral formula of ABO3 (isostructural with the mineral perovskite, CaTiO3). The general structure of perovskite is based on a cubic closedpacked assembly of composition AO3, each A ion is 12 coordinated with oxygen ions and B ion occupy the oct ahedral inter sticies (Figure 1 6 ) [17] Figure 1 6 Schematic of the unit cell of BaTiO3 [17] = Ba = Ti = O PAGE 28 28 Ferroelectric materials show certain common properties. At high temperatures they have high dielectric constants. As the temperature decreases the dielectric constant increases until a transitio n temperature known as Curie temperature [18] At this temperature the crystal changes to a lower symmetry and a spontaneous polarization is developed. In BT, this is manifested as a cubic to tetragonal phase transition. Experimentally X ray and neutron diffraction studies are performed to observe the phase transition [19] In ferroelectric perovskites there may also be lower transition temperatures present, e.g., BT shows transition from tetragonal to the orthorhombic and rhombohedral phases [9, 17] as the temp erature is decreased (Figure 1 7 ). Figure 1 7 Phase transition and spontaneous polarization variation with temperature in BaTiO3 [9]. Due to this spontaneous polarization, an apparent surface charge densi ty is observed accompan ied by a depolarizing field. As the minimum energy state is the most favorable state for a given condition, the ferroelectric phase will try to minimize the energy associated with polarization by twinning. This divides the crystal in to many oppo sitely polarized domains. The PAGE 29 29 type of domains formed in the crystal depends on the crystal structure. For example, the cubic to tetragonal phase transition generates both 90 and 180 domains in BT. Figure 1 8 shows a schematic representation of the formation of 180 domains of altern ating up and down polarization resulting in an overall polarization of zero However, by applying an electric field parallel to one of the polar directions thereby switching the dipoles the crystal can be transfor med from a multidomain state into a single domain state [17] 9 0 domains walls in BT separate the domains, which not only differ in spontaneous polarization but also in spontaneous deformation [20] Figure 1 8 Schematic of A ) surface charge associated with spontaneous pola rization, and B) formation of 180 domains to minimize the electrostatic energy [17] Figure 1 9 shows a prototypical polarization electric field hysteresis loop Initially an unpoled crystal behaves as a normal dielectric when a small electric field is applied (segment O A). As the electric field increases, the domains change orientation, thereby increasing the polarization of the crystal (segment A B). The state of saturation is reached by increasing the field, where the crystal will consist of a single domain (segment B C). On reducing (polarization path C D) and reversing the electric field (polarization path D F G) the dipoles switch. The value OD is called the remnant polarization (Pr). At G the dipoles are pointing in the opposite direction. Increasing the electric field again causes a reversal of the dipole direction and A B PAGE 30 30 polarization follows the G H C path. T he value OF is called the coercive field ( Ec). Extrapolation of the linear portion of B C gives the value of spontaneous polarization at E = 0. Figu re 1 9 Hysteresis behavior of the polarization, in relation to the applied electric field for a ferroelectric crystal [21] 1.5 Theoretical Approach to Describe Ferroelectricity The Landau Ginsburg Devo nshire (LGD) formalism gives an excellent description of the ferroelectric behavior of many materials The symmetry invariant free energy terms for the paraelectric phase of PbTiO3 (m 3 m) can be written as a Taylor series expansion of the Cartesian components (in the <100> directions) of pola rization, Pi ( i =1 3), and of the stress Xi ( i =1 6 in Voigt notation) as [22] : 2 1 6 1 3 5 3 2 4 44 2 2 2 1 3 2 1 2 3 2 2 3 2 2 1 12 2 3 3 2 2 2 2 1 1 11 2 6 2 5 2 4 44 1 3 3 2 2 1 12 2 3 2 2 2 1 11 2 3 2 2 2 1 123 2 2 2 1 4 3 2 1 2 3 4 2 2 3 2 2 4 1 112 6 3 6 2 6 1 111 2 1 2 3 2 3 2 2 2 2 2 1 12 4 3 4 2 4 1 11 2 3 2 2 2 1 12 1 2 1 P P X P P X P P X Q P P X P P X P P X Q P X P X P X Q X X X s X X X X X X s X X X s P P P P P P P P P P P P P P P P P P P P P P P P P P P GP P P X X X X X X (1 8 ) PAGE 31 31 Here X i X ij and X ijk are the dielectric stiffness and higher order stiffness coefficients at constant stress; P ijs is the elastic complian ce coefficie nt at constant polarization; ijQ is the electrostriction tensor The second order terms in the free energy function are isotropic The fourth order term has its minima along the <100> axes T he sixth order term has its minima along the (<110> direction. Hence, lowering the temperature will favor a change in the direction of polari zation from [001] oriented to [011] oriented The parameters for the LGD description of PT are well known and this approach has been widely used to de scribe the bulk ferroelectric properties [22] ferroelectric s in the presence of surfaces [23] and ferroelectri c nanostructures such as ultrathin films [24] nanoparticles [25] and nanowires [26] 1.6 Appl ication s of Ferroelec tric Materials The effect of temperature on the crystal structure, and hence on the properties, leads to many applications of ferroelectric materials Figure 1 10 illustrates the properties based on different temperature ranges of a typical ferroelectric material. The spontaneous polarization decreases with increase in temperature and goes to zero at the Curie temperature (Tc). Therefore, the material is ferroelectric below Tc and paraelectric above Tc. The permittivity ( ) diverges n ear Tc, while the reciprocal of permittivity (1/ ) is known to be linear with respect to temperature in the paraelectric region (Curie Weiss law). Around Tc, ferroelectric materials can be used for high permittivity capacitors (Figure 1 10a) due to the hi gh dielectric constant at the transition temperature. At relatively lower temperatures, they can be used for memory applications (Figure 1 10b). The strong temperature dependence of the spontaneous polarization below Tc is utilized for detection and imagin g in pyroelectric sensors (Figure 1 10c) The piezoelectric effect in the ferroelectric region can be PAGE 32 32 used for piezoelectric transducers (Figure 1 10d), whereas the piezoelectric effect in the paraelectric phase c an be used for electrostrictors in electrom echanical devices (Figure 1 10e) The high relative permittivity in the paraelectric region can also be used for electro optic devices (Figure 1 10f). Ferroelectric materials are also used as positive temperature coef ficient (PTC) thermistors. O f all these different possibilities [3 5] the memory application directly uses the ferroelectric p roperty of the crystal. Figure 1 10. Schematic of the temperature dependence of spontaneous polarization, permittivity and inverse permittivity in a ferroelectric material. The application of the ferroelectric material at different temperature ranges are indicated with A) F ) [9]. 1.7 Motivation for This Work 1.7 1 Ferroelectric Memories The current Si technology used for memory applications is limited in its ability to produce fine scale capacitors. Ferroelectric thin films, with high permittivity, have been considered as a possible solution to this problem [9]. The memory technology can be divided into two types: A Capacitor B Memory C Pyro sensor D Piezoelectric transducer E Electrostrictor F Electrooptic device PAGE 33 33 A B volatile and nonvo latile [8, 9, 27] Dynamic Random Access Memory (DRAM) is widely used for memory application and is a volatile memory. The necessity of nonvolatile memory in mobile devices, such as mobile phones, digital cameras and MP3 players, ha s be en successfully achieved by FLASH based technology. However, FLASH memories suffer major limitation due to their slow write access time, overall high power consumption, scalability and poor endurance [28] Ferroelectric materials have been considered as a leading candidate for the next generation non volatile memories. The first Ferroelectric RAM (FeRAM ) was introduced in 1988 using PbZrxTi1xO3 (PZT) [27] The subsequent development and projection of high density Fe RAM are illustrated in Figure 1 11. Figure 1 11. A ) The development of FeRAM technology with respect to the design rule and memory density [27] and B ) t he memory capacity roadmap for a standard FeRAM [29, 30] Commercially ferroelectric memories are used in integrated memory devices involving several other memory technologies, or used for some niche applications [28] However, a standalone non volatile ferroelectric memory for commercial application demands a better understanding of the material system. PAGE 34 34 One major limitation that hinders the widespread use of ferroelectric devices is need for a large vol tage across the film during each reading cycle [9]. This is necessary t o switch the polarization. However, repetitive application of large voltage s develops regions which cannot be switched with the external field due to the pinning of the domain walls [8]. This phenomenon is known as polarization fatigue [8, 9, 28, 31] The domain pinning is reported to occur at grain boundaries, at surfaces or at other defect sites. The effect of surfaces and grain boundaries on polarization will be discussed in detail in Chapters 4 and 6 respectively. D epinning of the trapped domain walls recovers the original switched polarization values and thereby the hysteresis loops [8]. F atigue can be further divided into three different categories: ageing, imprint and dielectric relaxation [31] F atigue resistance is proposed to be improved by improving the film fabrication process, by improving the electrode materials or by searching for new ferroelectric materials [9]. 1.7 .2 Possible Future Applications of Ferroelectric Materials There are a few potential future applicatio ns of ferroelectric materials which are in the earlier stages of development. Those include ferroelectric nanostructures, MFSFET (Metal Ferroelectric Semiconductor Field Effect Transistor), ferroelectric device fabrication using atomic force microscopy (AF M), and ferroelectric cooling devices [9, 28] Ferroelectric nanostructures are interes ting as ferroelectricity is retained over a very small size (~ 20 nm). The critical thickness is reported to be of the order of a few unit cells below which ferroelectricity is destroyed in thin films [32 36] The studies on ferroelectric nanostructures will help develop highdensity arrays of ferroelectric nanocapacitors that will be used for high capacity memory devices. Synthesis of the se high density arrays in a cost effective and timely manner is a key issue that is being addressed. For MFSFET application, the gate dielectric in the conventional FET is replaced by a ferroelectric material. The major issue reported for such devices is t he poor retention time. AFM PAGE 35 35 is also predicted to be used to write extremely high density arrays of ferroelectric domains for memory applications. Ferroelectric cooling devices are proposed based on the current report of substantial cooling achieved with moderate voltages applied to ferroelectric thin films. 1.8 Objective of This Work All these current and future applications require a better understanding of ferroelectric thin films. Efforts to understand the physics of ferroelectric thin films started a round 1970 [8]. The study of thin films is complicated by phenomena which are absent in the corresponding bulk materi als. In thin films, the substrate material mainly guides the microstructure development of the film through substrate texture. Presence of domain walls and grain boundaries can affect the properties of the thin film. In order to enhance the atomic level u nderstanding of ferroelectric thin films, this study fo cused on three major objectives, (i) t o establish the coupling of surface terminations and polarization for (100) term inated ferroelectric thin films ( presented in Chapter 4 ); (ii) t o characterize the effect of domain walls on ferroelectricity. A detailed discussion of the structur e a nd energetics of 180 domain walls is reporte d in Chapter 5 ; and (iii) t o examine the influence of grain boundaries on polarization. This investigation focused on the effec t polarization in the presence of special boundaries for paraelectric and ferroelectric materials and is addre ssed in Chapter 6 Along with the above Chapters, this dissertation covers the fundamentals of perovskites (ABO3), in particular the structure a nd properties of PbTiO3 (Chapter 2) A brief description of the simulation methods (first principles calculations and atomistic simulation) used in this stud y is described is in Chapter 3. T he effect of parallelization on the computational power of the in house simulation code to predict complex ferroelectric structures is reported in Chapter 7. T he conclusion and future work are presented in Chapter 8 . PAGE 36 36 A B X CHAPTER 2 AN OVERVIEW OF PEROVSKITE STRUCTURED MATERIALS 2.1 Introduction to Perovskites An ideal perovs kite is defined with a general formula of ABX3 [37] The A site cations are typically larger than the B site cations and are similar in size to the X site anions. Figure 21 shows a schematic of the perovskite structure, which indicates that the A site cations are surrounded by twelve anions, and the B si te cations are surrounded by six anions in an octahedral coordination. The perovskites with ideal structure a have space group m Pm 3 (#221) Not all ABX3 systems that are perovskites. Figure 2 1. Schematic representation of an ideal ABX3 perovskite structure. The structure shows the cubooctahedral and octahedral coordination of the A site and B site cations respectively [37, 38] PAGE 37 37 2.2 Perovskite Oxides Perovskite is the general name for a family of materials with ABX3 formula, which are named after the mineral perovskite CaTiO3 [28, 37, 38] In particular, the perovskite oxides with composition ABO3 are a very large family. All of A sites can be occupied by a single cation or by more than one cation; likewise, the B sites can also be occupied by only one cation or by multiple cations [28] Depending on the composition and cation ordering, the physical properties of the entire family can vary widely. Some of the applications are summarized in Figure 2 2 [38] Figure 2 2. Summarized list of various applications of perovskite materials emphasizing different functionality [38] The chart does not incl ude ferroelectric, pyroelectric or optoelectronic applications of perovskite materials. PEROVSKITE OXIDES ABO3 Insulator Sr+2Ti+4O3 2 Piezoelectric Pb+2(Zr,Ti)+4O3 2 Catalyst La+2(Co Mn)+ 2 3O3n 2 High k Capacitor Ba+2Ti+4O3 2 Metallic Conductor La+3Cr+3O3 2 Major Constituent of Earth Mg+2Si+4O3 2 Superconductor Ba2 +2Y+3Cu3 +2 3O78 2 High Melting Point Ba3 +2Mg+2Ta2 +5O9 2 Giant Magnetoresi s tance L a+2Mn+2 4O3 2 PAGE 38 38 O II B A O I O III In order to represent the ABO3 structure with a five atom basis, the st ru cture in Figure 2 1 is redrawn as in Figure 2 3. In this representation, the A atoms sit at the corners of the cube, while the B atom sits at the body center and the O atoms are located at the face centers. There are three different types of oxygen in the structure named OI and OII (the oxygens on the (100) BO2 plane, also referred as O_B) and OIII (the oxygen on the (100) AO plane, also referred as O_A in this dissertation). The A atoms are surrounded by 12 O atoms. The B atom has six first nearest neighbor oxygen forming a regular octahedron. The oxygens have a lower symmetry environment with 4 A atoms and 2 B atoms neighbors. Figure 2 3. Schematic representation of ABO3 perovskite structure with A atoms at the corners, B atom at the body center an d O atoms at the face center. The highlighted atoms represent the five atoms necessary to describe the unit cell. Most of the perovskites are considered to be ionic compounds. Therefore, considering a hard sphere model with ionic radii RA, RB and RO for A, B and O atoms respectively, Goldschmidt (1926) observed that the A cations can fit the 12 coordination precisely if the ionic radii of A and O atoms are equal [28, 37] So, taking the face diagonal into account, (2RA + 2RO) = 2 (2RB + 2RO) (2 1) Or (RA + RO) = 2 (RB + RO) (2 2) PAGE 39 39 In real perovskite structures, however, this relationship the deviation from this ideal situation can be measured by defining a tolerance factor (t), where t = (RA + RO) / [ 2 (RB + RO)] (2 3) For ideal perovskites, t = 1. For t > 1, the A O distance is large and the B atoms are too small for the oxygen octahedral. Therefore, the structure can develop a small distortion, resulting in a polar system, such as BaTiO3 (t = 1.0706 [39] ). On the other hand, for t < 1, the A atom is small in comparison to the s pace available within the oxygen octahedral, thereby resulting in a lower coordination number for A atoms. If t is slightly less that one, then rotations and tilting of the oxygen octahedral is typically observed, as in SrTiO3 (t = 1.0091 [39] ). If t is much lower than unity, the overall structure will be distorted favoring a six coordination for the A atoms, like in LiNbO3, which is a perovskite like compound [28] 2.2 The Use of Perovskite Oxides The use of perovskite oxides started around the World War II [38] At that time mica (with a dielectric constant of 10) was the most important dielectric materials used as an insulator in most capacitors. The huge demand for mica and its limited availability opened up the avenue for finding a su bstitute material with high dielectric constant. Between 1940 46, various researchers identified BaTiO3 as a very promising high dielectric constant material [38] However, it was in the later part of the 1950s, when the extraordinary high dielectric constant was correlated with the ferroelectricity in BaTiO3. The major synthesis, characterization and modifications of the BaTiO3 structure were performed by two groups at Penn State [38] 2.2.1 Structure Field Maps Following the previous work, Roy and his students developed the Goldschmidt approach to provide details of a structure field map for different perovskites. This map shows the stability of different ABO3 phases with variation in the A and B site atoms (change in ionic radii). Even with PAGE 40 40 the hard sphere model and each ion interacting with pure electrostatic inte raction approximation, the results are observed to be extremely reliable [38] Figure 2 4 represents the lattice energy contour map for A+3B+3O3 2 compounds [40] The results include how the crystal depends on the tolerance factor. Figure 2 4. Lattice energy contour maps for predicting the sta bility of different compounds with A+3B+3O3 2 structures [40] 2.2.2 Ferroelectricity BT is the most widely studied perovskite oxide show ing ferroelectricity [28, 38] With decrease in temperature, the phases in BT change from cubic ( m Pm 3 symmetry (#221) ), to tetragonal ( P4mm (#99) ), to orthorh ombic ( A mm2 (#38)) to rhombohedral ( R3m (#160)) (Figure 1 7) Each of these transitions takes place with only small atomic dis placements, mostly by the Ti atms relative to the oxygen octahedral. The direction of polarization for the tetragonal, PAGE 41 41 orthorhombi c and rhombohedral phases are along <100>, <110> and <111> respectively. The same sequence of phase transition is observed for KNbO3, which is a +1 (K) and +5 (Nb) type ferroelectric, compared to +2 (Ba) and +4 (Ti) type for BaTiO3. 2.3 PbTiO3 In order to understand the complex solid solution behaviors, it is important to understand the end members of the solid solutions. PbTiO3 is the material studied in this dissertation. Even though PbTiO3 is structurally similar to BaTiO3, PbTiO3 show only a cubic to t etragonal ferroelectric phase transition below the Curie temperature. The degree of tetragonality is much higher in PbTiO3 with a c/a ratio of 1.08 compared to 1.01 in BaTiO3. In addition to the Ti atom, the displacement of the Pb atom is substantial and c ontributes positively towards the overall polarization. 2.3.1 Crystallography of PbTiO3 PbTiO3 (PT) transforms from a cubic (paraelectric) to tetragonal (ferroelectric) crystal structure as the temperature is cooled below its Curie temperature. The perovskite crystal structure for the cubic and tetragonal ph ase of PT is shown in Figure 2 3 and Figure 2 5 respectively Figure 2 5 Schematic of the tetragonal p erovskite crystal structure of PT. The lattice parameters are represented as a and c a a c PAGE 42 42 A B C a a a a a a a a 2 .3.2 Structural Details for Cubic PT Cubic PT belongs to the space group mPm 3 ; in complete notation of this is m m P 2 3 4 The number of this space gr oup is 221 and the point group sy mmetry is m 3 m Going thr ough each individual operation, t he first letter, the centering type, represents the translational symmetry in three dimensions. Here, P represents a p rimitive symmetry operation. The following t hree operations characterize the symmetry operations along th e primary <001>, s econdary <111> and tertiary <110> direction of the cubic system. For cubic PT, the 4/m represents four fold rotational symm etry at the corners and the cell center along with a mirror perpendicular to the four fold rot ational axis ([100]) (Figure 2 6 (a)). The 3 symmetry represents the roto inversion along the body diagonal of the unit cell (Figure 2 6 (b)). The 2/m symmetry operation represents the twofold rotation along the edge s of the unit cell along with a mirror p lane perpendicular direction to the rotation (Figure 2 6 (c)) Figure 2 6 Schematics showing three different symmet ry opera tions represented by cubic PT. A ) A four fold rotation symmetry along [100]. B ) A three fold inversion symmetry along [111]. C ) A t wo fold rotation symmetry along [110]. PAGE 43 43 2.3.3 Structural Details for Tetragonal PT Tetragonal PT belongs to the space group mm P 4 The number of this space gr oup is 99 and the point group symmetry is mm 4 Again P repres ents a p rimitive symmetry operation. + Figure 2 7 Schematic two dimensional representation of P4mm symmetry operation [41] The P4mm symmetry is obtained by adding the 4mm point group to a tetragonal primitive (tP) Bravais lattice (Figure 2 7) [41] The 2 D tP Bravais lattice is represented as a square unit cell (Figure 2 7(a)). To represent the 4mm point group, the four fold rotation a xis (shown as a solid diamond) along with four mirror lines (shown as solid lines) are included (Figure 2 7(b)). For any arbitrary lattice site at (x, y), the four fold rotation will replicate it at positions at ( y, x), (x, y) and (y, x). Similarly, th e mirror lines will result in ( x, y), (x, y), (y, x) and ( y, x) tP Bravais lattice 4mm P4mm A B C D PAGE 44 44 a c a c a B A lattice points [41] This 4mm point group symmetry is applied to each lattice site in the square unit cell. The mirrors are aligned along the cell edges and cell dia gonals (Figure 27(c)). This indicates that the square lattice must have a four fold rotational axis at the center of the cell. Therefore, the mirrors along the cell diagonal pass through the center four fold rotation. In addition, the four fold rotation a t the center develops new mirrors which pass through the center of the unit cell edges. This implies that the atom at the center of the unit cell also posses 4mm symmetry. Since the edges contain two intersecting mirrors, the edges must contain a twofold rotation axis at the edge centers (shown with a solid ellipse) (Figure 27(d)). The presence of parallel diagonal mirror lines through the center and the origin of the square lattice indicate the presence of another mirror line centered halfway between the two original mirror lines. These new diagonal mirror lines (shown as dotted lines) are in fact glide lines. Figure 2 8 Schematic representation of the break of h igh symmetry in tetragonal PT, A ) (100) view of atomic displacements restricting th e presence of a horizontal mirror plane, and B ) the atomic displacements along with the tetragonal structure hindering the three fold rotoinversion. Comparing the symmetries for cubic and tetragonal structures, the tetragona l structure does not contain a mirror perpendicular to the four fold rotational axis. This is due to the relative PAGE 45 45 cation and anoin displacements along the [001] directions for the tetragonal crystal structure to generate polarization in the system (Figure 2 8(a)). Also, the three fold r oto inversion is absent in the tetragonal crystal structure due to the change in the lattice parameter along with the atomic displacements (Figure 2 8(b)). PAGE 46 46 CHAPTER 3 SIMULATION METHODS Simulation is an ideal tool for probing ferroelec tric phenomena as it can provide information with atomic scale resolution, much of which is not easily accessible from experiment. For the present study, Density Functional Theory (DFT) and Molecular Dynamics (MD) simulation methods a re used to investigate the effect of surfaces, domain walls and grain boundaries on the ferroelectric properties of PT. Electronic structure calculations [42 44] are utilized to study the ferroelectric behavior. These calculation methods treat the system at the level of the nuclei and electr ons. DFT and Hartree Fock (HF) theory are most commonly used to study ferroelectrics. The quantum mechanical interactions between electrons are given by an effective exchangecorrelation potential in DFT (for further details on DFT, please check S ection 3. 2.1) Zero temperature phase diagrams, phonon frequencies and elastic constants can be computed using DFT. One of the first electronic structure calculation of BT and PbTiO3 (PT) was performed by Cohen and Krakauer [45, 46] The method is used to calculate the total energy of the system, the electr omechanical response, the effect of domain boundaries and surfaces, etc. [47] Electronic structure calculations are often computationally expensive. The calculations are limited to only a small syst em with most of the calculations performed at 0 K. In order to overcome the limitations of electronic structure calculations, atomistic simulation me thod (MD) are also used (further details are available i n Section 3.2.2). In this approach e lectrons are n ot treated explicitly for describing the internal structure of the atoms. Therefore, the method can handle a larger system size compared to electronic structure calculations (>106 vs. ~103 atoms respectively) [48] Atomistic methods are capable of incorporating finite temperature effect and the dynamic evolution of the system with time. PAGE 47 47 Indeed, atomistic simulation methods have been successfully used in ferroelectrics to study bulk properties [49, 50] solid solutions [51] superlattices [52] thin films [53] surfaces [36, 48, 54] and nanodots [55] This chapter provides a general overview of the simu lation methods used in this study. 3.1 Periodic Boundary Conditions Atomic level simulations are limited by system size (i.e., the total number of atoms within the simulation cell). In particular, DFT methods can simulate ~103 atoms, while MD simulations can handle ~106 atoms on a routine basis. O rigin of these limitations in the number of atoms will become clear in the Section 3.2. Periodic boundary conditions (PBC) allow the prediction bulk pr operties of materials, even when the total number of atoms in the system is small. Figure 3 1. Schematic of a 2D periodic system. The simulation is performed with the highlighted box, which is the actual simulation system. One component of the system is selected to show how PBC is applied on it (follow the arrow) [56] Actual System Replica Replica Replica Replica Rep lica Replica Replica Replica PAGE 48 48 Perio dic boundary conditions (PBC) replicate the system in one, two or three dimensions as required for the specific materials structure under consideration. When applied in three dimensions, this infinite extension allows the system to resemble a bulk material In order to describe the implementation of PBC, it is instructive to consider a 2D simulation cell. Figure 3 1 shows the 2D replication process of the actual simulation system, which is shown by the highlighted square. During simulation, the replica images always follow the actual simulation system. Therefore, if any atom in the actual system moves during the simulation, then its periodic image in each neighboring imaginary cell moves the same way. When an atom leaves the actual system during simulation, one of its images enters the system through the opposite side. Therefore, the total number of atoms in the actual simulation cell is always conserved. It is important to note that PBC may also lead to unphysical interactions (self interactions) between an atom in the actual systems and its image(s) created by PBC. Therefore, care must be taken to carry out the simulations with large enough system s to avoid such self interactions. As discussed below, the interactions between atoms are typically cut off at s ome finite range. In general, the actual simulation system size to avoid self interaction should be larger than twice the cut off. 3.2 Simulation Methods Used in This Study Sections 3.2.1 and 3.2.2 briefly describe the DFT and MD approach respectively 3 .2.1 First Principles Calculations The fundamentals of first principles calculations are based on the interactions of the electrons and nucleus in a system, which governs all the properties of materials. Quantum mechanics describes these interactions throu gh the Schrdinger Equation [57]. The time independent, nonrelativistic Schrdinger Equation can be written as: PAGE 49 49 E H (3 1) where H is the Hamiltonian operator, is the eigenstate wave function and E is the eigenvalue (energy). Equation 3 1 can be solved exactly for only a very few specific cases such as the hydrogen atom However, for more complex systems, no analytic solution is possible. In order to understa nd the individual contributions, with Born Oppenheimer approximation, the Hamiltonian in Equation (3 1) is expanded to: E r r U r Vm hN i i j j iN i i N i i 1 1 1 2 2, ) ( 2 (3 2) where, m is the mass of electron and N is the total number of electrons in the system [58] The individual terms in the Hamiltonian in Equation 3 2 are, in order, the kinetic energy of each electron, the electron nuclei interactions, and the electronic electron interactions respectively. The kinetic and electron nuclei contribut ions are typically relatively easy calculate. The electron electron interaction is a big challenge to solve, since for each electron the wave function associated with all the rest of the electrons need to be considered. The total wave function is a functi on of the spatial coordinates of the N electrons, or Nr r r ..., ,......... ,2 1 (3 3) This full wave function can be approximated to be the product of individual wave functions known as the Hartree product, and is expressed as: r r rN ........2 1 (3 4) However, the Hartree approximation is not antisymmetric, i.e., the swi tching two electrons do not change the sign of the wavefunction as is required for Fermions This issue is addressed in the Hartree Fock approximation. There are other h igh level quantum chemical calculation ap proaches developed, including configuration interaction (CI), coupled cluster (CC), Moller  PAGE 50 50 Plesset perturbation theory (MP), and the quadratic configuration interaction (QCI) approach [58] In spite of all the se modifications the wave function for a particular set of coordinates is still difficult to obtain directly. An alternate formalism can be developed based on the probability of finding an electron at a particular set of coordinate system. This el ectron density at a particular position r in space is given as: i i ir r r n*2 (3 5) where, the asterisk indicates the complex conjugate of the wave function. The prefactor of 2 in Equation 3 5 takes into account the Paulis exclusion pri nciple that each individual electron wave function can be occupied by two electrons with opposite spin. Considering N electrons in the system, this density approach reduces Equation 3 2 to be a function of 3 N coordinates. This electron density is the main motivation towards the development of density functional theory (DFT). 3.2.1 .1 Density Functional Theory (DFT) The connection between the wave function and electron density was made possible by two fundamental theorems. Hohenberg Kohn t heorem : The theorem proposed in 1964 by Hohenberg and Kohn [59] states that the groundstate energy from Schrdingers Equation is a unique functional of the electron density, n(r) This means that a one to one mapping is possible between the groundstate wave functions and the ground state electron density. This theorem also reduces the problem of solving a function of 3 N variables to finding the function of electron density. However, this theorem does not define the functional of the electron density. The property of the fu nctional is explicitly defined by Kohn Sham in 1965 [60] PAGE 51 51 KohnSham f ormulation : The KohnSham theorem [60] states that the electron density that minimizes the energy of the overall functional is the true electron density corresponding to the full solution of the Schrdinger s Equation Based on HohenbergKohn theorem, the energy functional can be written as [58] : i XC i known iE E E (3 6) Where, ion i i i i knownE r d r d r r r n r n e r d r n r V r d m h E 3 3 2 3 3 2 22 (3 7) and i XCE is the exchange correlation functional wh ich defines all other interactions not included by Equation 3 7. The terms in Equation 3 7 are, from left to right, kinetic energy of electrons, Coulomb interaction between the electrons and nuclei, Coulomb interactions between electrons, and Coulomb inter action between nuclei. In order to solve for i XCE Kohn Sham showed that the correct electron density can be expressed as solving a set of Equation s involving single electrons only. Therefore, Equation 3 2 can be written in terms of the K ohn Sham Equation as: r r r V r V r V m hi i i XC H ion 2 22 (3 8) where Vion (r) is the electron ion potential, VH (r) is the Hartree potential, VXC (r) is the exchange correlation potential, i is the KohnSham eigenvalue and i is the wave function of state i The Hartree potential is given by: 3 2) ( r d r r r n e r VH (3 9) The exchange correlation potential is given as: PAGE 52 52 r n r E r VXC XC (3 10) where, EXC(r) is the exchange correlation energy. Even though Equation 3 8 looks very much similar to Equation 3 2, there is one very distinct difference: Equation 3 8 does not contain any of the summations that are present in Equation 3 2. This makes replacement of th e electron wave function by the electron density much more useful from a computational stand point. For example, increasing the number of electrons in the system only increases the number of single particle Equations. Since the exchange and correlation ter m depends on the electron density, KohnSham Equation is solved self consistently. The solution of Equation 3 8 proceeds as follows [58] : a) an initial guess of the electron density, n(r) is made. b) the Kohn Sham Equation is solved to find t he single particle wave functions r) c) the electron density is calculated from the single particle wave functions obtained from step b) using i i i KSr r r n*2 d) The initial guess n(r) and the calculated nKS(r) are compared. If the two densities are the same then the groundst ate electron density has been achieved; if not the trial electron density is updated to nKS(r) and the steps b d repeat between until self consistency is reached. 3.2.1.2 Exchange Correlation Functional The key idea behind the KohnSham theory was to m ake the unknown contribution (EXC[n( r) ] ) as small as possible. Nevertheless, it is important to choose a good exchange correlation approximation for the system under consideration. PAGE 53 53 The most common approximation used is known as the local density approxim ation (LDA). In this approximation the exchange correlation potential at any position ( r ) is assumed to have the exchange correlation potential from the homogeneous electron gas having the same density. Since a homogeneous gas approximation is not always plausible in real systems, a more sophisticated description of the functional is given by the generalized gradient approximation (GGA) [61] The GGA functional includes the local electron density a nd the local gradient in the electron density. One would expect that GGA with a more elaborate description of the exchange co rrelation functional would yield much more in accord with experiment than LDA. However, this is not always the case [5 8] Hence, care must be taken in choosing the exchange correlation functional for individual calculations. 3.2.1.3 Plane wave Implementation Although Equation 3 8 simplified the Schrdinger Equation solving many electron systems such as solids is still computationally very demanding. The basic idea of the plane wave approach is to improve the calculation for the wave function (). Since most of the solid state systems are periodic, according to Blochs theorem [62] the wave function can be described as r iK j k je r U r. .. (3 11) where, Uj is a lattice periodic component, eiK.r is the plane wave like component and the plane wave vector K is unique only up to the first Brillouin zone in the reciprocal space. Reducing Uj to a sum of discrete plane wave basis sets with wav e vectors G which are reciprocal lattice vectors of the crystal, Equation 3 11 can be rewritten as: r G K i G G K j k je C r). ( .. (3 12) PAGE 54 54 where, r iG G G j je C r U. ,. (3 13) Cj,G is the plane wave coefficient and m l G2 whe re, l is the lattice vector and m is an integral. The dimension of the plane wave basis set should be infinite for an exact solution. Thus, truncating the basis set to a finite number of plane waves substantially, as required to make this a computational scheme, can provide accurate results. This truncation is achieved by defining an energy cut off given by cutE G K 22 The accuracy of the calculation can be increased by increasing the energy cut off. 3.2.1.4 K Point Sampling Plane wave like wav e functions reduces the infinite electronic wave function calculations to wave functions of finite particles at special reciprocal space known as K points [63] The total electronic potential is calculated by the occupied electronic states at these K points. Since the wave functions at closely spaced K points are nearly identical, a few K points are n ecessary to describe the whole reciprocal space. Therefore, the integrations can be replaced by the summations over a finite number of K points given by BZ j j jK F w K d K F 1 (3 14) where, F(K) is the Fourier transform of integral function such as electron density or total energy, is the supercell volume, and wj is th e weighting factor. I t is important to test the number of K points needed to achieve convergence. PAGE 55 55 3.2.1.5 Pseudopotentials The plane wave description for solving the KohnSham E quation requires a larger number of plane wave functions to describe the large curvature of the wavefunction, arising from the rapid changes in the electron density, particularly near the core of an atom. This requires a high enegy cutoff ( Ecut), thereby i ncreasing the computational cost. This issue is handled with the pseudopotential approximation. The idea of pseudopotential is to divide the atoms into core and valence electrons and solve for the valence electrons explicitly while maintaining the accuracy in the core region. Pseudopotentials replace the strongly oscillating wave functions in the all electron wave functions with smoother wave functions in the core region. Therefore, pseudopotentials reduce the calculating time for each iteration. The two most widely used pseudopotentials are ultrasoft pseudopotentials (USPP) developed by Vanderbilt [64] and the more recent projector augmented wave (PAW) method devel oped by Blchl [65] In genera l PAW methods are more accurate, because of a better description of the core electrons than USPP. The core radii in the PAW description is s maller than in the USPP, it require a larger energy cutoff com pared to USPPs. Kresse and Joubert compared USPPs and PAW method for small molecules and extended solids [66] The results indicate that well constructed USPPs and PAW method give almost identical results and are in good agreement with the all electron calculations. However, for systems with strong magnetic moments or materials with large differences in electronegativity, the PAW method yields more reliable results than the USPPs. Both the pseudopotentials are used for investigating the effect of interfaces on ferroelectricity in this study. PAGE 56 56 3.2.1.6 Charge Calculation The properties of materials can be fundamentally described by the charge transfer between interacting atoms. The charges on atoms in molecules or solids are not defined by quantum mechanical theory, which provides a continuous electronic charge density. Different schemes have been proposed to partition electrons between atoms or molecules in the system. Mulliken and Bader [67 69] are the two most commonly used schemes Mulliken charge analysis: Mulliken charge analysis is based on the electronic orbital method. This scheme can be used when basis functions centered on atoms are used for the calculation of electronic wavefunction of the system. Therefore, Mulliken charge on any atom is the charge associated with the basis functions centered on the particular atom. Even though thi s method is fast in calculating charges on atoms, it is sensitive to the basis set used in the calculations. The charge assignment becomes arbitrary for an infinite basis. For plane wave basis functions, which are not associated with any particular atom, M ulliken charge analysis is not applicable. Bader charge analysis: This scheme is proposed by Bader which is based on the electronic charge density. The charge space in any calculation is divided into different regions defined by the minima in the charge density. At the dividing surfaces there is no normal component of the gradient of the electron density. Since this scheme is based on the charge density, it is not affected by the basis set used for the calculations. It can also be used for plane wave base d calculations. This scheme is used in this dissertation for charge calculations. Detailed information on Bader charge analysis is available at http://theory.cm.utexas.edu/bader In this study, for the su rface simulations of PT thin films, the USPP method is used, while for domain wall simulations both USPP and PAW methods are used PAGE 57 57 3.2.2 Molecular Dynamics (MD) Simulation Molecular dynamics (MD) simulation is a deterministic approach, in which the motion of atoms in a system is predicted by solving Newton s Equation of motion. According to Newton s second law i i im a F (3 15) where, Fi is the force on an atom i ; mi and ai are the mass and acceleration of atom i respectively. The force on each atom can be expressed in terms of the gradient of the potential energy with respect to position. Vi i F (3 16) w here, V is the potential energy. This energy is calculated through the interaction of the atoms in the simulation system. 3.2.2.1 Interatomic Interactions For ionic systems, these interatomic interactions are described as a lon g range and short range contribution. A detail of the interactomic interactions is presented here. Long range i nteractions : The longrange interaction between charged ions in a syst em is represented by a Coulomb s um. It is given as: N i N i j ij j i Coulr qq E1 1 ,2 1 (3 17) where, N is the total number of ions in the system, qi and qj are the charges on atom i and j respectively, and rij is the distance between atom i and atom j. The r1 dependence makes the sum condit ionally convergent. Traditionally the Ewald method is used to evaluate the electrostatic energy and forces. Also, the electrostatic energy of an ionic crystal (Madelung energy) is short ranged. The Madelung energy of an ionic crystal corresponds to the ele ctrostatic energy for a PAGE 58 58 charge neutral system [70] Therefore, the Madelung energy can be compared directly to the Ewald sum calculations. We use the more recent direct summation method. Ewald s um : The Ewald sum efficiently calculates the Coulomb interaction between ions, its neighbors and its periodic images. The summation is conditionally convergent, i.e., the result depends on the order i n which the terms are added [56] In the Ewald sum, each point charge is treated as if it were surrounded by a charge distribution of equal magnitude and opposite sign. This charge spreads out radially from the point charge and its distribution is normally considered as Gaussian. Therefore, the total sum is divided to a short range (screening interaction between the neighboring charges) and a long range (charge distribution of the same sign as the original charge) component. The short range term is calculated in the real space, while the lon g range term is calculated in reciprocal space (Fourier space). Both of the terms converge quickly in their respective spaces. Care must be taken to subtract the self term, which includes the interaction of the cancelling distribution centered at ri with itself. A cut off distance is used to increase the simulation efficiency with the loss of accuracy. Typically a cut off of 2 or 3 times the lattice parameter is enough to obtain a well behaved long range energy convergence. However, there are some disadva ntages to Ewald sum. The method is not physically transparent, is complicated to code and scales O(N2) for N ions. The optimal balance between the real and reciprocal space contributions lead to a O(N3/2) at best [71] The direct summa tion method was developed to be a comp utatio nally more efficient approach in calculating the Coulomb sum. Direct s um : The direct summation was developed by Wolf et al [70] a nd is widely used to calculate the long range interactions. The key feature observation underyling the direct sum PAGE 59 59 was that the lack of electroneutrality makes the r1 sum diverge over spherical crystal lattice shells. In the direct summation method, a sim ple method was developed to achieve charge neutralization for a spherically truncated r1 pair potential. Figure 3 2 2D Schematic of the charge neutralization scheme used in the d irect sum method For any atom i interacting with atom j an eq ual and opposite charge is placed at the perimeter of the circular cut off to make the system charge neutral [70] This approach performs the Coulomb sum with a spherical cutoff (Rc) and place countercharges at Rc such that the system as a whole is charge neutral (Figure 3 2). The Coulomb sum is thus modified as c c i i N R r i j ij j i c CoulR R q q r q q R Ec ij (3 18) w here, qi(Rc) is the net charge within the cutoff Rc. This method shifts the potential and its derivatives smoothly. Damping can be used to substantially improve the convergence of the energy to the Madelung energies. +q j q q j q q +q q +q q q +q q +q +q +q i R c x y PAGE 60 60 The direct summation method has been s uccessfully used to predict bulk material properties, surfaces, grain boundaries, etc. Full details of the approach, along with the case studies showing the comparison with Ewald are given in reference [70] This approach is used for all the interatomic potential calculations performed in this study. Short range i nteractions : The l ongrange potentials described above are not sufficient to pr edict any materials properties because the r1 tends to infinity as r goes to zero. This means that the interaction between a +q and q charge will have the lowest energy at r = 0. This is known as the C oulomb catastrophe. If the material system only had C oulombic interactions, the ionic system w ould collapse. In order to prevent this, a short range interaction is used which is mostly repulsive in nature. The Buckingham potential is widely used to describe the short range repulsion in ionic systems. It is expressed as 6expij ij ij ij ij ij Buckr C r A r V (3 19) where rij is the separation between two ions i and j and Aij, ij and Cij are interaction parameters. Th e first term in Equation 3 19 is the repulsive interaction, which decreases exponentially with incr ease in distance between atoms i and j The second term in Equation 3 19 is an attractive contribution due to the dipole interactions. T here are other short range interatomic potentials available to describe ionic materials one of which is known as the R ydberg potential: ij ij ij ij ij ij Rydr r C A r V exp (3 20) where again Aij, ij and Cij are potential parameters. Both the Buckingham and Rydberg short range potential description s are used in this study. PAGE 61 61 3.2.2.2 Shell Model for Describing Ferroelectricity In addition to the long range and short range interactions, a shell model interaction is utilized to stabilize ferroelectricity. In the shell model [72] description, e ach ion is further divided into a core and a shell. The ch arge on each ion is described by the sum of the respective core and shell charges. The core and shell of each ion interact with the cores and shells of other ions via Coulombic interactions (longrange). The short range interactions are described by the sh ell shell interaction between the ions in the system. This is equivalent to the repulsion between electron clouds between atoms. It is important to note that the shell model mimics the presence of a core (nucleus + core electrons) and valence electrons (s hell). There is no explicit description of electron in MD simulations. The core and shell of each atom are coupled by a spring. In the traditional shell model, this coupling is described by a harmonic spring. However, the current study describes the core s hell interaction with both harmonic and anharmonic terms. Therefore, this core shell coupling is given by: 2 0 4 4 2 2) ( 24 1 2 1 B k k V for 0 (3 21) where is the coreshell displacement, and k2 and k4 are the harmonic and anharmoni c spring constants respectively. B is a core shell penalty term used to keep the core shell within physically acceptable separation [73] 3.2.2.3 Potential Descriptions Used in This S tudy There are two different potential parameter sets used in this study. These are denoted as Potential 1 and Potential 2 respectively for the rest of the discussion. The potential parameters for both the potentials are summarized in Table 3 1 PAGE 62 62 Potentia l 1, developed by Sepliarsky et al. [74] and Potential 2, developed Asthagiri et al. [73] differ in their description of the short range interaction and total atomic charges. In Potential 1, the short range interactions for the oxygenoxygen (O O) shells are described by the Bucki ngham potential ( Equation 3 19), where as the Pb Ti Pb O and Ti O short range interactions are described by the Rydberg potential ( Equation 3 20). For Potential 2, all of the short ranged interactions are described by a Rydberg potential ( Equation 3 20 ). Table 3 1 Potential parameters used for describin g interatomic interactions [73, 74] Atom Core charge (e) Shell charge (e) k 2 (eV/ 2) K 4 (eV/ 4) A (eV) () C (eV 6) Potential 1 Pb +4.958 2.785 119.48 17968.5 Pb Ti 0.096 2.420131 12.5665 Ti +8.820 5.158 1428.59 36411.0 Pb O 6766.270 0.273805 127.7793 O +0.563 2.508 23.29 5514.7 Ti O 1130.010 0.359723 160.8363 O O 3634.8 61 0.314424 331.6058 Potential 2 Pb +5.146 4 3.3506 75.35 26896.59 Pb O 6291.337 0.290791 3.6840 Ti +9.729 7 6.8449 1937.88 961.16 Ti O 1416.395 0.265259 296.2822 O +0.705 7 2.2659 26.48 1324.55 O O 283.410 0.520557 103.2676 B = 0 for Pote ntial 1 and 250 eV2 for Potential 2 The core shell penalty term ( B ) in Equation 3 21 is zero for Potential 1, while for Potential 2, a combination of B = 250 eV2 and 0 = 0.2 is used [50] This pen alty term ensures that the core shell displacements remain within an acceptable range (<~ 0 ) during the simulation. Much larger values of B have been used for Potential 2 in a previous study [73] The cut off distances used for Potential 1 and Potential 2 are 6.5 and 10.0 respectively. 3.2.2.4 Integration Algorithms The evolution of the atoms in MD simulations is achieved by integrating Newtons Equation of moti on. This integration is performed at a small increment of time ( t ) and the PAGE 63 63 position, velocity and accelerations of the atoms are calculated for the system evolution. There are various integration schemes available [75, 76] Each method has its advantages and disadvantages based on the accuracy and computational time. The Verlet algorithm [75] is one of the most simple integration algorithms. A 5thorder Gear pr edictor corrector method [56] is use d as an integrating scheme for this study. PredictorCorrector a lgorithm : The predictor corrector algorithm is one of the most commonly used algorithms for MD simulations. The main concept of the predictor corrector is that the positions, velocities, acc elerations, and higher order derivatives of position at time t + t is estimated by a Taylor expansion about t and then corrected based on the predicted acceleration. For a fifth order predictor corrector, the p redictor forms are expressed as : 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 (3 22) where r v and a are the position, velocity, acceleration of each atom respectively. b c and d are third fourth and fifth derivative s of the position with time of each atom respectively T he superscript p corresponds to the predicted values. From the predicted positions, the total potential energy of the system is calculated. From this potential energy, the interatomic forces at time t + t are calculated using Equation 3 16. From these forces, the accelerations are calculated. Since, the predicted acceleration values are PAGE 64 64 not based on physics these are corrected with respect to the correct accelerations ( t t ac ). Therefor e, the error from the predictor step is given as: t t a t t a t t ap c (3 23) Using the predicted results and the error the corrected positions and other derivatives are calculated by: 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 (3 24) 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 [56, 76] These corrected values are then used to predict the positions and other higher order derivatives for the next iteration This procedure is repeated to predict the trajectory of each atom in the simulation The size of the t ) is one of the most important factor for sampling the simulation trajectory. The average phonon period of a ty pical atomic vibration is ~200 fs. A time step of 1 fs is generally sufficient to capture the atomic vibration. However, for core shell model s the internal vibrations of cores and shells are much faster than the vibration of the overall atom. Therefore, f or all the calculations performed here, a sampling time of 0.02 fs is used to sample the trajectory of the cores and shells accurately. 3.2.2 .5 Optimization Methods There are several methods available for geometry optimization. PAGE 65 65 Figure 3 3 Illustation of local minima and maxima of a function in an interval [a,b] [77] The basic idea behind geometry optimization is to find a minimum of an arbitrary function (such as energy) with many variables by manipulating the geometric degrees of freedom. Considering a function represented in Figure 3 3 in an interval [a,b], points B and D represent local maxima, whereas points C and E are local minima. For geometry optimization, the task is to find the global minima in the energy space for the entire structure. There are various methods available to find the minimum (or root), such as the bise ction, golden search, steepest descent, conjugate gradient and Newton Raphson method. The Bisection approach is the simplest approach, in which the value of the function is calculated at the mid point of the interval. The interval is then updated and used as the new interval for the function evaluation. This process is repeated until the function stops changing. PAGE 66 66 The Golden Search method is analogous to the bisection approach. The basic idea in this approach is to bracket a minimum and the shrink the brack et. Figure 3 4 represent a schematic of the golden search approach, where the original function is bracketed by points [1, 3, 2]. During the process the function is evaluated at 4, which replaces 2, then the function is evaluated at 5 and replaces 1. In th e next step 6 replaces 4. Therefore, after three steps the minima of the function is now bracketed by [5, 3, 6]. This process can be repeated to achieve the minima of the function. Care must be taken to make sure the point at the center is the lower than t he two outside points. Figure 3 4 Schematic showing the golden search approach for a function brackted originally by [1, 3, 2] [77] Steepest descent, unlike bisection and golden search methods, uses the first derivative of the function to achieve optimization of any function. It is well known that the first derivative of a function at a minima or maxima is always zero. In order to find the minima for a one dimensional curve, the first derivatives are used to guide towards the minima. The general rules are: move to the right if the first derivative is less than zero; move t o the left if the first derivative is great than zero; and stop if the first derivative is zero. In higher dimensions a negative gradient PAGE 67 67 of the function is utilized to decide the path towards the minimum. This is known as steepest descent, where the algor ithm continuously calculates one dimensional minimization of the function along the direction of the steepest descent. Figure 3 5 Schematic of steepest descent method showing consecutive steps perpendicular to each other [77] Figure 3 5 shows the graphical representation of a steepest descent approach. The optimizatio n scheme moves in a steepest direction to a minimum along that direction. This approach converges slowly to the minimum. Further improvement to this resulted in the conjugate gradient method, in which the search direction is given by a combination of force and previous search direction. Therefore, the direction of motion is determined in succession at each step of the iteration. The number of iterations to achieve optimization compared to steepest descent is less, where as the time per iteration is much lar ger than steepest descent. PAGE 68 68 Newton raphson method is a bit more sophisticated. This method uses a gradient of the function to identify a direction and uses the curvature (second derivative) to predict where the function passes through a minimum. This approach is very useful when the simulated structure is close to the minimum. Due to these typical qualities of the optimization schemes, a combination of different optimization methods is usually a very useful approach. Simulated annealing is also another opti mization scheme, where the system is heated to release all the metastable phases and then cooled down slowly. If necessary, this process is repeated until optimization is achieved. In this research steepest descent, conjugategradient and Newton Raphson me thods are used at some point to optimize various PT and SrTiO3 structures. 3.2.2 .6 Ensemble Details The ensemble in which a simulation is performed describes the specific thermodynamic conditions. The state variables used to describe the thermodynamic cond ition of the system are pressure (P) and volume (V) temperature (T) and the internal energy (E) chemical potential ( ) and total number of atoms (N) In general, the ensembles are categorized as three types microcanonical, canonical and grand canonical ensembles [56] Systems with fixed N V and E during simulation are described as isochoric microcanonical ensemble (NVE). Similarly NPE is also the isobaric microcanonical ensemble. However, these ensembles do not correspond to most of the experimental conditions. Therefore, other ensembles are more useful for experimental comparison. The canonical ensemble is defined by fixing the N and T of the system. The temperature of the system is maintained at an average value throughout the simulation and the total energy (E) is allowed to fluctuate. The ensemble with fixed N, V and T is known as the isochoric canonical ensemble. The ensemble with fixed N, P and T is known as the isobaric canonical ensemble. In PAGE 69 69 NPT ensemble the pressure and temperature are maintained to a specified average value and t he volume of the system is allowed to change. For grand canonical ensembles, the total number of atoms (N) in the system is allowed to change so as to maintain a fixed uniform chemical potential, There are a few different grand canonical ensembles possible. For t he VT grand canonical ensemble the of the different species is maintained at a specified value along with the volume and temperature. However, the system allows the exchange of parti cles (N) with a surrounding system. A summary of the thermodynamic variables are mentioned in Table 3 2. Table 3 2 Summary of the eight thermodynamical ensembles and the ir corresponding fixed and dependent variables. P ressure (P), temperature (T ) and che mical potential for all species ( ) are the i ntensive variables ; and volume (V), total number of atoms (N), and the internal energy (E) are the e xtensive variables For specific conditions the internal energy E is modified as the enthalpy (H = E + PV), the Hill energy (L = E iNi), and the Ray enthalpy (R = E +PV iNi). Fixed Variables Dependent Variables Ensemble N V E P T Microcanonical N V T P E Canonical N P H V T Isoenthalpic isobaric N P T V H Isothermal isobaric V L N P T Grand microcanonical V T N P L Grand c anonical P R N V T Grand isothermal isobaric P T N V R Generalized Since NPT and NVT ensembles are the most widely applied, and will be used in this study, it is necessary to briefly discuss the constant pressure and constant temperature algorithms used during MD simulation. PAGE 70 70 Simulation System 3.2.2.7 Pressure Control There are various algorithms developed to maintain an average pressure in the system. The Andersen [56, 78] scheme uses a hydrostatic approach to maintain the pressure. Figure 3 3 shows the schematic of the Andersen scheme for a 2D simulation system. A set of pistons are attached to the system with a certain mass. It can be imagined that if the inst antaneous pressure in the simulation system increases than the targeted average pressure (which is on the pistons), then the system expands and vice versa. Figu re 3 6 Schematic of a 2D simulation system showing the piston setup to maintain pressure. Considering a 3D system with volume V and a piston of mass Q, the potential and kinetic energies of the piston are described as: 22 1 .) ( .) ( dt dV Q E K Energy Kinetic PV E P Energy Potential (3 25) Where, P is the pressure of the system. Equation 3 25 can be reduced to 2 3 / 2 3 / 12 1 . dt s d m V E K s V U E Pi i i (3 26) PAGE 71 71 Where, is is the reduced coordinate of atom i at time t and is represented as t s V t xi i 3 / 1 (3 27) Using Lagrangian formalism of classical mechanics, Equation 3 26 can be solved to obtain the following Equation s of motion: Q P p dt V d dt dV dt s d V mV t F dt s di i i 2 2 3 / 1 2 23 2 (3 2 8 ) Here Fi is the force on atom I, and p is the instantaneous pressure, which is given as i i i i i iF r V v m p 3 1 3 12 (3 2 9 ) Therefore, the Andersen scheme deals with hydrostatic pres sure, i.e., the pressure is isotropic. T his method allows the simulation system to expand or contract the same amount in each direction. For a more complex scheme, which can allow different expansions or contractions in different directions, Parrinello and Rahman [79, 80] extended the Andersen scheme to allow changes in both shape and size of the simulation system. This approach is used for this work. 3.2.2.8 Temperature Contro l In order to maintain the temperature during simulation, a thermostat algorithm wa s implemented. Similar to the constant pressure scheme, instantaneous temperature of the system is calculated during simulation. Based on Boltzmanns relationship, the macro scopic temperature (T ) of the system can be calculated from the average kinetic energy as: For the system having, the average internal kinetic energy < K > of the system is related to its macroscopic temperature through PAGE 72 72 T N k mv Kdf B N i i2 1 2 11 2 (3 30) where N is the total number of atoms in the system, m is the mass of each atom, v is the velocity of atom kB is the Boltzmann s constant, and Ndf is the number of internal degree s of freedom of the syste m. Therefore, the average instantaneous temperature < Tins> at any time is expressed as N i i df B insmv N k T1 22 1 2 (3 31) Th is temperature < Tins> describes the temperature of the simulation system. The instantaneous temperature is compared with the target temperature ( T0) to adjust the temperature of the system. There are various constant temperature schemes available, including velocity rescaling, Berendsen [81 ], Langevin [82 84] and Nos Hoover [85, 86] The simplest approach to maintain the temperature is simple velocity rescaling. Since temperature of the system is related to the velocity of each atom pres ent in the system, the velocities of the atom s can be adjusted in order to control the system temperature This is given by ins old newT T v v0 (3 32) where vnew is the rescaled velocity and vold is the velocity prior to rescaling. Most of the MD simulation results discussed in this report are performed with the velocity rescaling to maintain the temperature of the system. 3.3 Model Validation (Bulk PbTiO3) Before performing surface, domain and grain boundary simulations, it is important to show the reliability of the interatomic potentials in previously determined bulk properties. There are two different shell model potentials used for this study, namely Potential 1 and Potential 2 PAGE 73 73 (p arameters are mentioned in Table 3 1). These simulations are performed on a 5 x 5 x 5 unit cell supercell size (i.e., for 625 ions) with periodic boundary condition. The simulations are performed by a fifth order predictor corrector algorithm to integrate Newtons Equation of motion in an NPT ensemble with a 0.02 fs time step. The temperature is maintained in the system using velocity rescaling [56] The effect of temperature on lattice parameter in bulk PT is shown in Figure 3 7 As reported previously in the literature [50, 74] both interatomic interactions predict the tetragonal to cubic phase transformation with increase in temperature. The change in lattice parameter and c/a ratio with temperature shows qualitative agreement with experimental results [87 89] Figure 3 7 Lattice parameter as a function of temperature for PT. The experimental lattice parameter values (dashed lines) are from Shirane et al. [87] The lattice parameter s for Potential 1 an d Potential 2 are shown with filled symbols with open symbols respectively. Inset is the variation of c/a ratio compared with the experimental data. 3.80 3.85 3.90 3.95 4.00 4.05 4.10 4.15 4.20 0 200 400 600 800 1000 1200 Temperature (K) Lattice Parameter () 1.00 1.02 1.04 1.06 1.08 0 250 500 750 1000 1250 Temperature (K) c/a ratio Experiment Potential 1 Potential 2 PAGE 74 74 In addition, the ferroelectric to paraelectric phase transformation is confirmed with the evaluation of po larization within the bulk PT structure for both potentials. Figure 3 8 shows the comparison of polarization predicted by these potentials with experimental values [89 91] The variation of polarization with temperature is in very good qualitative agreement with experiment. The polariz ation value decrease s with increasing temperature and goes to zero, indicating the ferroelectric to paraelectric phase transition in PT. The transition temperature is predicted to be ~500 K and ~900 K for Potential 1 and Potential 2 respectively. Experimen tal ly the Curie temperature is observed around 763 K [18, 87] Potential 2 has also been reported to reproduce the elastic constants to within 2030 GPa (within 20 %) of the experimental values for tetragonal PT at 300 K [92] A quantitative comparison of both the potentials with experimental lattice parameters and polarization for PT is presented in Table 3.3. Figure 3 8 Polarization as a function of temperature fo r PT obtained with both the potentials. Experimental results are obtained from ref [91] 0 10 20 30 40 50 60 70 80 0 200 400 600 800 1000 1200 Temperature (K) Polarization (C/cm2) Potential 1 Potential 2 Experiment PAGE 75 75 Table 3 3 Quantitative analysis of bulk PT predicted from both interatomic potentials. Experiment Potential 1 Potential 2 Curie Temperature (Tc) 763 K 500 K 900 K Lattice parameter (a) at room temperature (~ 300 K) 3.904 3.886 3.867 Lattice paramete r (c) at room temperature (~ 300 K) 4.152 3.974 4.116 c/a ratio at room temperature (~ 300 K) 1.063 1.023 1.064 Lattice parameter at Tc 3.970 3.909 3.943 Polarization at room temperature (~ 300 K) 74.80 35.16 57.91 Comparison of both the p otential results c learly indicates Potential 2 displays much better quantitative agreement with experiment than Potential 1. This is expected since Potential 2 was developed after Potential 1 and included additional DFT derived data to fit the tetragonal P T bulk structure. Nevertheless, the degree of agreement between Potential 2 and experiment for bulk PT is somewhat fortuitous The DFT derived data for fitting tetragonal PT was performed with local density approximation. It is well known that the local de nsity approximation ( LDA ) underestimates the experimental equilibrium volume and therefore potentials derived from the DFT data should usually underestimate the transition temperatures [50, 54] However, Potential 2 was fit ted to a DFT database containing both PT and Pb(Mg1/3Nb2/3)O3 (PMN) information and was designed to simulate the compositional phase diagram of Pb(Mg1/3Nb2/3)O3xPbTiO3 [50] A s a side effec t of the inclusion of the PMN data wa s that the quality of the fit to the DFT data of PT wa s actually poorer, but in the case of Potential 2 this error wa s in t he direction of the experimental values, resulting in a more accurate potential than would be expected. Since the degree of fid elity of the description of bulk PT is different for the two potentials, a comparison of results obtained with each potential will allow the identification of generic system behavior that does not depend sensitively on the potential. PAGE 76 76 In addition to the empirical potentials, firstprinciples calculations based on DFT calculations were perform ed with V ienna abinitio S imulation P ack age (VASP) [93 95] The calculations were performed with Vanderbilt ultrasoft pseudopotentials (USPP) [64] at the level of LDA. The pseudopotential treats Pb 5d a nd Ti 3p states explicitly as valence states. All the calculations were carried out with a cutoff energy of 395 eV (29 Ry). A 6 x 6 x 6 k point Monkhorst Pack [63] mesh wa s used. The accuracy of the pseudopotentials and the effect of cutoff energy on convergence for PT has already been established in the literature [96 98] The bulk optimized lattice parameter for PT was predicted to be a = 3.86 and c/a = 1.05 (experimental values are a = 3.904 and c/a = 1.063 [19] ). Overall, atomic level simulation methods are e xtremely useful in studying ferroelectric materials. The DFT method yields high material fidelity but is limited to system size. Atomistic simulations, on the other hand, are capable of simulating large systems Both these simulation methods are wisely use d to address the surface, domain wall and grain boundary effects on ferroelectricity. PAGE 77 77 CHAPTER 4 EFFECT OF SURFACES O N FERROELECTRICITY 4.1 Introduction Ultrathin ferroelectric (FE) films have attracted considerable attention due to their numerous potential applications [3 5] including nonvolatile memory components [6] These devices demand the establishment of the thickness property relationship, and a complete understanding of the surface termination and domain effects. This chapter focuses on the structure and en ergetics of the (001) surfaces in PT using atomic level simulations. Previous simulations on ferroelectric thin films of PT [36, 54] have focused on the trends of the ferroelectric properties with system size. By contrast, this study primarily describes the effect of possible individual (001) surface terminations on ferroelec tricity in PT. Table 4 1. List of simulation details for PT thin films Simulation method Potential description Method addressed in the text as System size (unit cells along X, Y and Z) Other details DFT Ultra soft pseudo potential (USPP) with LDA approximation DFT 1 x 1 x 7.5 1 x 1 x 8 Vacuum size ~ 12 Constant volume Optimization Atomic level Sepliarsky potential Potential 1 6 x 6 x 32.5 6 x 6 x 32 Vacuum size ~ 40 Constant volume Optimization Atomic level Asthagiri potential Potential 2 6 x 6 x 32.5 6 x 6 x 32 Vacuum size ~ 40 Constant volume Optimization Symmetric slab Asymmetric slab Table 4 1 shows the simulation details used for this study. The predicted surface properties using two different empirical interatomic potenti als are compared to equivalent DFT results. PAGE 78 78 4.2 (001) Surface Terminations in PbTiO3 Free standing films with (001) surface terminations in PT can form two distinct (001) surface terminations, a PbO layer or a TiO2 laye r. Based on the nominal valence mode l of Pb2+, Ti4+ and O2 , both of these terminations are charge neutral In contrast, the (110) and (111) surface terminations yield non charge neutral surface terminations, and hence not included in this study. Combinations of all possible surface terminations result in three different (001) terminated thin films, PbO/ PbO PbO/ TiO2 and TiO2/ TiO2. The films with different terminations on the two sides (PbO/ TiO2) are known as asymmetric slabs, while films with same terminations (PbO/ PbO and TiO2/ TiO2) are know n as symmetric slabs. The asymmetric structures result in stoichiometric films while the symmetric structures are not stoichiometric. For thick enough films, the surface structures and energies predicted for each termination (PbO or TiO2) should be indepe ndent of the type of surface on the other side of the film i.e., symmetric or asymmetric slabs. The direction of polarization in ferroelectric thin films results in additional degree of freedom for (001) PT surfaces. For a thin film system, when the fer roelectric polarization is perpendicular to the surface, one of the terminations has polarization pointing into the film, whi le the polarization in the other termination points out of the film. P olarization pointing into the surface is termed In polarization, while polarization pointing out of the surface is termed Out polarization. The additional complexity of the direction of polarization with different (001) surface terminations results in four possible configurations for free sta nding films, as illustra ted in F igure 4 1 Simulations are performed to establish the effect of each individual surface terminations ( PbO In, PbO Out, TiO2In and TiO2Out ) on the ferroelectric properties of PT. PAGE 79 79 Figure 4 1 Schematic of all four possibl e free standing ferroelectric (001) P T films along with an atomistic view of an asymmetric slab In order to describe the surface with simulation within a three dimensionally periodic simulation cell, a block of atoms is placed in the simulation box, vacuu m is added in the through film direction ( z for our schematics in Figure 4 1 ) while the block fits perfectly in the x y planse such that periodic boundary condition in the x and y directions result in a defect free perfect crystal in that x y plane. Als o, to represent the effects of the presence of an infinitely thick film, the in plane ( x y plane) lattice parameters are fixed to the predicted bulk crystal values in the ferroelectric phase for each potential (a = 3.866 for Potential 1, a = 3.843 for Potential 2 and a = 3.860 for DFT). All the surface simulations are performed at 0 K by optimizing all the ions x y z PbO termination TiO 2 termination vacuum vacuum PAGE 80 80 at constant volume. The other simulation details specific to interatomic potentials or DFT calculations are given at the appropriate points be low. All the surface simulations with interatomic potentials are performed with a 6 x 6 unit cells in the x y (001) plane and 32 or 32.5 unit cells thick (in the z direction) for symmetric and asymmetric films. The vacuum size used for these surface simula tions are kept ~ 40 which is much larger than the cutoff of the interatomic potential thereby ensuring that the two surfaces do not interact through the vacuum.. As we shall see, these films are sufficiently thick for the middle layers of the film to ha ve the properties of bulk PT. Optimization is achieved by relaxing all the ions (both cores and shells) in the system until all of the forces are less than 0.0001 eV/. For DFT calculations, a 1 x 1 x 7.5 unit cell system is considered for all the s ymmetr ic slabs and a 1 x 1 x 8 unit cell is considered for the asymmetric slab calculations. These result in 15 and 16 cation oxygen layers along the z direction for the symmetric and asymmetric slabs respectively. The surfaces are created by adding a ~ 12 th i ck vacuum in the z direction A 6 x 6 x 1 k point Monkhorst Pack [63] mesh wa s used. In the absence of experimental results, the DFT results are used as benchmarks against which the interatomic potentials are compare d for the surface results Characterization of surfaces typically focuses on atomic relaxation, interlayer spacings, surface rumpling (s hown schematically in Figure 4 2 ) and surface energies. The displacement of a metal and oxygen ion along the z direction relative to the ferroelectric reference structure are given as z(M) and z(O) respectively. The interlayer distance (dij) is then defi ned as the difference between the cation displacements z(M) between layer i and j In particular, d12 is the cation spacing between the first (i.e., outermost) layer and the second layer, while d23 is the spacing between t he second and third layers (F ig ur e 4 2 (b) ). For these calculations, the cations PAGE 81 81 are taken as the reference because they scatter electrons more strongly than the oxygen ions, making them easier to detect experimentally [99] By convention, a negative value corresponds to an inward displacement, while a positive value corresponds to an outward displacement. Surface r umpling is defined as the amplitude of relative displacements of the cations and oxygen in the same crystallographic plane (F ig ure 4 2 (c) )., s = z(M) z(O) [100] For TiO2 planes, an average displacement of the two oxygen atoms is considered for calculation. Figure 4 2 <010> view of A ) P bO terminated PT with vacuum, B ) change in the interlayer distance with the first two interlayers illustrated in the figur e, and C ) s urface rumpling obtained after relaxation due to the movement of the individual atoms from the centro symmet ric position. The vacuum is not shown in B) and C) for clarity. 4. 3 Critical Issues Before examining the predicted surface energetics and structure, it is important to assess the effects of the limitations of the atomistic methods. A B C Vacuum d23 Interlayer Distance d1 2 s Surface Rumpling = Pb = Ti = O PAGE 82 82 4.3 .1 Use of Fixed Charge Simulations for Surfaces The most significa nt approximation of the empirical potentials used in this study, is that the charges on the ions are always fixed. Thus, there is no mechanism for the surface charges to change from those of the bulk. Figure 4 3 Comparison of the charge v ariation in each ion in surface and bu lk with Bader charge analysis A ) PbO te rminated, and B ) TiO2terminated cubic PT. The maximum deviation (~ 5 7%) is observed to be for the Pbions on the PbO and TiO2terminated surfaces. Similar result on charge variation is reported by Piskunov et al. ([43]) with Mulliken charge analysis. Cell 1 and Cell 2 represents the surface unit cell and the next sub su rface unit cell in the thin film respectively. B 1.5 1.0 0.5 0.0 0.5 1.0 1.5 2.0 2.5 Ti O_Ti Pb O_PbCharge (e) Cell 1 Cell 2 Bulk A 1.5 1.0 0.5 0.0 0.5 1.0 1.5 2.0 2.5 Pb O_Pb Ti O_TiCharge (e) Cell 1 Cell 2 Bulk PAGE 83 83 To assess the impact of this approximation, a Bader charge analysis [67] have been performed for the cubic PT surface determined from the DFT calculations. As shown in Figure 4 3 the surface charges are only slightly different from the bulk values. Cells 1 and 2 are the surface unit cell and the next sub surface unit cell in the thin film for each termination respectively. A maximum deviation of ~5 7% from the bulk charge value is observed for Pb ions on the PbO and TiO2terminated surfaces. The charge deviations for all the other atoms are less than 3%. Similar analysis performed with Mulliken charges [43] on cubic PT surface reported a charge variation of less than 10% for all the atoms at or near the surface. These small variations in the charge values of atoms at the surface compared with the bulk provided a posteriori justification for using fixed charge potentials for PT surface simulations. 4.3 .2 Short circuit vs. Open circuit Electrical B oundary Condition Another important issue to address for simulating ferroelectric thin films is the electrical boundary conditions. Thin film can be simulated under short circuit (film electrical boundary condition is complete with electrodes) or open circuit (free standing film with no electrodes to complete the electrical circuit) electrical boundary conditions. The short circuit electrical boundary condition is achieved by applying electrodes on both side of the thin film. The electrodes compensa te the charges on the surface of the thin film, thereby screening the depolarizing field. However, for opencircuit electrical boundary condition there are no electrodes present along with the thin film. The free standing film by itself does not complete the electrical circuit and hence exhibits an open circuit electrical boundary condition. There have been a number of DFT studies considering short circuit boundary conditions [101103] in symmetric slabs. This boundary condition can be achieved by (i) simulating the effects of metal electrodes, or (ii) by applying an external electric field in order to cancel the field inside the slab. The idea behind both the approaches is to reduce the effect of the depolarizing field created due PAGE 84 84 to the diff erence in the charge distribution of the top and bottom surface of the ferroelectric slabs. For approach (i) Sai et al [103] and Kolpak et al. [101] have examined the effect of Pt and SrRuO3 electrodes on PT thin films. The Pt electrodes stabilized the polarizatio n in PT thin films by screening the depolarizing field. On the other hand, SrRuO3 electrodes partially compensated the surface charges, resulting in a polarization that is only half that of the bulk value in the interior of the thin film. Moreover, for the PT system with SrRuO3 electrode, the presence of ferroelectricity depended on whether the structure of SrRuO3 electrode was equilibrated or not [103] The second approach to achieve short circuit electrical boundary condition in ferroelectric thin films was investigated by Mayer and Vanderbilt [102] The results show that the application of an external electric field to compensate the depolarizing field led to reasonable BaTiO3 thin film structures, whereas it destroyed the ferroelectricity in PT. With thes e significant technical issues associated with the shortcircuit boundary conditions being unresolved even for electronic structure calculations, the current study focused on the opencircuit electrical boundary condition for ferroelectric PT thin films In this boundary condition, the film is not attached to any electrodes. Hence, the electric field of the film is not compensated as is observed during short circuit electrical boundary conditions. The open circuit electrical boundary condition can al so be alternatively mentioned as free standing film in the report. 4.3 .3 Dipole Effect in Symmetric and Asymmetric Slabs In PT thin films, the combination of surface terminations can result in symmetric and asymmetric slabs. Ideally for thick enough thi n film system, symmetric and asymmetric slabs should result the same surface characteristic for individual terminations. On the one hand, the symmetric slabs (PbO PbO and TiO2TiO2) have no net dipole, but their nonstoichiometry can PAGE 85 85 potentially influence the atomic displacements, thereby affecting the ferroelectric property. On the other hand, asymmetric slabs are stoichiometric, but contain a dipole moment perpendicular to the slab due to the presence of charged surfaces [43] This dipole moment can lead to spurious dipole dipole interactions in the direction normal to the surface due to the periodic boundary conditions. For the DFT calculations, the dipole correction [104, 105] available in VASP was used to correct for these interactions. In order to quantify the differen ces between symmetric and asymmetric slab configurations, DFT calculations with 1 x 1 x 7.5 unit cells for symmetric slabs and 1 x 1 x 8 unit cells for asymmetric slabs were performed and characterized (see T ables 4 2 and 4 3). For all the ter minations, th e symmetric and a symmetric slab calculations show similar qualitative and semi quantitative trends for surface structure. The effect of dipole correction in asymmetric slabs for all the terminations show better structural agreement with symmetric slabs in interlayer distance and surface rumpling description (Table 4 2). Comparable interlayer distances and surface rumpling for symmetric and asymmetric slabs in cubic PT have been reported by Piskunov et al [43] using hybridDFT method The interatomic int eraction results, discussed in S ection 4.3, show for sufficiently thick films the choice of symmetric versus asymmetric slabs does not affect the re sults. Table 4 2 Surface characterization for PbO terminations using DFT (in % of c lattice parameter) Slab Type PbO In termination PbO Out termination d 12 d 23 S d 12 d 23 S Symmetric 4.00 0.01 5.31 10.60 6.80 11.97 Asymmetric 2.73 0.24 3.85 10.09 5.36 11.03 Asymmetric, with dipole correction 3.09 0.27 4.37 11.31 6.05 12.20 Cubic slab 6.89 3.07 3.51 6.89 3.07 3.51 Performed for Cubic PT with PbO termination [43] PAGE 86 86 Table 4 3. Surface characterizations for TiO2terminations using DFT (in % of c lattice parameter) Slab Type TiO 2 In termination TiO 2 Out termination d 12 d 23 S d 12 d 23 S Symmetric 8.98 5.89 2.13 6.70 2.38 10.75 A symmetric 7.67 5.12 1.72 6.25 1.20 8.57 Asymmetric with dipole correction 8.60 5.74 2.06 6.99 1.35 9.60 Cubic slab 8.13 5.32 3.12 8.13 5.32 3.12 Performed for Cubic PT with TiO2termination [43] 4.4 Characterization of Surfaces 4.4 .1 Atomic Displacements The atomic relaxation in each layer is calculated by examining the displacement of each ion from its equilibrium position in the b ulk ferroelectric. In the bulk, the cations and anions are displaced from their centrosymmetric position in each layer to develop ferroeletricity. Thus, the Ti and O in a TiO2 crystallographic plane ; and the Pb and O in a PbO crystallographic plane do not actually lie in a single physical plane. In order to characterize the modifications of the crystal structure at the surfaces, the surface displacements of the ions are determined with respect to the positions they would have in a bulk terminated ferroelect ric phase. D ue to the planar symmetry of these films t he displacements of all the atoms of a single species are the same in each layer (i.e., each atomic x y plane) Figure 4 4 shows the displacement of each atom along the z direction for all the four sur face terminations. Cells 1, 2 and 3 represent the surface and next two sub surface unit cells in the PT thin film respectively. The maximum displacement is observed at the surface unit cell (Cell 1) for all the terminations. Also, all the atoms at the surf ace relax into the bulk (positive values in displacement) for all surface terminations with the exception of oxygens on the TiO2 plane for PbO out terminations PAGE 87 87 Figure 4 4 Cell by cell atomic displacement of each individual atom type from t heir original relaxed polarize d condition with Potential 2. A) PbO In, B) PbO Out, C ) TiO2In, and D ) TiO2Out surfaces. Cell 1, Cell 2 and Cell 3 represents the surface unit cell and the next two sub surface unit cells in the thin film respectively. Posit ive and negative values of displacement correspond to relaxation into and out of the surface respectively. 4. 4 .2 Surface Rumpling and Interlayer Distance Figure 4 5 illustrate the change in the interlayer distance (see Figure 4 2 (b)) in the out ermost planes and the change in the planar rumpling at the surface for PbO In and Out polarizations. For the surfaces, the DFT and both empirical potentials c alculations show that the surface relaxes inward, as evidenced by the negative values for d12. The predict ed magnitude of the change in the first interlayer distance ( d12) is similar for both the potentials and the DFT 0.20 0.10 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 Cell 1 Cell 2 Cell 3 Displacement () Pb O_Pb Ti O_Ti 0.20 0.10 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 Cell 1 Cell 2 Cell 3 Displacement () Pb O_Pb Ti O_Ti 0.20 0.10 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 Cell 1 Cell 2 Cell 3 Displacement () Ti O_Ti Pb O_Pb 0.20 0.10 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 Cell 1 Cell 2 Cell 3 Displacement () Ti O_Ti Pb O_Pb PbO In TiO 2 Out Ti O 2 In PbO Out A B C D PAGE 88 88 calculations. A similar level of agreement also holds for the surface rumpling (s ). Similar analyses were done for the TiO2terminated surf ac es and are shown in Figure 4 6 The worst agreement is observed for d23. Only in the case of TiO2Out termination do t he potentials predict an opposite sign from the DFT results. H owever, in all the cases d23 is rath er small. Thus, the surface structures p redicted by the potentials are qualitatively very similar to the DFT results. Figure 4 5 Interlayer distance and surface rumpling calculated for f erroelectric PT with surface, A ) PbO In termin ation and B ) PbO Out termination. The relaxation is giv en in percent of c lattice parameter. Figure 4 6 Interlayer distance and surface rumpling calculated for ferroelectric PT with surface, A) TiO2In termination and B ) TiO2Out termination. The relaxation is given in percent of c lattice param eter. A 15 10 5 0 5 10 15 Relaxation (in % of clattice parameter) Potential 1 Potential 2 DFT d12 d23 s PbO In B 15 10 5 0 5 10 15 Relaxation (in % of clattice parameter) Potential 1 Potential 2 DFT d12 d23 s PbO Out A 15 10 5 0 5 10 15 Relaxation (in % of clattice parameter) Potential 1 Potential 2 DFT d12 d23 s TiO 2 In B 15 10 5 0 5 10 15 Relaxation (in % of clattice parameter) Potential 1 Potential 2 DFT d12 d23 s TiO 2 Out PAGE 89 89 For both terminations, the surface rumpling for Out polarization surfaces are higher than the corresponding In polarization cases. This is due to the larger cation displacement into the bulk compared to the oxygens for the Out polarization surfaces (F igures 4 4 (b ) and 4 4 (d )) than the In polarization surfaces (F igures 4 4 (a) and 4 4 (c)). 4.4 .3 Surface Polarization In addition to the structural details, the atomic level simulations can provide valuable information on the change in polarization fr om the bulk to the surfaces of the se films. Since polarization calculation requires a definite unit cell volume to be defined, the unit cell s used for the polarization calculation s depend on the nature of the surface. For PbO terminated surfaces, a Ti cent ered unit cell is analyzed (Figure 4 7 (a)), and for TiO2terminated surfaces, a Pb centered cell is analyz ed (Figure 4 7 (b)) The polarization values at the bulk like interior of the film are expected to be independent of the choice of the Ti centered or Pb centered unit cell In this way it is possible to de termine the polarization in a physically sensible manner through the entire thickness of the film. Figure 4 7 Schematic of a A) Ti centered unit cell and B) Pb centered unit cell used for p olarization calculation. = Pb = Ti = O A B PAGE 90 90 4.4 .3.1 Polarization Calculation for DFT In order to further validate the empirical potential predictions, it is necessary to calculate the cell bycell polarization for the DFT simulations In order to extract polarizations with u nit cell resolution, the method of Meyer and Vanderbilt [106] is followed In this method, the final optimized output structure and the Born effective charges *Z for tetragonal PT, which we obtained from Zhong, KingSmith and Vanderbilt [107] are considered. The polarization of each unit cell is then calculated as i c iu Z w e P* (4 1 ) w here, e is the electron charge, C is the volume of one unit cell, u is the displacement of atom from its centro symmetric position, w is the corresponding weighting factor, corresponding to the number of atoms of th at type in the unit cell (for example, wfor Pb atoms in a Ti ce ntered unit cell (see Figure 4 7 (a)). Before characterizing the thin films, a polarization calculation was conducted on an optimized bulk PT ferroelectric structure. From the relaxed a tomic positions, a bulk polarization of 81.0 C/cm2 was calculated using Equation 4 1, which matches extremely well with the published DFT value of 81.2 C/cm2 [106] To fully characte rize the surface polarization, F igure 4 8 provides a cell by cell profile o f the polarization normal to the PT thin film surfaces for thick fil ms with PbO Out/PbO In and TiO2Out/TiO2In surfaces using both the interatomic potentials For comparison, t he results obtained from DFT calculations with a thinner film (1 x 1 x 7.5 ) are presented in the same plot Since the system size calculated with D FT was small, polarization values from the top three unit cells nearer to the each surface are plotted for the respective termination in Figure 4 8. PAGE 91 91 Polarization values for four unit cells closest to the surface are shown for respective terminations. Figure 4 8 Cell by cell c direction polarization obtained for a 6 x 6 x 32.5 tetragonal PT system with Potential 1 and Potential 2 where, A ) represents the polarization variation for PbO Out/PbO In and B ) represents the pola rization variation for TiO2Out/TiO2In symmetric slabs. The 1 x 1 x 7.5 unit cell DFT results for respective terminations is presented for comparison The empirical potential results indicate that the films are thick enough to show bulk like behavior in t he middle of the film. The results indicate that the polarization in the Out cases are lower compared to the corresponding In cases, and for thick films its independent of the surface terminations on the other side of the film. A PbO Out PbO In 0 10 20 30 40 50 60 70 80 0 10 20 30 Layer # Polarization (C/cm) Potential 1 Potential 2 DFT B TiO 2 Out Ti O2In 0 10 20 30 40 50 60 70 80 0 10 20 30 Layer # Polarization (C/cm) Potential 1 Potential 2 DFT PAGE 92 92 The interatomic potential re sults show that the polarization profile for each surface is independent of the type of surface (PbO or TiO2) present on the other side of the film which is clear by achieving the bulk polarization away from the surface for all the thin films. For all of these relatively thick films, the polarization in the interior of the film converges to the bulk single crystal values (60 C/ c m2 for Potential 1 and 69 C/cm2 for Potential 2). However, for the films studied with DFT, the polariz ation in the center layer only reaches 90% of the bulk value (~73 C/ cm2 vs. 81.0 C/cm2) It is interesting to note that the magnitude of the polari zation is reduced at all of the surfaces, with different surfaces show ing decreases of different magnitudes (Figure 4 8) The polarization is more strongly suppressed for Out than In polarization cases, which is consistent with the first principles predic tion for PT thin films [102] Simi lar results were observed for Pb O Out/ TiO2In and TiO2Out/ PbO In surfaces. 4.4 .3.2 Polarization Contribution due to Surface Rumpling The analysis of the surface polarization is made more complicated by the surface rumpling. Even for a conventional charge neutral ionic surface, such as the NaCl (001) surface, rumpling leads to a dipole moment developing polarization, in the su rface region [108 110] Although this rumpling leads to an electric polarization, it is of a qualitatively different nature from the ferroel ectric polarization because it is part of the inherent structure of the surface and, thus, cannot be switched by the application of an electric field. The presence of such a rumplingassociated polarization at the PT (001) surface was reported by Sepliarsky et al. [54] in the simulation of very thin films. In particular, for a film with four unit cell thickness th ere was no polarization observed in the interior of the film However, there is an inward polarization of ~3 C/cm2 for PbO terminated films at least some of which presumably arises from the rumpling. PAGE 93 93 Figure 4 9 Polarization induced due to surface rumpling for tetragonal PT with no ferroelectricity for PbO termination (left hand side) and TiO2ter mination (right hand side) with both the potentials. Positive values of polarization indi cate total polarization along + z (into the surface, i.e., In polarization for PbO termination), while negative values indi cate total polarization along z (into the surface, i.e., In polarization for TiO2termination). Inset shows the schematic of the polarization developed due to rumpling. To est imate the effects of rumpling for the PT films, simulations were carried out wit h the in plane lattice parame ters of the f erroelectric phase. For these simulations, the shell model was turned off, thereby suppressing the ferroelectricity. As shown in F ig ure 4 9 this surface relaxation results in a polarization pointing into the film irrespective of the surface terminations. The magnitudes of these polarizations are less than 10 C/cm2 for all surfaces. It is reasonable to assume the polarization in the ferroelectric films that can be associated with their surface rumpling is of a similar magnitude. These inward polarizations due to the rumpling contribute to the total polarizat ion in the surface regions of the ferroelectric films in different ways, depending on the direction of the ferroelectric polarization. The polarization due to rumpling tends to reinforce the ferroelectric polarization for surfaces with an Inpolarization. 6 4 2 0 2 4 0 5 10 15 20 25 30 35 Layer # Polarization (C/cm) Potential 1 Potential 2 PbO TiO 2 y z x PbO TiO 2 PAGE 94 94 Where as, the polarization due to rumpling tends to counteract the ferroelectric polarization for surface with an Out polarization. Thus, the In polarization surfaces are expected to show a smaller deviation of the polarization from the bulk values than th e Out polarization surfaces. This is exactly what is observe d in F ig ure 4 8 4.4 .4 Surface Energy In addition to the above surface characterizations, it is of interest to calculate the surface energies in order to measure the relative stability of a surfa ce. The surface energies for all possible surf ace terminations have been calculated For the empirical potential calculations, the surface energies wer e obtained by a direct comparison between simulations of bulk single crystals and systems with surfaces. However, it is not that straightforward for DFT calculations. T he procedure developed by Heifets et al. is followed [111] t o calculate the surface energies using DFT The energies are divided into two parts, (i) cleavage energy and (ii) relaxation energy. First, the cleavage energies for unoptimized symmetric slabs ( PbO In/PbO Out and TiO2In/ TiO2Out ) are calculated. Since cleaving a PbO In surfac e simultaneously generates a TiO2Out surface and a PbO Out surface generates a TiO2In surface, the relevant cleavage energy is distributed equally between the respective surfaces. Therefore, the cleavage energy can be given as [111] : bulk unrel slab unrel slab unrel sE N TiO E PbO E E 2 ) ( ) ( ) (4 1 (4 2 ) where, PbO Eunrel slab ) ( and 2 ) (TiO Eunrel slab includes the unoptimized In and Out polarization slab energies for PbO and TiO2 respectively, bulkE is the energy per bulk unit cell and N is the number of bulk cells. The relaxation energies for PbO ( In and Out) and TiO2 ( In and Out) are calculated by comparing the relaxed and unrelaxed energies: PAGE 95 95 0 1 2 3 4 5 Potential 1 Potential 2 DFT Surface energy (eV/Surface cell) PbO TiO2 A E A EA Eunrel slab slabrel ) (2 1 (4 3 ) w here A is either the PbO or TiO2 termination Thus, t he surface energy for each termina tion is then calculated as the addition of the cleavage and relaxation energies. A E E A Erel unrel s s ) ( (4 4 ) This surface energy calculated with Equation 4 4 for each termination represents the average energy of the In and Out polarization. Figure 4 1 0 compares the average surface energies for each termination ( In and Out) as calculat ed w ith the potentials and DFT. In each case, it is seen that there is little difference between the Pb O In/PbO Out average energy and the TiO2In/TiO2Out surface energies. Potential 2 and the DFT predict similar average surface energies with opposite order o f stability. However, since the differences are rather small this is not a major concern. Potential 1 predicts the same trend in the average energies for both terminations as DFT, but with a considerably higher value. Figure 4 10. Comparison of a verage In and Out polarization surface energies for PbO and TiO2termination in tetragonal PT. PAGE 96 96 Sepliarsky et al. [54] observed an inplane reconstruction that leads to a polarization parallel to <110> i n ultrathin films (fewer than 8 layers) using Potential 1, with an energy decrease of 0.002 eV/unit cel l. This energy difference is significantly smaller than the energy differences between the In and Out polarization surfaces seen in this study, and therefore not considered Figure 4 11. Cell bycell excess energy per unit cell compared to the bulk PT energy for asymmetric slabs calculated with Potential 2, where A) PbO In/TiO2Out, and B) PbO Out/TiO2In terminations. Similar energy profile s were obtained with Potential 1. The excess energies are added to obtain the surface energy for each te rmination. The average energies provide an overall contribution of both In and Out polarizations. The individual energy contributions can give more details regarding the stability of the In and Out polarized surfaces. The DFT approach can only provide the average energy for the In and Out polarization of a specific termination. However, the interatomic potential simulations can separate the individual surface energy contributions for the In and Out polarizations. Figure 4 11 shows the cell bycell e xcess energy (in comparison to the bulk unit cell energy) along the z direction for asymmetric slabs (Figure 4 11 (a) shows PbO In/TiO2Out and Figure 411 (b) 0.5 0.0 0.5 1.0 1.5 2.0 0 10 20 30 Cell # (EEbulk) (eV/unit cell) 0.5 0.0 0.5 1.0 1.5 2.0 0 10 20 30 Cell # (EEbulk) (eV/unit cell) PbO In PbO Out TiO 2 In Ti O 2 Out A B PAGE 97 97 0 1 2 3 4 5 Potential 1 Potential 2 Surface energy (eV/surface cell) PbOIn PbOOut TiO2In TiO2Out shows PbO Out/TiO2In terminations) using Potential 2. The excess energy in the interior of the film is zero, indicating the energy per unit cell converges to that for the bulk single crystal. Therefore, the surface energies for each individual termination is calculated by summing the excess cell bycell energies from the relevant surface to a region in the interior where the single crystal energy is recovered. Similar results are observed for symmetric slabs. Similar results are observed for Potential 1. Figure 4 12. Individual In and Out polarization surface energies predicted from both t he empirical potentials. Both the potentials predict Out polarized surfaces to be more stable than the In polarized surfaces. The surface energies calculated for each termination using the potentials are shown in Figure 4 12. The trends in the surface e nergies for each termination are the same for both potentials. In particular, both potentials predict that Out polarization surfaces are considerably lower in energy than the corresponding Inpolarization surfaces The stable surfaces predicted by Potentia l 1 and Potential 2 are the TiO2Out and PbO Out surfaces respectively. Unfortunately, as discussed above the se energies cannot be extracted separately from the DFT calculations for PAGE 98 98 validation. However, the good agreement for the average energies shown in f ig ure 4 1 2 gives the results for Potential 2 significant credibility. This raises an interesting question that why the Out polarization surface s are so much lower in energy than the In polarization configurations. This can be justified with the relative signs of the atomic displacements due to the surface rumpling and its contribution to the polarizat ion. As discussed in Section 4.4 .3 the coupling of surface rumpling with surface relaxation in the ferroelectric films reduces the polarization in the Outpolarization surfaces more compared to the In polarization surfaces for each termination. Also, the atomic relaxations indicate that the larger cation relaxation towards the bulk (PbO Out an TiO2Out terminations in Figures 4 4 (b) and 4 4 (d)) i s more favo rable in stabilizing the surfaces than the larger anoin relaxation towards the bulk (PbO In an TiO2In terminations in Figures 4 4 (a) and 4 4 (c)) for the individual terminations. 4.5 Discussion on Stability of Free Standing Films In order to predict the energetically favorable free standing films, different surface termination and polarization combination s have to be considered. T he energies of the Out polarization surfaces are calculated to be much closer to each other with both potential s (Figure 4 12) (PbO Out and TiO2Out surfaces are 2.55 and 2.40 eV/surface cell for Potential 1 respectively, and 0.77 and 0.93 eV/surface cel l for Potential 2 respectively). Therefore, only asymmetric slab configurations (PbO TiO2) are used in this discussion. Three di fferent distinct polarization conditions are considered for the surface energy calculations. 4.5 .1 Normal Polarization (In Out Polarization) This is the polarization condition calculated so far. The idea is to find the minim energy surface terminations w ithin an asymmetric slab configuration. From Figure 412, t he total surface energies of the most stable asymmetric slab configuration are 6.25 eV/surface cell for PAGE 99 99 Potential 1 ( with PbO In/TiO2Out configuration ) and 3.06 eV/surface cell for Potential 2 ( wi th PbO Out/TiO2In configuration ) For DFT calculations the minimum total energy is considered as the sum of the average PbO and TiO2terminated energies (2.96 eV/surface cell). The other possible polarization conditions considered for the asymmetric sla b are the Tail Tail and Parallel configurations. Schematic of both the configurations are shown in Figure 4 13. Figure 4 13. Schematic of A ) asymmetric TailTail ( Out Out ) polarization condition, and B ) asymmetric parallel polarization co ndition in free standing ferroelectric films. For symmetric slabs the terminations on both sides will be either PbO or TiO2, which are not shown here for clarity. 4.5 2 Tail Tail Polarization (Out Out Polarization) From the surface energy calculations in S ection 4.4.4 it is clear that for free standing films, the PT system would strongly favor having the polarization point ing out of each surface. In order to attain such a structure, the film would have to have an interface between an up and down polarize d domain running parallel to the surface (Figure 4 13 (a)). Physically, this is extremely unlikely, as it corresponds to dipole moments of opposite directions being end to end (defined as the Tail Tail structure for discussion) Therefore, the total surfac e energy for these configurations will be the total surface energies calculated for each Out terminations plus the A B PAGE 100 100 0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 Potential 1 Potential 2 DFT Energy (eV/Surface Cell) Normal TailTail Parallel domain wall energy. The domain wall energies are calculated t o be 2.02 and 2.28 eV/surface cell for Potential 1 and 2 respectively. Thus the to tal energy required for two Out polarization surfaces plus an interface in an asymmetric slab are 6.97 and 3.9 8 eV/surface cell for Potential 1 and 2 respectively. Figure 4 14. Total surface energy of PbO TiO2 terminated PT thin film for different termination polarization combinations with empirical potentials and DFT simulations. Parallel polarization condition shows the lowest surface energy for all the simulations. The DFT result for Tail Tail configuration is not calculated. 4.5 .3 Parallel Polarization Another possibility for the polarization configuration is that the polarization in such systems does not lie normal to the film surface. In particular, the films considered for the calculation contained polarization parallel to the s urface ( F ig ure 4. 1 3 (b) ). This case avoids both the possibility of In and Out polarizations as well as development of interfaces as in the Tail Tail cond ition described in Section 4.5.2. The surface energies for the PbO TiO2 surface terminations were cal culated to be 5.32 and 2.49 eV/surface cell for Potential 1 and 2 respect ively. The DFT calculations predict a surface energy of 1.95 eV/surface cell These PAGE 101 101 calculations are consistent with previous DFT results for polarization parallel to the surface [100] The surface energies for different polarization configurations are summarized in Figure 4 14. It is clear that for free standing films, polarization parallel to the surface is the most favorable configuration. Even though the Tail Ta il configuration has two Out polarization surfaces, this configuration is the least favorable due to the creation of an interface. It is important to note that experimentally it is difficult to obtain a free standing film. This is because the thin film is synthesized on a substrate and the interfacial strain state guides the stability of the polarization direction. 4.6 Conclusions The interatomic potentials and DFT calculations provided an atomic level understanding of the combined effect of surface relaxation and polarization in PT thin films. Both the potentials lead to qualitatively similar results for the structure, polarization and energetics of these thin films. The predicted results are consistent with DFT calculations. The surface characterization s predict that the surfaces with Out polarization configuration are more stable than the In polarization surfaces. However, for free standing films polarization parallel to the surface is observed to be more favorable under open circuit electrical boundary conditions A nother important objective of this study was to assess the fidelity of these interatomic potentials for surfaces, which will be used to examine the interaction of domain walls and surfaces. The consistent qualitative prediction for both t he potentials gives some confidence that the potentials can be used successfully to elucidate more complex behaviors in these ferroelectric materials. PAGE 102 102 4 .7 Future Work All the calculations mentioned in S ections 4.4 and 4.5 were performed at 0 K. Therefore, attempts were also made to predict the effect of temperature on ferroelectric behavior in the presence of surfaces with both the potentials. However, both the potentials failed to maintain the surface structure and the system was unstable. The most impor tant observation made for the cause of this failure was the relative core shell displacement. In particular the shells of the Pb atoms near the surface always move to unacceptable ranges (more than 0.2 ), thereby increasing the energy of the system drasti cally, which eventually resulted in unstable structures. Potential 2 has a penalty function to maintain the core shell displacement within 0.2 for each type of atom in PT. This penalty function has been changed to at least three orders of magnitude to m aintain the core shell displacement with the application of temperature. However, all the attempts so far remain unsuccessful. A detailed analysis is necessary for surfaces with temperature to guide develop better sets of interatomic interactions, which wi ll be capable of predicting the effect of temperature. In addition, the surface simulations performed in this work are under open circuit electrical boundary conditions. For a better description of the experimental set up, it is necessary to investiga te the effect of substrates and electrodes. PAGE 103 103 CHAPTER 5 EFFECT OF (100) 180 DOMAIN WALLS IN LEAD TITANATE 5.1 Introduction After investigating the surface energetics and polarization variation in ferroelectric thin films, this chapter focuses on the effec t of 180 domain walls in bulk PT. Domains and domain walls are central to applications of ferroelectric materials for memory applications [112] which demands a clear understanding of domain wall structure, polarization switching and domain dynamics, which are still unclear. There have been considerable experimental [113120] theoretical [121123] and computational [14, 106, 124128] efforts to characterize domains in ferroelectric materials. Among computational approaches, first principles calculations have been widely used to addr ess the atomistic nature of the domain walls in ferroelectric materials [106, 124127] In order to investigate the influence of 180 domain walls in PT, (100) and (110) domain wall terminations are considered. The structural characterization and energetics are addressed using DFT. Table 5 1 shows the simulation details use d for this study. Two different pseudo potentials are used to describe the domain wall characteristic in PT. Table 5 1. Details of DFT simulation approaches for PT domain wall structures Simulation method Potential description Method addressed in the tex t as System size (unit cells along X, Y and Z) Other details DFT Ultra soft pseudo potential (USPP) with LDA USPP 1 x Y x 1 Constant volume Optimization DFT Projector Augmented Wave (PAW) with LDA PAW 1 x Y x 1 Constant volume Optimization Y = 6, 8, 10, 12, 14, 16 5.2 Types of Domain Walls in PbTiO3 The cubic to tetragonal phase transition at the Curie temperature typically leads to the development of different types of domain in PT. Based on the crystal structure, the tetragonal PAGE 104 104 phase has a poss ibility to develop polarization in six different ( x y and z ) directions. A combination of these polarization directions leads to the presence of either 180 or 90 domains in PT (Figure 5 1) Figure 5 1 Schematic representation of the generation of 90 and 180 domain walls during the cubic to tetragonal phase transition below Curie temperature. Based on the structure, the 180 domain walls does not have any lattice mismatch at the domain wall, where as the 90 domain walls have a lattice mismatch. In PT the experimental = Pb = Ti = O 180 Domain P P P Tetragonal (ferroelectric) Cubic (paraelectric) < T c 90 Domain x y z PAGE 105 105 lattice mismatch is nearly 6%. This creates a lot of strain in the system. DFT studies of 90 domain walls in PT have been carried out by Meyer and Vanderbilt [106] This study focuses on the (100) and (110) oriented 180 domain walls in PT. 5.2.1 Pb centered and Ti centered 180 Domain Walls Figure 5 2 Schematic representation of A) Pb centered and B) Ti centered (100) 180 domain wall. The middle plane in each schematic indicates the domain wall. A P Pb centered 180 Domain Wall x y z P x y z P Ti centered 180 Domain Wall P B = Pb = Ti = O PAGE 106 106 DW x y z a c A Y X B Y Z C X Z D Figure 5 3 Schema tic of A ) (01 0) oriented Pb ce ntered 180 domain wall showing the positions of the Pbatoms. Only Pb atoms are shown for clarity. B) X Y projection, C) Y Z projection, and D ) The X Z projection of the domain wall (DW) The bold line in each schematic represents the location of the doma in wall. Further analysis of the domain wall configurations show that there are few different possibilities for the domain wall to sit. It could be centered at PbO planes (Pb centered, Figure 5  PAGE 107 107 2 (a)), TiO2 planes (Ti centered, Figure 5 2 (b)), or betwee n these two planes. The structure of such domain walls has been previously studied using first principles calculations [106, 124127] 5.2.2 Domain Wall Set up For both Pb centered and Ti centered domain walls, supercells of size 1 x Y x 1 with Y = 6, 8, 10, 12, 14, and 16 are used for this study. In the initial structure, h alf of the PT unit cells are oriented with polarization in the + z direction (up polarization), while the other half are oriented with polarization in the z direction (down polarization). Since periodic boundary conditions are applied in all three spatial directions, this system set up results in a pair of domain walls in the x z plane (Figure 5 3 (a)). In order to maintain the inversion symmetry, the atoms in the domain walls are initially placed at their centrosymmetric positions. 5.3 Simulation Methodolog y All the domain wall simulations are performed with DFT at the level of LDA using the Vanderbilt ultrasoft pseudopotentials (USPP) [64] and PAW (Projector Augmente d Wave) [65] method as incorporated in VASP [93 95] The USPP used for domain simulations is the same as that discussed in Section 4.1 in the context of surface simulations. PAW potentials are considered for this study, since thes e are currently state of the art for DFT calculations and that were not available at the time of the earlier studies. Therefore, the DFT studies of surfaces were limited to the USPP. The PAW pseudopotential s treat Pb 5d and Ti 3p states explicitly as valen ce states. All the cal culations are performed with 550 eV ( ~ 40 Ry) cutoff energies. Since the supercells consist of 1 x Y x 1 perovskite unit cells, a 6 1 6 Monkhorst Pack k point [63] mesh is used for the domain simulations. Convergence is achieved when the force on each atom reaches 0.01 eV/. A conjugate gradient method is used to optimize the dom ain structures. The effect of the pse udopotentials, cutoff energy and optimization method on convergence was been established PAGE 108 108 based on previous studies for the PT systems [96, 124, 129] The zero stress lattice parameters for the bulk tetragonal ferroelectric phase using PAW potentials are a = 3.8 6 6 9 and c / a = 1.0436, and are used for all the subsequent PAW DFT calculations. Experimental and USPP lattice parameters for the tetragonal phase are a = 3.904 and c / a = 1.063 and a = 3.86 and c / a = 1.05 respectively [87] 5.4 Motivation for the Domain Wall Simulations The literature shows that Pb centered (100) oriented 180 domain walls are more stable than the Ti centered domain walls [106, 125] However, the first principles studies on domain walls to date have only considered atomic relaxation parallel to the direction of the bulk polarization (along z direction for Figure 5 3). Even though it might seem appropriate to optimize the structure only along one direction, the presence of domain wall breaks the local tetragonal symmetry. The projection of the domain wall structure on an x y plane shows similar lattice projecti on for all the atoms (Figure 5 3 (b)). The presence of the domain wall does not affect the local symmetry. The yz projection of the domain wall structure shows a local break in symmetry along the z dir ection ( Figure 5 3 (c)), which is parallel to the domain wall (along the direction of the bulk polarization). Therefore, optimization parallel to the domain wall along the z direction is expected. This was exactly the focus of the previous first principles calculat ions. However, the x z plane projection also shows a local break in symmetry, which is perpendicular to the domain wall (along the y direction) (Figure 5 3(d)) There fore, it is important to investigate the possibility of local atomic relaxation in the tw o Cartesian directions perpendicular to the bulk polarization ( x transverse direction and y normal direction). PAGE 109 109 A B C B N Fig ure 5 4 Schematic representation of A) Ising type, B) Bloch type and C ) Nel type walls. A mixed wall will have contributio n from a combination of two or all the three kinds of walls. The Bloch and Nel rotations are represented by B and N respectively. The relaxation of atoms along the transverse and normal directions would change the nature of the polarization across the d omain wall. It is well known that the 180 ferroelectric domain walls show Ising character for polarization variation (Figure 5 4(a)). However, the local relaxation along the transverse and normal direction around the domain wall may develop a mixed Bloch like (Figure 5 4(b)) or Nel like (Figure 5 4(c)) character [112] In particular, polarization in the normal direction ( Pn) indicates a Nel like rotation of the polarization vector PAGE 110 110 (N) in the n z plane. The presence of polarization i n the transverse direction ( Pt) indicates a Bloch like rotation (B) in the t z plane. Also, this may change the predicted domain wall width, thereby changing the threshold filed required for domain wall motion [130] In order to systematically understand the effect of atomic relaxation around the domain wall, thre e different optimization schemes are utilized. 5.5 Structural Optimization Schemes There are three different distinct relaxation schemes considered for the (010) 180 PT domain walls: (i) R Z, (ii) R XYZ, and (iii) A XYZ. In the R Z scheme, the dimensio ns of the simulation supercell are held fixed (R) at values corresponding to the lattice constants calculated from bulk PT. All the atoms are allowed to relax only in the z direction i.e., parallel to the bulk uniaxial polarization. This is the scheme prev iously used in the literature. The calculations are performed for both USPP and PAW potentials. The comparison of the results with the previ ous USPP literature and with each other will provide a great deal of confidence for the other optimization schemes. For the R XYZ scheme, the supercell is held rigid with dimensions kept to the lattice constants calculated from bulk PT. All the atoms can relax in all three spatial directions. Only PAW potentials are used to for this study. This scheme will provide ins ight into the effect of in plane relaxation perpendicular to the bulk polarization direction. Finally the A XYZ is an extension of the R XYZ scheme, in which the dimensions of the supercell are allowed to adjust to minimize the stress along with the all t hree spatial relaxation. For large enough separation between the domain walls, the results predicted by R XYZ and A XYZ schemes should converge. PAGE 111 111 5.6 Model Validation (Characterization of Domain Walls with R Z Scheme) In order to determine the effect of ato mic relaxation in the presence of a domain wall in PT, it is necessary to validate the basic results with the reported values with R Z optimization scheme. This will provide confidence in both the simulation set up and characterization techniques, which wi ll be used for the R XYZ and A XYZ schemes. The work by Meyer and Vanderbilt [106] on PT domain walls used the Vanderbilt ultrasoft pseudopotentials (USPP) [64] Since, this potential is available with the VASP package, the USPP with LDA is used for model validation. 5.6 1 Domain Wall Energies The very first comparison performed for the domain wall simulations are the domain wall energies. The domain wall energies are calculated by A E E EBulk Tot Domain Tot Wall Domain 2 ) (_ (5 1 ) Where, ETot_Domai n is the energy of the system with the domain wall, ETot_Bulk is the energy of the same system without the domain wall, and A is the area of the domain wall. The factor of 2 in the denominator accounts for the two domain walls present in the system (see S ection 5.2.2). Figure 5.5 shows the domain wall energies calculated, using Equation (5 1) are comparable to the published results. Even though the energies for both the terminations are slightly higher (which is not surprising, since the Vanderbilt group use their in house simulation package to carry out all the simulations) the Pb centered 180 ferroelectric domain walls are energetically more favorable than the Ti centered domain walls. This trend is consistent with the literature [106, 125] The Pb and Ti centered domain wall energies are reported to be the minima and maxima in the energy curves for domain wall motion [106] Therefore, the energy barrier for the PAGE 112 112 0 20 40 60 80 100 120 140 160 180 200 PbCentered TiCentered Domain Wall Energy (mJ/m2) Literature USPP domain wall motion is predicted to be 25 mJ/m2 against the literature value of 37 mJ/m2 [106] The work by Meyer and Van derbilt [106] has shown that Ti centered domain wall is the transition state between two Pb centered domain walls. Figure 5 5 Comparison of domain wall energies of Pband Ti centered domain walls using USPP with R Z scheme. The literature values are obtained from Meyer and Vanderbilt [106] Since, Pb centered domain walls are more stable than the Ti centered domain walls, the rest of the model validation characterization focuses on Pbcentered domain walls only. 5.6 .2 Atomic Displacements The z displacement of atoms around the Pbcentered domain walls are calculated with respect to the Pb atom at the domain wall (Figure 5 6). The relat ive displacements are calculated in terms of the z lattice parameter. The largest deviation around the domain wall along [010] is observed for the Pbatoms, with a displacement of 0.45 (~ 11.5% of the c lattice parameter). The literature value for this o ffset for the USPP potential is ~14.5% of the c lattice parameter (0.6 ) PAGE 113 113 Figure 5 6 Relative z displacement of atoms with respect to the Pb atom in the Pbcentered domain wall as a function of the distance from the domain wall for Y = 6 and 8 A) Pb and OIII atoms, and B) Ti, OI and OII atoms The integer distances indicate PbO planes and the half integer dustances indicate TiO2planes. 5.6 3 Ferroelectric Distortion By calculating the change in bonding environment for each atomic pla ne parallel to the domain wall, the ferroelectric distortion ( Rz) around the domain wall provides an estimate of the ferr oelectric order. In general, ferroelectric distortion is calculated as the ratio of the displacement of the metal atoms relative to an oxygen atom in the PbO or TiO2 plane [97, 125] 0.44 0.46 0.48 0.50 0.52 0.54 3 2 1 0 1 2 3 Distance from domain wall [units of a] z [units of c] 6 unit cells 8 unit cells B Ti O II O I 0.10 0.08 0.06 0.04 0.02 0.00 0.02 0.04 0.06 0.08 0.10 3 2 1 0 1 2 3 Distance from domain wall [units of a] z [units of c] 6 unit cells 8 unit cells A Pb O III l ocation of domain wall PAGE 114 114 for the domain structure and the undistorted bulk ferroelectric phase. The formulation of Rz for TiO2plane is explained for clarity. Figure 5 7 Schematic of tetragonal ferroelectric PT unit cell highlighting the Ti bonding. The short ( bs) and long ( bl) bond lengths of Ti with oxygen is utilized to calculate the ferroelectric distortion around the domain wall. The experimental bond lengths are collected from Shirane et al [19] Figure 5 7 shows the bonding environment of Ti in the ferroelectric PT. In order to calculate the Rz for TiO2plane, the short and long Ti O bonds are considered as the reference. Thus, the ferroelectric distortion is defined as s aver new s aver zb b b b R (5 2) Where, bs is the short Ti O bond length in bulk ferroelectric phase (Figure 5 7), bs,new is the short Ti O bond length at a distance d from the domain wall. baver is the average Ti O bond length and is given as ( bs + bl)/2. The value of Rz changes from 0 (for paraelectric phase, i.e., bs,new = baver) to 1 (for bulk ferroelectric phase, i.e., bs,new = bs). b s = 1.78 1.98 b l = 2.38 P x y z = Pb = Ti = O PAGE 115 115 This ferroelectric distortion is shown for the TiO2planes along the z direction for the Pbcentered 180 domain walls. The ferroelectric lattice distortion is observed to be 0. 83 for the first nearest neighbor Ti O2planes around the domain wall and reaches the bulk ferroelectric state beyond it (a value of 1.0 means there is not distortion) No ferroelectric distortion was observed for the PbO planes across the domain wall. A similar characteristic is reported in the literature with the lattice distortion of 0.80 [106] and 0.73 [125] for the first nearest neighbor TiO2planes for Pb centered 180 domain walls. Figure 5 8 Ferroelectric lattice distortion of TiO2planes in the z direction for a Pb centered 180 domain wall. Only Pb atoms are shown for clarity. The location of the domain wall is on the PbO pl ane shown as a thicker line. 5.6 .4 Summary of USPP Calculations All the domain calculations performed with USPP using R Z scheme match with the literature very well. Using the same USPP with LDA description, there were slight quantitative differences betwe en the calculated and Meyer and Vanderbilt [106] results. This is not sur prising, since the calculations performed for this work used the VASP package, while Meyer and Vanderbilt used their in house first principle package. Therefore, the sources of the small differences between these calculations are not easy to identify, though most likely lie in slightly different choices in cut off energy and/or k point mesh. z TiO 2 plane y PbO plane R z : 1.00 0.83 0.83 1.00 PAGE 116 116 0 20 40 60 80 100 120 140 160 180 200 PbCentered TiCentered Domain Wall Energy (mJ/m2) Literature RZ The domain results discussed so far provides a great deal of confidence to carry out further calculations on different relaxation schemes. Since PAW potentials are curr ently the state of the art for DFT calculations, a consistent set of data is calculated for comparing R Z, R XYZ and A XYZ relaxation schemes for 180 domain walls in PT. 5.7 Characterization of Domain Walls with PAW Potential (R Z Scheme) The domain wall energies and cell bycell polarization across the domain wall are calculated for t he (010) 180domain walls in PT by allowing only uni axial optimization (R Z Scheme) of atoms. 5.7 1 Domain Wall Energies Figure 5 9 Comparison of domain wall ene rgies of Pb and Ti centered domain walls using PAW potential with R Z scheme. The literature values are obtained from Meyer and Vanderbilt using USPP [106] Figure 5 9 shows the calculated energies for Pb and Ti centered domain walls for Y = 10, calculated using the R Z scheme with the PAW method. For both walls, the PAW method predicts slightly lower domain wall energies than the literature (USPP) values [106] However, in commont with the USPP calculations, the calculations predict the Pb centered (100) 180 PAGE 117 117 domain walls have lower energy (128 mJ/m2) than the Ti centered domain walls (163 mJ/m2), (Fig. 3). The energy barrier for domain wall move ment is ~ 35 mJ/m2, which is very similar to the value of ~ 37 mJ/m2 previously reported [106] Thus the overall energy comparison is very good for the R Z scheme with PAW potential. 5.7 .2 Polarization Calculation The polarization calculation was performed i n the same way as discussed in Section 4.4 .3. However, for PAW pseudo potentials, the Born effective charges required to calculate the polarization were not available in the literature. Therefore, these values were calculated for each crystallographically distinct ion in the ferroelectric and paraelectric PT, usi ng the Ber ry phase calculation [131] and are reported in Table 5 2. The polarization in the bulk ferroelectric PT is calculated to be 86.7 C/cm2. Table 5 2. Calculated Born effective charges of the optimized Cubic and Tetragonal PT structures as calculated with a PAW potential, and compared with literature values. Atom Z xx Z yy Z zz Cubic Pb 3.91 (3.90 a ) 3.91 3.91 Ti 7.10 (7.06 a ) 7.10 7.10 O I 2.60 ( 2.56 a ) 5.84 ( 5.83 a ) 2.60 O III 2.60 2.60 5.84 Tetragonal Pb 3.84 (3.92 a 3.85 b 3.74 c ) 3.84 (3.85 b 3.74 c ) 3.49 (3.63 b 3.52 c ) Ti 6.37 (6.71 a 6.37 b 6.20c ) 6.37 (6.37 b 6.20 c ) 5.47 (5.41 b 5.18 c ) O I 2.65 ( 2.56 a 2.66 b 2.61 c ) 5.26 ( 5.51 a 5.32 b 5.18 c ) 2.22 ( 2.22 b 2.16 c ) O III 2.24 ( 2.24 b 2.15 c ) 2.24 ( 2.24 b 2.15 c ) 4.64 ( 4.60 b 4.38 c ) a Vanderbilt ultrasoft pseudopotential (LDA) Reference [107] b Norm conserving plane wave pseudopotential (LDA) Reference [132] c Full potential ab initio LAPW method with local orbital extension (GGA) Reference [126] PAGE 118 118 100 80 60 40 20 0 20 40 60 80 100 20 10 0 10 20 Distance from wall () ZPolarization (C/cm2) Y = 06 Y = 08 Y = 10 Bulk Figure 5 10. Cell bycell z polarization variation across the Pb centered domai n wall for different system sizes (Y = 6, 8 and 10) with PAW pseudopotential. For all the system sizes the bulk polarization is achieved in the second unit cell away from the domain wall. Figure 5 10 shows the cell bycell z polarization variation along [010] for different system sizes for Pb centered domain walls. For all supercell sizes, the polarization reaches the bulk value within two lattice parameters from the domain wall, suggesting a very narrow domain wall width [106] The polarization in the first nearest neighbor unit cells around the domain wall reach~ 68.5% of t he bulk polarization.. Fig. 5 10 shows that the polarization profile is fairly independent of the system size, indicating that there is negligible interaction between the domain walls in the simulation supercell. Similar polarization profiles were obtained for the Ti centered domain walls with a polarization value of ~ 56.1% of the bulk polarization in the unit cells proximate to the domain wall. 5.7 .3 Polarization Fitting to Calculate Domain Width The Ginsburg Landau Devo nshire (GLD) formalism is a succes sful phenomenogoical method to describe the ferroelectric behavior of many materials In particular, it has been used successfully used to describe PT systems [22] Considering only the z direction polarization PAGE 119 119 0 20 40 6080 100 120 140 160 180Pbcentered TicenteredDomain Wall Energy (mJ/m2) RZ RXYZterms ( Pz), the solution to a one dime nsional fourth order GLD model for the domain wall profile is given by [106, 133] : Pz Pz otanh( xn/o) (5 3) w here Pz,o is the bulk polarization, xn is the distance from center of the domain wall, and wo is the domain wall width. The domain wall profile in Figure 5 10 was fitted to Equation 5 3 with Pz,o = 86.7 C/cm2 (the bulk polarization value ) and wo = 2.2 0.02 Therefore, the domain wall width is 2 wo = 4.4 0.04 5.8 Characterization of Domain Walls with PAW Potential (R XYZ Scheme) In order to achieve optimization with R XYZ, the atoms in the supercell were allowed to relax in all the thr ee crystallographic directions. 5.8 .1 Domain Wall Energies Figure 5 11. Comparison of domain wall energies of Pband Ti centered domain walls using PAW potential with R Z and R XYZ schemes. The domain wall energies are observed to be lower for the R XYZ scheme than R Z scheme. As for the the R Z relaxation, the simulations predicted that the 180 Pb centered domain walls are more stable in PT than the Ti centered domain walls (Figure 5 11) for R XYZ scheme. PAGE 120 120 100 80 60 4020 0 20 40 60 80 100 40 20 0 20 40 Distance from wall () ZPolarization (C/cm2) Y = 06 Y = 08 Y = 10 Y = 12 Y = 14 Y = 16 The relaxation of the in plane const raint reduced the domain wall energy by 4.7 mJ/m2 (3.7 %) for Pb centered and 3.4 mJ/m2 (2.1%) for Ti centered walls compared to their respective R Z schemes As will be seen below, although this energy difference is small, it results in a significant stru ctural relaxation that breaks the uniaxial symmetry of the ferroelectric phase. 5.8 .2 Polarization Calculation The polarization profile along the z direction through [010] is essentially identical to that observed for the R Z relaxation scheme (Figures 5 1 0 and 5 12) for all the system sizes studied (Y = 6, 8, 10, 12, 14 and 16). Figure 5 12. Cell bycell z polarization variation across the Pb centered domain wall for different system sizes with PAW pseudo potential within R XYZ scheme. Despite t he relaxation in all three spatial directions, there was no polarization in the transverse direction (Pt = 0) to the domain wall, suggesting no atomic movement parallel to the domain wall (the x direction ). However polarization was present in the directio n normal to the domain wall (Pn), which is consistent with the structural explanation described in Figure 5 4. Since both Pb and Ti centered domain walls showed the presence of Pn, each of the domain terminations is discussed here. PAGE 121 121 5 4 3 2 1 0 1 2 3 4 5 60 40 20 0 20 40 60 Distance from wall () YPolarization (C/cm2) Y = 06 Y = 08 Y = 10 Y = 12 Y = 14 Y = 16 5.8 .2.1 Pb Centered Do main Walls For Pb centered domain walls, the polarization calculation along the normal to the domain wall direction developed a polarization (Pn) with a maximum value of 2.4 C/cm2, which is ~ 2.5% of bulk Pz. This polarization points away from the domain wall. The profile of Pn is observed to be independent of the system size and goes to zero away from the domain wall. Figure 5 13. Cell bycell y polarization (Pn) variation across the Pb centered domain wall for different system sizes with PAW pseudo potential within R XYZ scheme. The polarization values are larger near the domain wall and reduce to zero at the bulk like regions of each domain. Domain w idth with Pn p olarization : Two different fitting schemes are used for calculating the domain wall width. Scheme 1: The one dimensional solution of the fourth order GLD model was modified to describe the effect of Pn and Pt. Scrymgeour et al. [133] expressed the polarization components Pn and Pt as a polynomial expansion in the odd powers of the primary order parameter, Pz. Since the polarization in the normal ( Pn) and transverse ( Pt) directions are expec ted to be small, a solution form for the polarization was proposed based on the primary polarization Pz as [134] : PAGE 122 122 ) ( ~ ) (1 2 1maxn k z ik k k n ix P x P (5 4) where, i=n, t, z and k =1, 2, 3. For kmax=1, Pz reduces to the 1 D solution as shown in Equation 5 3. Since Pt = 0 for R XYZ schemes, only Pz and Pn are considered for fitting. For kmax = 3, Equation 5 4 can be written as: ) / ( tanh ) / ( tanh ) / tanh( 0 0 1 / ) ( / ) (0 5 0 3 0 3 2 1 max max ,w x w x w x P x P P x Pn n n n n z z (5 5) Simultaneous fitting of Pz and Pn for Pb centered domain walls based on Equation 5 5 yields wo = 4.660.01 (Figure 5 14), which is almost twice the value ( wo = 2.20.02 ) obtained by fitting Equation 5 3. Thus, considering only a 1D solution for Pz resulted in two different domain wa ll widths for the same wall. This indicated that both Equations 5 3 and 5 5 are not suitable for calculating the domain wall width. From this fitting experience of the domain wall width, the expression in Equation 5 5 is modified to include twice the domai n wall width along the normal direction ( Pn). This is acceptable since the Pn polarization profile reaches a maximum and decreases to zero away from the domain wall, unlike Pz which reaches the maximum value and maintains that throughout the domain. Theref ore, only the Pn contribution to the fitting is modified as: 5 0 5 3 0 3 0 1 max2 tanh 2 tanh 2 tanh w x w x w x P Pn n n n (5 6) Thus, a simultaneous fitting of Equations 5 3 and 5 6 resulted in a domain wall width of ~2.33 which is exactly half of the domain wall width obtained using Equation 5 5 for Pbcentered domain wall. PAGE 123 123 Scheme 2: The other modification of the fitting scheme was performed by including the higher order terms for the Pz component. Thus, the modified expression for the individual polarizations are normalized to the to their corresponding maximum values and can be well fitted by the general GLD prediction as ) / ( tanh ~ / ) (1 2 max o n k ik i n iw x P x P (5 7) A simultaneous fitting of Pz and Pn resulted in a domain wall width of wo~ 4.660.01 for Pb centrered domain wall. The key difference between Equation s 5 7 and 5 5 is that the value of z2 and z3 are non zero in Equation 5 7 The matrix in Equation 5 7 is given by: 025. 0 92 3 030 0 53 7 006 0 61 3 009 0 43. 0 015 0 36 1 006 0 93 1XYZ R ik (5 8 ) confirming that higher order tanh terms are indeed required to describe the Pz. It is important to note that the z k = 1, which is consistent with the constraints that the 1, as z nk = 0 is consistent with the constraint that n 0 as z The presence of the nonuniaxial polarization indicates the presence of a Bolch or Nel like charac ter for the Pb centered domain walls [134] For this boundary, Pt = 0 indicates the absence of Bloch like rotation (B) for the 180 Pb centered domain walls. The Pn polarization on the other had leads to a rotation of the polarization vector (N) in the n z plane. For Pb centered domain walls, the rotation angles are calculated to be B = 0 and N ~ 1.43, which indicates a combination of Ising and Nel type characteristic for these domain w alls. PAGE 124 124 Figure 5 14. Cell bycell normalized polarization fitting performed for both Pz and Pn (system size Y = 16) for Pbcentered domain wall using R XYZ scheme The open symbols are the data points obtained for polarization in each unit ce ll from the optimized structure and the solid lines are the fitting curves. 5.8 .2.2 Ti Centered Domain Walls The polarization around a Ti centered domain wall is observed to be similar to that for the Pb cent ered domain walls discussed in Section 5.8.2. 1. The optimization with R XYZ scheme developed Pn = 1.52 C/cm2 (~1.75% of bulk Pz) at ~ 3.87 away from the domain wall. Using fitting scheme 1 the domain wall width of wo ~ 3.22 is predicted for Ti centered domain walls. However, A simultaneous fitting of Pz and Pn for Ti centered domain walls using Equati on 5 7 resulted in wo = 6.440.41 (Figure 515). Thus, the wall width for Ti centered domains is 12.880.41 This is substantially larger than the P b centered domain wall, which i s 4.66 It is important to note that polarization values around 1.93 from the center of the domain wall are neglected for the fitting. The matrix is given by: PAGE 125 125 471 0 960 4 550 0 070 9 178 0 120 4 381 0 760 0 554 0 950 1 191 0 180 2XYZ R ik (5 7) Figure 5 15. Cell bycell normalized polarization fitting performed for both Pz and Pn (syste m size Y = 16) for Ti centered domain wall using R XYZ scheme. The open symbols are the data points obtained for polarization in each unit cell from the optimized structure and the solid lines are the fitting curves. The rotation angles for Ti centered dom ain walls were calculated to beB = 0 and N ~ 1.00. These analyses indicate that the presence on non uniaxial polarization components increase the effective domain wall width. 5.8 .3 Atomic Displacements The atomic displacements for each type of atom normal to the domain wall ([010]) are calculated with respect to their respective positions in the bulk ferroelectric phase. PAGE 126 126 5.8 .3.1 Pb centered Domain Walls Figure 5 16 schematically represents the atomic relaxation around the Pb centered domain wall. T he [001] view show that the oxygens ( OIII) present around 1.93 away from the domain wall on the PbO plane relax the most, by 0.03 The Ti atoms, also 1.93 away from the domain wall, on the TiO2plane move ~ 0.01 Both these atoms relax towards the d omain wall. Fig ure 5 16. Schematic of the atomic displacement of each type of ion around the Pb centered domain wall. The <001> vie w of A) PbO and B ) TiO2 planes are shown with the arrows indicating the direction of atomic movement. The oxygens on the PbO planes (OIII) moves the maximum (~ 0.03 ), while the Ti atoms move ~ 0.01 Using Equation 4 1, the polarization contribution from each individual atomic relaxation can be separated. The Ti displacement contribute to +1.25 C/cm2 pointing towards the domain wall while the OIII displacement contribute 3.61 C/cm2 to the polarization pointing away from the domain wall. Therefore, the net polarization is 1.25 3.61= 2.35 C/cm2, which is consistent with the magnitude (2.40 C/cm2) and direction (away from the domain wall) reported above. The displacements of the rest of the atoms (Pb, and the oxygens (OI and OII)) are less than 0.005 and therefore did not included for the discussion. = Pb = Ti = O x y PbO plane O III A O I TiO 2 plane O II B PAGE 127 127 1.000 1.010 1.020 1.030 1.040 1.050 Bulk Y = 06 Y = 10 Y = 16 c/a Ratio 5.8 .3.2 Ti centered Domain Walls The relaxations of atoms around the Ti centered domain walls are similar to the Pb centered domain walls. All the atoms move towards the domain wall. The Pbatoms move ~ 0.013 and the oxygen (OIII) atoms on the PbO plane move ~ 0.023 around the domain wall. The displacements are less than 0.006 for the rest of the atoms (i.e., all the atoms on the TiO2 plane). 5.9 Characterization of Domain Walls with PAW Potential (A XYZ Scheme) Fig ure 5 17. Change in tetragonality with increase in system size for 180 Pbcentered domain wall systems using A XYZ scheme. It is clear that increasing system size increases the accuracy to follow the bulk system. In this scheme both the dimensions of the supercell and the atomic relaxation in all possible crystallogra phic directions are allowed. The optimized lattice parameter approaches the bulk ferroelectric structure with increase in the system size (Figure 5 17). The largest system studied for the A XYZ schemes are Y = 16, which is still not a large enough system t o achieve bulk tetragonality. However, the qualitative polarization description is still predicted accurately. PAGE 128 128 5 4 3 2 1 0 1 2 3 4 5 40 30 20 10 0 10 20 30 40 Distance from wall () YPolarization (C/cm2) RXYZ AXYZ Fig ure 5 18. Change in Pn using R XYZ and A XYZ are the the same for Pb centered domain wall with Y = 16. For the Pb centered domain wa ll with Y =16, the converged lattice parameters are a = 3.8731 and c / a = 1.0325 (Figure 5 18), which show smaller tetragonality compared with the bulk values of a = 3.8 6 6 9 and c / a = 1.0436. The polarization calculated in the interior of these domains r eached a value of 80.4 C/cm2, compared with the bulk value of 86.7 C/cm2. The domain wall energy is calculated to be 114 mJ/m2, which is somewhat lower than the value of 122 mJ/m2 determined using R XYZ scheme. Most importantly, however, these calculations also show a profile for Pn that is essentially indistinguishable to that determined for fixed supercell dimensions. 5.10 Discussion Having established the presence of polarization normal to the domain wall in (100) 180 domain walls in PT, the study is extended to character ize (110) 180 domain walls in PT. Based on the stacking sequence of (110) planes, the 180 domain wall can be either PbTiO centered PAGE 129 129 0 20 40 6080 100 120 140 160 180 Pbcentered Ticentered PbTiOcentered OOcentered Domain Wall Energy (mJ/m2) and OO centered. All the calculations were performed with a 1 x 6 x 1 supercell size with PAW potential. 5.10.1 Domain Wal ls Energy Comparison between (100) and (110) Domain Walls Figure 5 19. Comparison of domain wall energies of (100) and (110) domain structures using R XYZ scheme. The Pb centered and OO centered domain walls have comparable energies. The domain wall energies were calculated using both the R Z and R XYZ schemes. For RZ optimization, the energy of the OO centered (110) 180 domain wall is 124 mJ/m2, whereas the energy of the PbTiO centered domain wall is much higher at ~156 mJ/m2. The domain wall energies are reduced by ~3 mJ/m2 for both the (110) domain walls using R XYZ optimization, which is ~2.3 % for OO centered and ~ 1 .9 % for Pb Ti O centered walls From the difference in these energies, the barrier for domain wall motion in the [110] direction is ~32 mJ/m2, which is comparable to the ~35 37 mJ/m2 for [010] direction [106] Figure 5 19 shows the comparison of the domain wall energies between the (100) and (110) 180 domain structures. The (100) Pbcentered and (110) OO centered domain walls have comparable energies and are more stable than the other two domain terminations. PAGE 130 130 100 80 60 40 20 0 20 40 60 80 100 15 10 5 0 5 10 15 Distance from wall () ZPolarization (C/cm2) (100) Pbcentered (110) OOcentered 5.10.2 Polarization Variation Comparison between (100) and (110) Domain Walls Figure 5 20. Cell bycell z polarization variation across the (100) Pbcentered and (110) OO centered domain wall with PAW pseudopotential within R XYZ scheme Figure 5 20 illustrates that the cell bycell z polarization profiles ( Pz) for (100) Pb centered and (110) OO centered domain walls are very similar. Bulk polarization values are achieved approximately two unit cells away from the domain wall. The (110) OO centered domain walls also show a polarization perpendicular to the domain wall ( Pn) with a maximum value of Pn ~ 1.8 C/cm2 (2.1% of bulk Pz) pointing away from the domain wall. Overall, the domain characteristics are predicted to be very similar for b oth (100) and (110) 180 domain walls. 5.11 Critical Assessment of the Inter atomic Potentials All the simulation results discussed so far in Chapter 5 are using first principles calculations. Interatomic potentials can also be used to describe the domain struct ures and energetics. Attempts we re made to simulate the (100) 180 domain wall structures. Interestingly, the interatomic potentials were able to describe the presence of the in plane normal polarization for all the domain wall calculations. However, the interatomic potentials are not sensitive enough PAGE 131 131 0 50 100 150 200 250 300 Literature Potential 1 Potential 2 Domain Wall Energy (mJ/m2) Pbcentered Ticentered to capture the direction of these small in plane polarizations correctly. Also, Potential 1 does not distinguish between (100) Pbcentered and Ti centered domain walls. It predicts the same energy for bo th the domain walls. Potential 2, on the other hand, predicts Ti centered domain wall to be more favorable This does not indicate that the potentials are not useful for domain studies. To provide a better prospective, a difference in 0.0033 eV/PbTiO3 resu lts in 10 mJ/m2 change in the domain wall energy for Potential 2. This precision is very difficult to probe with the current interatomic potentials. Some details of the results obtained through the interatomic potentials are reported here. 5.11.1 Domain Wa ll Energies Figure 5 21. Comparison of domain wall energies of (100) domain walls using Potential 1 and Potential 2. The literature values are obtained from Meyer and Vanderbilt [106] Figure 5 21 compared the energies of (100) 180 domain walls predicted by the two interatomic potentials (Potential 1 and Potential 2) with the literature first principles calculations [106] The interatomic potential description with Potential 1 fails to predict the relative stability of the domain walls. There is an equal probability for both Pb centered and Ti centered domain walls possible with Potential 1. Potential 2 on the other hand predicts the oppos ite stability for PAGE 132 132 the domain walls. The Ti centered domain wall is stable over the Pb centered domain walls by ~ 25 mJ/m2. Thus, care should be taken while using these interatomic potentials for domain wall characterizations. Current work by Shimada et al [14] developed a new set of interatomic interactions to describe the domain walls in PT. The potential can reproduce the domain wall energy of the (100) 180 Pbcentered domain wall energies of 132 mJ/m2, which is the same as the first principles calculation results. However, the Shimada potential fails to stabilize the Ti centered domain walls, hence no energy comparison is possible. Also, this potential works with the DFT derived lattice parameter, and not with the optimized lattice parameter predicted by the potential. 5.11.2 Polarization Analysis Figure 5 22. Cell bycell polarization normal to the domain wall (Pn) predicted with A) Potential 1, and B) Potential 2. Both the results show polarizations pointing towards the domain wall. It is interesting to note that both the interatomic potentials predict the presence of polarization normal to the domain wall (Pn) a nd the absence of polarization in the transverse direction (Pt), which is consistent with the DFT calculations. However, for Pn, the predicted direction of the induced polarization is towards the domain wall, which is exactly the opposite 0.6 0.4 0.2 0 0.2 0.4 0.6 40 30 20 10 0 10 20 30 40 Distance from wall () YPolarization (C/cm2) Y = 06 Y = 08 Y = 10 Y = 12 Y = 14 0.8 0.6 0.4 0.2 0 0.2 0.4 0.6 0.8 40 30 20 10 0 10 20 30 40 Distance from wall () YPolarization (C/cm2) Y = 06 Y = 08 Y = 10 Y = 12 Y = 14 A B PAGE 133 133 direction of what the DFT calculations predicted (Figure 5 22). Further analyses of the polarization values indicate that even for Y = 14 unit cells, the polarization never reaches the bulk value. It passes through zero at the center of each domain to maintain the symmetry. Moreover, the Pn values predicted by both the interatomic potential are nearly 3 times smaller than the DFT values. 5.12 Conclusions First principles calculation of the 180 domain wall in PT show the presence of nonuniaxial polarization normal to the domain wall ( Pn). A simultaneous fitting of Pn and Pz using modified GLD theory indicate a domain wall width almost twice than what is normally calculated with uniaxial polarization only. This modified domain wall width can directly affect the threshold filed required for domain wall motion [130] The comparison of the energies indicate an equal probability of the (100) 180 Pbcentered and (110) 180 OO centered domain walls in PT. 5.13 Future Work The presence of the in plane polarization due to the domain wall is predicted by both the atomic level simulation m ethods. However, the interatomic potentials predict an opposite in plane polarization direction compared to the first principles calculations. Hence, better interatomic potential description is necessary to accurately describe the domain walls in PT. The e ffect of temperature can be included to study the domain dynamics using interatomic potentials. PAGE 134 134 CHAPTER 6 STRUCTURE AND PROPERTIES OF 5 (310) TILT GRAIN BOUNDARIES IN STRONTIUM TITANATE A ND LEAD TITANATE 6.1 Effect of Grain Boundaries on Ferroelectrici ty Investigations of the effects of grain size on ferroelectric behavior show significant deviations from the bulk properties, followed by the suppression of ferroelectricity below a critical size [135, 136] Strukov et al [135, 137] characterized the grain size effect by synthesizing various BT thin films w ith grain sizes varying between 35 165 nm, which displayed a decrease in the cubic to tetragonal phase transition temperature with decreasing grain size. Other experimental studies on polycrystalline PT thin films also exhibited a smaller spontaneous polar ization and a larger coercive field value when compared with a PT single crystal [90, 138] The increase in the coercive field was attributed to the grain boundaries which restrict the motion of the domains u nder an external electric field. TEM analysis of PT thin films established that the width of the domains decreases with decreasing grain size. Investigations on unsupported thin PT films indicated that the dielectric constant decreases abruptly after some critical size (200 nm) rather than increasing monotonically with decreasing grain size. Also the coercive field increases abnormally below the critical size [115] Decreasing the grain size below the critical size (120180 nm), only a few grain showed the existence of lamellar domains with a domain width in the range of 5 20 nm. Further decreases in the grain siz e (70 90 nm) mostly resulted in single domain grains. Bending of domain walls due to decrease in size was also reported, and had been attributed to the distribution of the strain field, interactions among the domain walls or the charging of the domain wall s [115] In order to address the atomistic details of the effect of grain boundaries on ferroelectricity, th is study focused on simple tilt boundaries. Grain boundaries (GBs) are generated from the mismatch in lattice orientation in adjacent grains. These boundaries can be a combination of tilt PAGE 135 135 or twist. A twist boundary can be generated by rotating one half of the crystal with an angle over the other half about an axis perpendicular to the interface (Figure 6 1(a)). Thus, twist GBs consists of a series of screw dislocations. Tilt GBs can be achieved by the relative rotation of two adjacent lattices along the p lane of the interface. The tilt GBs are symmetric if the atomic planes parallel to the GB on either side is the same (Figure 6 1(b)). Tilt GBs can be described as a series of edge dislocations located in the plane of the GB. Figure 6 1. Schematic o f the relative rota tion of the grains to generate A) twist and B ) symmetrical tilt boundaries [139] GBs are characterized by the energy of the boundaries, Egb as A E E Egb ) ( ) ( (6 1 ) wh ere, E() is the energy o f s tructure with GB E is the energy of the crystal without a GB, and A i s the area of the GB The Egb depends on the rotation angle (). Figure 6 2 illustrates the theoretical calculations of the effect of tilt angle on the Egb in aluminum [139, 140] Both experimental measurement and theoretical calculations predict similar Egb behavior with th e angle of rotation. From Figure 6 2 it is clear that the Egb increases substantially with increase in the rotation angle at = 0, and decreases rapidly at = 90. However, there are additional cusps A B PAGE 136 136 in energy at intermediate rotation angles, which are st ructurally related to some special boundaries. These boundaries are described by the coincident site lattice (CSL) model. Figure 6 2. Computed Egb as a function of rotation angle of symmetrical <110> tilt boundaries in Al [139] 6.2 Coincident Site Lattice (CSL) The CSL model represents the occurrence of coincident atomic positions in the GB with ideal lattice positions in both adjacent grains. Since both grains have periodic lattices, the coincident sites should also be periodic. Therefore, the repeating lattice at the GB is defined as the CSL. The density of the coincident sites or the size of the repeating cell for the CSL is defined as lattice crystal of cell elementary of Volume CSL of cell repeating of Volume (6 2 ) Since the elementary cell of CSL has the same lattice parameter along the normal to the GB, Eq. 6 2 can be modified as PAGE 137 137 lattice crystal of cell elementary of Area CSL of cell repeating of Area (6 3 ) This provides comparison of the two dimensional overlap of the CSL with the original crystal lattice. Along with the value, a rotation ang le and the rotation axis is necessary to define a twist boundary with CSL, whereas an additional GB plane is defined for tilt boundaries. Figure 6 3 illustrates a 36.9 <100> 5 tilt and 36.9 5 twist GB. Further information on CSL formalism is availabl e in [141] Figure 6 3. Schematic of a 36.9 <100> 5 GB structure of A ) tilt boundary, where the GB plane is perpendicula r to the plane of the page and B ) twist boundary, where the GB plane is parallel to the plane of the page [139] A B PAGE 138 138 6.2.1 CSL Boundaries in Perovskites Previ ous work on CSL boundaries in perovskites focused mostly on SrTiO3. Most of the CSL boundary studies addressed the GB structure and defect energetics. There has been significant amount of experimental (transmission electron microscopy (TEM) and STEM), and atomic level simulations to address various tilt GBs, e.g., [001] (310) 5, [001] (210) 5, [110] (111) 3 in SrTiO3 [142145] There is also limited literature on [110] (111) 3 tilt boundaries in BaTiO3 [144] However, there is no report on the CSL boundaries in PT in the literature. Since the main focus of this study is to address the effect of GB on the polarization in PT, it is important to investigate the polarization induced in a paraelectric structure, like SrTiO3. This result will provide the measure of the non switchable part of the polarization contribution due to the presence of the domain wall. The [001] (310) 5 grain boundary in the bicrystal system is investigated first for SrTiO3 and then for PT. 6.3 Bicrystal Simulation S etup All the GB ststems are simulated using a periodic supercell containing two GBs. The simulation set up is schematically represented in Figure 6 4. The GBs are parallel to the xy plane. Rigid body translation should be characterized for the relative mot ion of Grain 1 with respect to Grain 2 along x and y directions. Figure 6 4. Schematic of a bicrystal simulation set up. Due to periodic boundary condition there are only two grains in the system (Grain 1 and Grain 2). GB1 and GB2 represent the two grain boundaries. The GBs are perpendicular to the plane of the paper. Grain 2 Grain 1 Grain 2 GB 2 GB 1 z y PAGE 139 139 Similar to the different relaxation schemes for domain walls, two different optimization schemes are employed to find the stable GB. In the one scheme, only z relaxation of atoms and t he supercell is carried out with fixing the xymovement (R Z scheme). In the other scheme, the whole structure is optimized with change in the supercell dimension (A XYZ scheme). 6.3.1 Geometry of the (310) Planes in SrTiO3 In order to generate the [001] (310) 5 symmetrical tilt GB (STGB), it is necessary to characterize the geometry of the (310) planes in SrTiO3. Figure 6 5. Schematic of the two types of (310) surface terminations. The top panel represents the SrTiO plane with a +4 charge, w hile the bottom panel shows the OO terminated surface with a net charge of 4. Figure 6 5 shows the possible (310) surface terminations in SrTiO3. The terminations in (310) follow the SrTiO/OO/SrTiO/.sequence. Considering formal charges of +2, +4 and 2 on Sr, Ti and O respectively, the surface terminations show a +4 charge per SrTiO termination and a 4 charge per OO termination. Therefore, care must be taken to generate a charge neutral Sr O Ti O a = 3.905 3.162a = 12.35 [ 1 3 0] [3 1 0] [0 0 1] PAGE 140 140 stoichiometric GB supercell structure which will be repeated perio dically. The GBs are constructed using GBstudio software. 6.3.2 Geometry of [001] (310) 5 STGB in SrTiO3 Figure 6 6 Schematic of the four possible twinned structures for (310) twist boundary in S rTiO3 [145] A) cation twin, B) anion twin, C) reverse cation twin, and D ) reverse anion twin. The solid line indicates the GB plane and the dotted lines indicate the repeat units for each structure. The arrows indicate the cations that are in close proximity of each other in each structure. A B C D GB GB PAGE 141 141 3.86 3.88 3.90 3.92 3.94 3.96 3.98 Experimental Grimes Thomas Sekiguchi CrawfordCS McCoyCS UDKumarCS AkhtarCS Akhtar TinteCS Tinte Lattice Parameter () Considering a positive negative stacking sequence and eliminating all improbable stackings (such as two like charged planes facing each other at the b oundary), four possible distinct structures can be obtained (Figure 6 6) [145] The structures illustrated in Figures 6 6(a) and 6 6(b) have either a cation rich or a n oxygenrich plane as the GB mirror plane. Hence, these are termed as cation twin (also mentioned as Ctwin) and anion twin (also mentioned as Atwin) respectively. Figures 6 6(c) and 6 6(d) are named as reverse cation twin (also mentioned as RevCtwin) and reverse anion twin (also mentioned as RevAtwin) respectively, indicating a reversal of the cation sublattice about the mirror plane (opposite type across the GB plane), while maintaining the oxygen sublattice twinning. 6.4 Characterization of Bulk SrTiO3 Before simulating the above structures to predict the stability of the [001] (310) 5 STGB, it is required to characterize bulk SrTiO3 with the interatomic interactions. Eight potential interactions [74, 146151] are considered to simulate SrTiO3. Figure 6 7 represents the calculated lattice parameter of SrTiO3 using for the available potential interactions. Figure 6 7. Prediction of cubic lattice parameter of SrTiO3 using various interatomic interactions. The names following CS indicate the core shell description of the ions in the models. PAGE 142 142 From Figure 6 7, it is cle ar that the core shell interaction given by Udayakumar [151] greatly over estimates the lattice parameter. Similar is the case for Crawfo rd [147] which is also another core shell model. Even though the interaction description of Sekiguchi [149] seems to over predict the lattice parameter, this rigid ion model description is more reliable in describing surface phenomena in SrTiO3 than all other potentials [152] Further analysis revealed that antiferrodistortive (AFD) structure is predicted to be more favorable for Akhtar [146] and Tinte [74] potentials in bulk SrTiO3. 6.4.1 Antiferrodistortion in SrTiO3 It is well known that SrTiO3 tranforms from paraelectric to an AFD trasition at ~105 K. Thus, the first task was to characterize the stability of bulk AFD phase with respect to the paraelectric phase. In order to characterize AFD, and also to check the effect of core shell on the phase transformation, all the calculations are performed with core shell and reduced rigid ion model for both the potentials. The core shell calcula tions are denoted by CS after the potential name. The AFD in SrTiO3 predicted by both Akhtar and Tinte potentials are characterized with respect to the energy grain for the transformation of the paraelectric to AFD phase (Figure 6 8). The energy gain pre dicted by Tinte potentials are much closer to the first principles calculations reported in the literature. However, the Akhtar potential shows almost 6~7 times larger energy grain compared to the LDA calculation [153] In general, the core shell model shows a larger grain in energy compared to the corresponding rigid ion m odel. The change in the lattice parameter (both a and c ) are characterized and compared with experiment and first principles calculations. Even though the first principles calculations predict a large variation in the a lattice constant (Figure 6 9(a)), t he interatomic potentials did a nice job in getting the values close to the experiment. Similar level of agreement is observed for the c/a PAGE 143 143 0.12 0.10 0.08 0.06 0.04 0.02 0.00 LDA PBE PBEsol HSE B3PW Akhtar AkhtarCS Tinte TinteCS Energy Gain (EAFDE) (eV) ratio (Figure 6 9(b)). Comparison of the shell model to the rigid ion description shows a lower value of predicted a lattice constant, thereby a larger c / a ratio for the shell model. Figure 6 8. Comparison of energy grain due to the formation of the AFD phase in SrTiO3 for Akhtar and Tinte potentials with first principle calculations [153] Figure 6 9. Comparison of A) a lattice parameter and B ) c/a ratio for the AFD p hase in SrTiO3 [153] Figure 6 10 shows the relative rotation angle c omparison for the AFD s tructure. Akhtar potential over predicts the rotation, while Tinte potential predicts similar angle predicted by 3.80 3.82 3.84 3.86 3.88 3.90 3.92 3.94 EXP at 65 K EXP at 78 K LDA PBE PBEsol HSE B3PW Akhtar AkhtarCS Tinte TinteCS Lattice Parameter (a in ) 0.98 0.99 0.99 1.00 1.00 1.01 1.01 EXP at 10 K EXP at 65 K LDA PBE PBEsol HSE B3PW Akhtar AkhtarCS Tinte TinteCS c/a ratio A B PAGE 144 144 0.00 2.00 4.00 6.00 8.00 10.00 12.00 14.00 EXP at 4.2 K EXP at 77 K LDA PBE PBEsol HSE B3PW Akhtar AkhtarCS Tinte TinteCS Rotation Angle (Theta) 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 EXP at 4.2 K EXP at 77 K LDA PBE PBEsol HSE B3PW Akhtar AkhtarCS Tinte TinteCS AFD Order Parameter q (x1011 m) some of the first principles calculations. However, the potentials overestimate the distortion compared to the experimen t. There is no trend between the coreshell and rigid ion description for the rotation angle. It is important to point out that even though there is distortion in the octahedra, there is no polarization induced in the system due to the AFD. Fig ure 6 10. Comparison of angle of rotation of the octahedral during AFD in SrTiO3 [153] Figure 6 11. AFD order parameter comparison between experiment, first principles calculations [153] and interatomic potential predictions in SrTiO3. Both first principles and interatomic potentials overestimate the order parameter. PAGE 145 145 The order parameter for AFD is characterized by q = (q1, q2, 0), w hich is proportional to the rotation angle ( ) of the oxygen octahedral. The magnitudes of q1 = (a/2)tan( ), and q2 = (b/2)tan( ), where a and b are the lattice parameters of SrTiO3 [154] Since the a and b lattice parameters predicted by both the potentials are same, the AFD is characterized by q = (a/2)tan( ) where, a is the lattice parameter of the AFD SrTiO3 structure [155] Figure 6 11 compares the AFD order parameters obtained from experiment with the first principle s and interatomic potential calculations. Both the simulation methods overestimated the AFD order parameter. Therefore, from all the above potential descriptions only five rigid ion descriptions, namely Akhtar and Tinte for AFD structures and Grimes, Thom as and Sekiguchi are used for cubic SrTiO3 GB calculations. The last three potentials are selected due to their differences in the description of the charge states. For Grimes potential the charge on each ion is described by full charge model (Sr +2, Ti +4 and O 2). Both Sekiguchi and Thomas are partial charge models. For Sekiguchi (Sr +1.331, Ti +2.662 and O 1.331), thus each (001) layer in bulk SrTiO3 is charge neutral, whereas for Thomas (Sr +1.84, Ti +2.36 and O 1.40) each (001) layer in bulk SrTiO3 is charged. Additionally, the core shell description of the potentials makes it difficult to optimize the geometry of the GBs. Ravikumar et al. [145] used the rigid i on Akhtar potential to study GBs. 6.5 Characterization of Cubic SrTiO3 Grain Boundary For the [001] (310) 5 STGB, the GB energies are calculated using both R Z and A XYZ schemes, as discussed in Section 5.5 for domain walls. All the descriptions are dis cussed comparing the relaxed Egb predicted from the potentials along with polarization and structural analysis around the GB. All the results are discussed for cubic SrTiO3 and then for the AFD systems. PAGE 146 146 6.5.1 R Z Relaxation Scheme The basic idea behind thi s relaxation scheme was to maintain the GB in plane lattice parameter and allow the relaxation only normal to the GB. The effect of system size is simulated to find out the minimum system size necessary to avoid the interaction between the GBs and achieve bulk behavior in the middle of each grain. The results on system size from the R Z scheme are used as a guideline for the A XYZ relaxation scheme. Figure 6 12. Egb for different twinned STGBs w ith respect to system size for A) Grimes, B) S ekiguchi, and C ) Thomas potentials in the R Z scheme. All the Egb results calculated with R Z scheme show a consistent trend for all the potentials. The Egb values decrease with increase in system size. The energy change is appreciable (~ 10%) for the syst ems with GB distance between ~12 and ~24 and shows a 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 0 10 20 30 40 50 60 70Distance between GBs (in )Relaxed GBE (mJ/m2) CTwin ATwin RevCTwin RevATwin 0 1000 2000 3000 4000 5000 6000 7000 0 10 20 30 40 50 60 70 Distance between GBs (in ) Relaxed GBE (mJ/m2) CTwin ATwin RevCTwin RevATwin 0 1000 2000 3000 4000 5000 6000 7000 0 10 20 30 40 50 60 70 Distance between GBs (in )Relaxed GBE (mJ/m2) CTwin ATwin RevCTwin RevATwin Grimes Sekiguchi Thomas A B C PAGE 147 147 0 500 1000 1500 2000 2500 3000 3500 4000 Ravikumar Grimes Thomas Sekiguchi Relaxed GBE (mJ/m2) Ctwin Atwin RevCtwin RevAtwin ~1% deviation for the systems with GB distance between ~50 to ~60 with all the potentials and structures (Figure 6 12). The Egb stability is reverse anion twin > cation twin > reverse cation twin > anion twin. The reverse anion twin is still the favorable Egb and anion twin is the least stable structure with the R Z scheme. 6.5.2 A XYZ Relaxation Scheme In this scheme all the atoms are allowed to fully optimize along with the ability to change t he supercell shape and volume. Most of the calculations in the literature allow a complete relaxation of the atoms while maintaining a constant volume of the supercell. Figure 6 13. Relaxed Egb predicted from four different potential descriptio ns for a fully dense (310) STGB in SrTiO3. The literature values are obtained from Ravikumar et al. [145] All the potentials predict the reverse twins to be more fav orable than the regular twins. The relaxed Egb (Figure 6 13) predicted by all the three potentials are in the comparable range for all the GBs as mentioned by Ravikumar et al. [145] using Akhtar potential. These energy values are substantially lower (~3 4 times) than the calculations from the R Z scheme, which indicates that the relaxation of atoms in all the three directions leads to a more stable GB PAGE 148 148 configuration. Al l the potentials predicted the anion twin to be the least favorable, while the regular twins show an overall larger Egb compared to the reverse twins. All the three potentials predicted the reverse anion twin to be more favorable. The reason for this diff erence in the Egb of the regular and reverse twins can be determined from the structure. First the extreme structures are considered. In the case of anion twin, the STGB contains cations (both Sr and Ti shown in arrows in Figure 6 6(b)), of the same kind a t close proximity to each other. In the perfect structure the first nearest similar kind of cations sit at a unit cell (3.905 ) apart. However, the closest Sr Sr and Ti Ti distance in the anion twinned structure are 1.235 Therefore, the electrostatic r epulsion between these cations generates a much open GB structure with comparatively larger energy. Considering the reverse anion twinned structure, there are no similar cations (both Sr and Ti shown in arrows in Figure 6 6(d)), or oxygens in close contact in the nearest neighbor cell. The arrows in Figure 6 6(d) emphasize the closest Sr and Ti atoms around the GB. The average Sr Ti bond distance in the perfect structure is ~3.382 and the nearest Sr Ti bond length in the initial reverse anion twinned s tructure is ~ 2.310 This comparatively lesser deviation from the bulk environment maintains a dense boundary with lower energy for reverse anion twins. For both cation twin (Figure 6 6(a)) and reverse cation twin (Figure 6 6(c)), the GB structure contai ns oxygens very close to each other (1.235 compared to the bulk value of ~2.76 ). Along with the anions, cations of the similar kind are close to each other (2.47 shown in arrows in Figure) at the GB for cation twins. Therefore, the electrostatic rep ulsion from anions is strong for both types of cation twins, while the regular cation twin also experiences additional cation repulsion, which is much weaker than the anion twin. S ince the charge on oxygens is PAGE 149 149 either equal or less than the cations for all the potentials, the cation cation interaction for anion twin results in a relatively higher Egb than the rest of the twinned structures. 6 .5.3 Structural Analysis Figure 6 14. Relaxed GB structure s with Sekiguchi potential for A) c ation twin, B) anion twin, C) reverse cation tw in, and D ) reverse anion twin. Visual inspection of individual GBs after relaxation shows a larger open volume for the anion twin (Figure 6 14(b)) than the rest of the GBs. However, the lattice distortion is significantly more for the reverse cation twin (Figure 6 14(c)). Small amount of distortion is also noticeable for the reverse anion twin (Figure 6 14(d)). Surprisingly even with a much higher energy, the anion twinned GB does not show any bulk distortion. Further ana lysis is necessary to provide more information about the GB structures. Polarization analysis was carried out for all the cases using results obtained from Sekiguchi potential. In order to be consistent with the domain wall notations, polarizations perpendicular to the boundary are denoted as Pn, while those parallel to the boundary (the axis shown along the A B C D Sr Ti O PAGE 150 150 plane of the paper for Figure 6 14) is given by Pt. Polarization parallel to the uniaxial direction (in and out of the page) for Figure 6 14 is repres ented by Pz. The maximum polarization variation is observed at the boundaries, which is as expected. The cation twin the structure generated a Pn polarization (headtail configuration) throughout the structure, with ~0.8 C/cm2 in magnitude. The maximum va lue of Pn is observed to be ~ 1.5 C/cm2 near the GB. Anion twin structures show a very high Pn near the GB, while this contribution is negligible in the bulk of each grain. A maximum of ~ 3.6 C/cm2 polarization is observed near the grain boundaries, indi cating substantial displacement of atoms around the boundaries. This is consistent with the visual analysis presented in Figure 6 14(b). A maximum Pt of ~1.2 C/cm2 polarization was observed near the GBs. Figure 6 15 represents the Pn and Pt variation in t he GB bicrystal system. The polarization at each boundary is independent of the other boundary, representing a different local structure attained at the boundaries after relaxation for the anion twin. At each boundary, both the polarization components are symmetric. The Pz values are calculated to be zero for all the cases. The reverse cation twin resulted in a large polarization throughout the structure. The maximum value of Pn and Pt calculated are ~9.8 C/cm2 and 6.5 C/cm2 respectively. This observation agrees well with the distortion observed in Figure 6 14(c). The reverse anion twin did show the presence of a very small amount of polarization near the GBs. These polarization values decays towards the interior of each grain. The maximum value of Pn an d Pt calculated are ~0.9 C/cm2 and 0.9 C/cm2 respectively. The perpendicular component of the polarization ( Pn) is observed to be symmetric at each boundary in the system, while the Pt component represented inversion symmetry in terms of polarization. PAGE 151 151 5.0 4.0 3.0 2.0 1.0 0.0 1.0 2.0 3.0 4.0 5.0 0 40 80 120 160 Distance (in ) Polarization (C/cm2) Pt Pn 1.0 0.5 0.0 0.5 1.0 0 40 80 120 160 Distance (in ) Polarization (C/cm2) Pt Pn Figure 6 15. Cell bycell polarization variation for anion twin GB using Sekiguchi potential. Both the boundaries in the bicrystal behave differently around the GB. Figure 6 16. Cell bycell polarization variation for the reverse anion twin GB using Sekiguchi potential. Pn values for both the boundaries show similar polarization variation while Pt values show inversion in polarization. PAGE 152 152 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 Ravikumar Akhtar Tinte Relaxed GBE (mJ/m2) Atwin RevAtwin 6.6 Characterization of AFD SrTiO3 Grain Boundary Although Grimes and Thomas potential follow the simi lar trend as Sekiguchi potential, the analysis of the GBs predicted by Akhtar and Tinte potentials needed special attention. Analogous to the cubic SrTiO3 characterization, energetics, polarization and structural analysis is performed for the AFD systems. Only A XYZ relaxation scheme is considered for this study. 6.6 .1 Energetics of GBs The energies of the anion twinned structures are tabulated for both Akhtar and Tinte potentials. The energy comparisons show the reverse anion twinned structures to be more favorable (Figure 6 17). The Tinte potential predicted an unusually high Egb for anion twinned structures. The cation twinned structures were unstable with the AFD and did not converge. Figure 6 17. Relaxed Egb predicted for AFD SrTiO3 for a full y dense (310) STGB. The literature values are obtained from Ravikumar et al. [145] The reverse anion twinned structures are more favorable than the anion twins. 6.6 2 Polarization and Structure Analysis The polarization analysis was carried out for the reverse anion twinned structures. Figure 6 18 represents the polarization variation for the Akhtar potential. It is interesting to notice the PAGE 153 153 presence of the uniaxial c omponent of the polarization ( Pz) along with the in plane polarization contributions ( Pn and Pt). The polarization generated at each GB show no correlation for each individual component. Figure 6 18. Cell bycell polarization variation f or the reverse anion twin GB us ing Akhtar potential, A ) the in plane polarizations and B ) in uniaxial polarization ( Pz). The polarization values do not show any correlation at the GBs. The comparison of the magnitude of Pn and Pt with the cubic structure s revealed that the AFD distortion generated much larger in plane polarizations. The presence of Pz hinted towards 20 15 10 5 0 5 10 15 20 0 20 40 60 80 Distance (in ) Polarization (C/cm2) Pt Pn 6 4 2 0 2 4 6 0 20 40 60 80 Distance (in ) Polarization (C/cm2) Pz A B PAGE 154 154 possible distortion along the uniaxial direction of the lattice. The maximum magnitude of Pz is observed to be nearly four times lower than t he maximum for Pn and Pt. All the polarization contributions reach zero at the bulk of each grain. Thus, structural analysis is necessary to confirm the distortion of lattice in all the three direction to support the presence of Pz, Pn and Pt. Figure 6 19. Structure of a relaxed reversed anion twin in (310) STGB of SrTiO3 with Akhtar potentia l showing octahedral rotation, A ) <00 1> view showing the GB region, and B) view of one grain showing octahedral rotation on the Pz Pn plane, and C ) view of one grain showing octahedral rotation on the Pz Pt plane, where the GB plane is parallel to the page. Oxygen atoms are not included for clarity. GB indicates the location of the grain boundary in the structure. The evidence of octahedral ro tation about the uniaxial direction (Figure 6 19(a)) is shown, which is expected for the AFD. However, additional octahedral rotation is also observed about the Pt (Figure 6 19(b)) and Pn (Figure 6 19(c)) directions. The distortion of the octahedras contri butes to the polarizations along all the three directions. A GB P n P t P n P z GB P t P z C B PAGE 155 155 15 10 5 0 5 10 15 0 10 20 30 40 50 60 70 80 Distance (in ) Polarization (C/cm2) Pt Pn Following the same characterization for the Tinte potential, the polarization profile demonstrated the presence of only inplane contributions ( Pn and Pt). Unlike Akhtar potential observations, the Pt profile shows inversion around the center of each grain, while Pn profile is symmetric around each boundary. The maximum value of Pn is calculated to be ~11 C/cm2, whereas the Pn reaches a maximum of ~ 7 C/cm2 near the GB. Even though Pn values go to zero away from the GB, Pt values show a polarization fluctuation (<1 C/cm2) in the bulk of each grain. Figure 6 20. Cell bycell in plane polarization variation for the reverse anion twin GB using Tinte potential. The Pt profile shows inversi on around the center of each grain, while Pn profile is symmetric around each boundary. The absence of Pz indicated no octahedral distortion around each axes. Structural analysis of the octahedral orientation is shown in Figure 6 21. Octahedral rotation is demonstrated only along the Pz axis (Figure 6 21(a)), while small amount of rotation is observed near the GB for the Pt (Figure 6 21(b)) and Pn (Figure 6 21(c)) axes. PAGE 156 156 Figure 6 21. Structure of a relaxed reversed anion twin in (310) STGB of SrTiO3 with Tinte potential showing octahedral rotation, A ) <00 1> view showing the GB region, and B) view of one grain showing octahedral rotation on the Pz Pn plane located near the GB, and C ) view of one grain showing octahedral rotation on the Pz Pt plane only near the GB, where the GB plane is parallel to the page. Oxygen atoms are not included for clarity. GB indicates the location of the grain boundary in the structure. Since the bulk AFD phase is much better explained by the Tinte pot ential compared with the Akhtar potential, the results calculated through the Tinte potential has more credibility for the GB structure and energetics. Other independent measurements are necessary to provide confidence in the predicted results. 6.7 Study o f 5 Tilt Boundary in PbTiO3 After successfully characterizing the structure and energetics of the SrTiO3 [001] (310) 5 grain boundary, the same GB is studied for PT bicrystal system. The additional challenge is the A GB P n P t B P n P z GB P t P z C PAGE 157 157 change in structure from cubic to tetragonal and the addition of polarization. Therefore, it is necessary to establish the probable GB systems for PT bicrystal. Considering the (310) plane in PT the (310) planes in SrTiO3. Figure 6 5 represents the (310) plane of SrTiO3. The most important structural variation is PT is tetragonal. Therefore, the tetragonal axis can be along a; along [ 1 3 0], or along the [310] direction, which is going in and out of the paper. This leads to three distinct [001] (310) 5 grain boundaries. However, the primary polarization along the tetragonal axis can orient in +z (Up) or z (Dn). Therefore, the structure and polarization combination results in six distinct GBs for PT. The Up Up or Up Dn polarization conditions setup for PT grain boundaries is similar to the domain wall setup discussed in Figure 5 3. The Dn Dn polarization condition is the same as the Up Up condition and is represented as Up Up for the rest part of the discussion. In this setup the primary polarizatio n is pointing parallel to the GB (pointing in and out of the page). All the GB calculations for PT are performed with Potential 2. Both R Z and A XYZ schemes are utilized for optimization. 6.7 .1 R Z Relaxation Scheme In this method the atoms are allowed to relax perpendicular to the GB. Comparing this scheme to the domain wall simulations it is apparent that the R Z approach restricts any in plane polarization for domain walls, whereas in GB cases this approach can develop polarization perpendicular ( Pn) to the primary polarization ( Pz). Prior analyzing the polarizations, the energy of each twinned structures are calculated for Up Up and Up Dn polarization configuration. Figure 6 22 represents the energy comparison for all the structures. The reverse anion t winned structures are predicted as the lowest energy structure for both Up Up and Up Dn polarization configurations. For all the cases the UpUp structures show comparatively lower energy than the PAGE 158 158 0 1000 2000 3000 4000 5000 6000 UpUp UpDn Grain Boundary Energy (mJ/m2) Ctwin Atwin RevCtwin RevAtwin corresponding UpDn cases. Using the experience of domain w all, it is clear that the Up Dn polarization configuration generates a domain wall at the GB. Therefore, the Egb for Up Dn cases is predicted to be comparatively higher than the corresponding Up Up cases. Figure 6 22. Relaxed Egb for (310) STGBs in PT with R Z scheme using Potential 2. Figure 6 23. Cell bycell in plane polarization variation for the Up Up reverse anion twin in PT using Potential 2. Due to the R Z scheme, the contribution of Pt is zero, while Pn profile is symmetric aro und each boundary. 15 10 5 0 5 10 0 20 40 60 80 100 Distance (in ) Polarization (C/cm2) Pt Pn Up Up Reverse Anion Twin PAGE 159 159 The analysis of the polarization contributions along the three directions showed no change in primary polarization ( Pz). The relaxation perpendicular to the boundary resulted in Pn in the unit cells near the GB (Figure 6 23). The devel oped inplane polarization is symmetric around each GB. The maximum value of Pn (~ 6 C/cm2) is observed near the domain walls and the contribution goes to zero inside each grain. The direction of these polarizations is difficult to present since it fluctu ates between positive and negative values near the boundary, however in general the direction of the polarization is either headtail or tail head (inout or out in) at each boundary. Figure 6 24. Cell bycell in plane polarization variation for the Up Dn reverse anion twin in PT using Potential 2. The contribution of Pt is zero, while Pn is non zero for any unit cell within the simulation system. Similar polarization analysis on the Up Dn system resulted in a finite amount of Pn throughout the reverse anion twinned structure. This characteristic may indicate either the system size considered for the Up Dn system is not large enough to avoid the long range interaction developed due to the presence of a domain wall at the GB, or the characteristi c is physical. 15 10 5 0 5 10 15 0 20 40 60 80 100 Distance (in ) Polarization (C/cm2) Pt Pn Up Dn Reverse Anion Twin PAGE 160 160 0 500 1000 1500 2000 2500 3000 3500 4000 4500 Atwin RevAtwin Grain Boundary Energy (mJ/m2) RZ AXYZ 6.7 .2 A XYZ Relaxation Scheme Approaching along the same line as the domain wall cases, it was obvious to investigate the effect of complete atomic relaxation in all the three directions. In this scheme the volume of the simulation cell was also allowed to change. Even though the method was attempted for all the structures with both Up Up and Up Dn polarizations, only the anion twinned structures are successfully optimized. Figure 6 25. Comparison of Egb with R Z and A XYZ relaxat ion schemes for (310) STGBs in PT with Up Up polarization using Potential 2. The relaxation of atoms along all the three directions resulted in much lower Egb for both the anion twins with Up Up polarization in PT (Figure 6 25). Considering the lowest ene rgy reverse anion twinned structure, cell by cell polarization is characterized. The in plane polarizations are determined to be independent of each boundary and go to zero within each grain. The uniaxial polarization is observed to be reducing at the GBs However, the bulk polarization achieved in each grain is ~2.3% larger than Pz calculated for bulk PT without GBs. PAGE 161 161 A closer look at the lattice shows a ~0.45% increase in the lattice parameter along the uniaxial polarization direction ( c lattice parameter ). Figure 6 26. Cell bycell polarization variation for the reverse anion twin using Potential 2 w ith A XYZ optimization scheme, A ) the in plane polarizations, and B ) in uniaxial polarization ( Pz). The polarization values do not show any correlation at the GBs.. The analysis was extended for anion twins. The converged uniaxial polarization in the bulk of each grain is ~6% less than the bulk Pz in PT single crystals. Thus, the relaxation of atoms in 15 10 5 0 5 10 15 0 20 40 60 80 100 Distance (in ) Polarization (C/cm2) Pt Pn 0 10 20 30 40 50 60 70 80 0 20 40 60 80 100 Distance (in ) Polarization (C/cm2) Pz A B PAGE 162 162 all directions and optimization of the s imulation volume resulted in structures with different bulk Pz values than observed for PT single crystal. 6.8 Conclusions Effect of symmetrical [001] (310) 5 tilt grain boundaries in SrTiO3 and PT are examined using interatomic potentials. All the results indicate the presence of in plane polarization due to the presence of grain boundaries in cubic and AFD SrTiO3 and ferroelectric PT. The values of the max imum in plane polarization is observed to be nearly ~5 times larger than maximum polarization values reported for the domain wall cases in Chapter 5. 6. 9 Future Work Due to the lack of literature on PT GBs, this systematic study will help explore the eff ect of GBs on polarization. First principles calculations on the similar systems will provide more details on the results calculated with the interatomic potentials. Along with the Up Up and UpDn polarization conditions, the other four possible structural and polarization combination will provide a complete understanding of the energetically favorable structures. Since development of polarization is very sensitive to the type of input structure, as discussed for SrTiO3, it will be interesting to establish if there is any direct relationship with the magnitude of polarization variation around the boundary with the stability of it. The effect of charge n eutral defects, especially the S chottky defects and their distribution around the GBs will provide insight into the interaction of defects with ferroelectricity. PAGE 163 163 CHAPTER 7 PARALLELIZATION OF T HE SHELL MODEL MOLEC ULAR DYNAMICS CODE 7.1 Introduction The simulation results discussed for the thin film surf aces and domain walls in Chapter s 4 and 5 respectively we re for single crystalline PT structures. The effect of grain boundaries were considered in Chapter 6 for a PT bicrystal. The natural choice to see the effect of grain boundaries and domain walls on the ferroelectric properties will be to simulate polycryst alline PT samples (with both bulk and surface configurations). The simulation of polycrystalline PT system will be taking a step towards a direct comparison between the simulated and experimental structures. However, this task is challenging with the simul ation approach considered for Chapters 4, 5, and 6, which uses a serial code running on a single processor to calculate the desired properties. Even though it is possible to simulate the complex polycrystalline system with the existing simulation approach, some systems are just too large for one processor. I t will take a very long time to achieve an y useful results. Figure 7 1 illustrates the effect of change in the total number of atoms in the simulated system (bulk PT) on the central processing unit (CPU ) time for the existing serial code. It is important to note that the average materials properties are calculated with a total 100,000 1.000,000 MD steps depending on the system energy profile, (the time step for each calculation is ~0.04 fs). The 1000 M D steps are considered as a reference for the discussion of th e simulation time in this Chapter The calculations are performed for interatomic interactions with the shell model and with the rigid ion model (by turning off the shell model) for Potential 2. As expected, the shell model description of the interactions took longer (~2.5 times) per MD step compared to the corresponding rigid ion model (Figure 7 1). The total CPU time increases linearly with the increase in the total number of atoms in the simul ated system, exhibiting the PAGE 164 164 expected O(N) scaling, where N is the total number of atoms. A typical polycrystalline s imulation system contains up to 106+ atoms, which indicates the time to simulate the same number of MD steps will take ~50+ times than the s ystem with 1080 atoms. This would lead to prohibitively a long waiting time for results. Therefore, special simulation approaches are developed to reduce this waiting time to a manageable range. Parallel computing is one of the approach in which the total simulation task is divided among more than one C PU to reduce the total time that one must wait for results. The larger the number of CPUs for a given system size, the lesser the actual wait time to finish a simulation. Figure 7 1. The effect of the total number of atoms in the simulation system on the total time for 1000 MD runs for a bulk PT system using Potential 2 using serial code. The shell model results took ~2.5 times longer time than the respective rigid ion models. The results show a linear increase in the total time with the increase in the total number of atoms. Inset shows the same data on a linear scale. The variation in the slope for the linear plot suggests that the shell model code is not properly optimized. 10 100 1000 1000 10000 Number of Atoms Time for 1000 MD Steps (sec) Shell Model Rigid Ion Model 0 100 200 300 400 500 600 700 800 1000 1500 2000 2500 3000 3500 4000 4500 5000 Number of Atoms Time for 1000 MD Steps (sec) Serial Code PAGE 165 165 7. 2 Parallel Comput ing The basic idea of parallel computing is to divide the total simulation task among multiple CPUs (referred as nodes), where each of them performs approximately the same amount of work. This is achieved by the use of a massage passing interface (MPI). In this process each CPU involved in the calculation function as an individual node where the only information each node knows about the other nodes is through the passing of information between them. This distribution of a simulation task among multiple nu mber of CPUs reduces the workload in each node and thereby reducing the wait time for the user. However, as the total number of nodes is increased for a fixed system size, the communication between the nodes can result in an increase of the total CPU time. 7.2. 1. Link Cell Algorithm and Spatial Decomposition A neighbor list is a basic approach to reduce the computational load by saving the information of atoms that are interacting or may interact within a relatively short simulation time with any atom. Use of a neighbor list approach scales as O( N2), where N is the total number of atoms in the system. A link cell algorithm along with the neighbor list is used to scale the process down to an O( N) calculation. The basic idea of link cell is to divide the syst em into cubic cells of length a little larger than the potential cut off. Therefore, each cubic cell (3D) is surrounded by 26 cells, and any atom within a cell can interact with atoms within the cell and its neighboring cells only. Considering Figure 7 2(b ) for a 2D case, Node 7 can interact with its eight neighboring cells (Nodes 2, 3, 4, 6, 8, 10, 11 and 12). This process of identifying the location of any atom reduces to an O( N) calculation, thereby resulting in an O( N) process for the link cell algorith m. The communication between the nodes is defined by the types of spatial decomposition scheme used for the MPI. The basic idea of spatial decomposition is to divide the simulation PAGE 166 166 system into several blocks such that each block is assigned to a particular node. This approach can ideally reduce the cost of the calculation by an order of 1/ M where M is the total number of processors used for the decomposition. Thus, for a fixed system size, the calculation time would ideally be reduced by half when the tota l number of nodes is doubled. Figure 7 2 shows the schematic of a 1D and 2D spatial decomposition. Figure 7 2. Schematic representation of A) 1D spatial decomposition, and B ) 2D spatial decomposition. The arrows show one part of the passi ng step. The red box around Node 7 represents the buffer region for this node, while the red arrows indicate the communication of Node 7 in the three direction s ). B A 1 2 3 4 5 6 7 8 9 10 11 12 1 3 14 15 16 1 2 3 4 PAGE 167 167 The 1D or 2D decomposition is defined by the division of the simulation system in either one or two directions respectively. Care must be taken while passing information from one node to another. Since the interactions considered in this thesis are defined by pair potentials, the atoms at the edge of each decomposition nodes have to carry info rmation from its adjacent nodes. This is made possible by defining a buffer layer around each neighboring nodes. Also, during simulation there is a possibility that the evolution of the system will result in the migration of atoms from one node to another. Thus, it is necessary to pass all the information of the migrated atom (like position, velocity, acceleration, and other higher derivatives and its global index) to the neighboring node. Special care should be taken to make sure that the core and shell of each atom are moved together in these conditions for shell model descriptions. The communication between the nodes is critical, and avoiding unnecessary communications between the nodes can improve the efficiency of the parallelization. The passing of inf ormation is performed in a specific way. In every passing step, first the even number nodes pass the information to the odd number nodes, then the odd number pass the information to the even number nodes. Therefore, at each passing step, a particular node can either only send or only receive information. This implementation requires an even number of nodes for the spatial decomposition. In order to have the necessary information of a node completely shared with its neighbors, the information is sent in two directions for a 1D spatial decomposition. However, for a 2D spatial decomposition it is done in three directions. First the nodes send information to the left and right, then top and bottom, and finally along the diagonal direction to pass all the buffer region atoms information. 7.2.2 Dynamic Memory Allocation Most of the variables with atom specific information are defined as arrays. Fixing the size of the arrays most of the time uses more space than required for each variable. This can be PAGE 168 168 avoided by using a dynamic memory allocation, which enables the array sizes to adjust automatically according to the system size. This also allows to allocate and deallocate the arrays, thus using the memory efficiently. With an effort to simulate polycrystalline fe rroelectric material, the existing rigid ion model parallel code with 2D spatial decomposition is modified to include the core shell description for each atom. The major tasks include updating the position for each atom to include the core and shell; mass and charge of core and shell; longrange interaction of corecore, core shell; short range interaction of shell shell. 7.3 Energy Comparison After the implementation of the shell model, the code is tested for energies, forces and stresses for bulk PT st ructures. Two different energy characterizations are performed to make sure the code is predicting t he correct structure. Table 7 1 shows the total energy contributi on from the core and shell of each atom type compared between the serial code and the modif ied parallel code The results show an exact match in all the three atom types for both core and shell energy contributions. The energies for Pb shel and Ti shel vary at the fifth significant figure between the two codes. Individual longand short range c ontributions also matched with the serial code for each atom type, indicating proper implementation of the charge and shell model. Table 7 1 Comparison of energy/atom for core and shell of each atom type in PT between the existing serial code and modified parallel code Atom Type Serial Code Energy (eV/atom) Parallel Code Energy (eV/atom) Pb core 44.64325 44.64325 Pb shell 30.77267 30.77268 Ti core 157.42347 157.42347 Ti shell 129.97852 129.97853 O core 6.51691 6.51691 O shell 18.88514 18.88514 PAGE 169 169 In order to verify the correctness of the code, the energy contributions from each atom is calculated with the sh ell model turned off. Table 7 2 demonstrates the energy con tributions from each atom type matched well for the rigid ion model. The Pb core and Ti core atoms energy contributions differ at the fifth significant figure. This test confirms that the updated code is now useful in using for systems with pure rigid ion, pure shell model or mixed rigid ion and shell model interactomic interac tion descriptions. Similar levels of agreement are also achieved for the forces on each atom and the stresses in the system. Table 7 2 Comparison of energy/atom of each atom type in PT by turning off the shell contributions between the existing serial code and modified parallel code Atom Type Serial Code Energy (eV/atom) Parallel Code Energy (eV/atom) Pb core 14.25710 14.25718 Ti core 43.65772 43.65773 O core 12.39091 12.39091 7.4 Performance of Parallel Code In order to calculate the effect of the total number of processors on the efficiency of the simulation, a large bulk PT system (6 x 24 x 60 unit cell, 43200 atoms) is considered. All the calculations in this section are performed with both shell model and rigid ion model (switching off th e shell model) to check the relative time used with both the potential description for Potential 2. The calculations are performed with 1, 4, 8, 12, 16 and 20 processors. The 1 processor jobs are obtained from the serial code calculations. Figure 7 3 sho ws the 2D spatial decomposition of the 6 x 24 x 60 system for 4 (F igure 7 3(a)) and 20 (Figure 7 3 (b)) processors. All the decomposition was performed along the yz plane. The decomposition resulted in a 2 x 2 and 2 x 10 division of the total system for the 4 and PAGE 170 170 20 processors respectively. The unit cells in each processor were carefully selected to take care of the any cut off issues. Figure 7 3 2D spatial decomposition of the 6 x 24 x 60 system with A ) 4 (2 x 2) and B ) 20 (2 x 10) processors. Due to the 2D spatial decomposition multiple s of 4 processors are required for simulations. For a fixed number of atoms in the total simulation system, the speed increased ~2.9 times for the core shell model, and ~2.4 times for the rigi d ion model between one and four processors. However, since both the codes are different, the detailed quantitative comparison is not very useful. Therefore, simulation time calculations are carried out between 4 20 processors for a fixed number of atoms i n the simulation system. Figure 7 4 represents the time taken for 1000 MD steps for the same system size with variable number of processors with and without shell A B PAGE 171 171 0 500 1000 1500 2000 2500 3000 0 4 8 12 16 20 24 Number of Processors Time for 1000 MD Steps (sec) Shell Model Rigid Ion Model 0.0 0.2 0.4 0.6 0.8 1.0 1.2 0 4 8 12 16 20 24 Number of Processors Normalized (Time/Node/Step/Atom) model using Potential 2. For shell model descripti on, the speed is increased by ~2.6 times fo r 20 processors in comparison to the four processors, whi le it is increased by ~1.8 times for the rigid ion model. In general, the shell model description is ~1.4 times slower than the rigidion model for the 20 processor calculations. Figur e 7 4 Total simulation CPU time for running 1000 MD steps with number of processors for a 2D spatial decomposition using 6 x 24 x 60 system size (43200 atoms) at 100 K with 0.04 fs time step. Both shell model and rigidion model show a substantial reducti on in time with increase in the number of processors. Inset shows the normalized time (Time/Node/Step/Atom) for both shell model and rigidion model with respect to the four processor time. Similarly, another set of calculations are performed in which the total number of atoms in each node is kept constant (2160 atoms). The idea is to check the effect of communication between nodes and how it affects the simulation time. The calculations are performed with 2160 atoms with one processor, 8640 atoms with 4 pr ocessors, 17280 atoms with 8 processors, 25920 PAGE 172 172 0 100 200 300 400 500 600 700 800 900 1000 0 4 8 12 16 20 24 Number of Processors Time for 1000 MD Steps (sec) Shell Model Rigid Ion Model atoms using 12 processors, 34560 atoms with 16 processors and 43200 atoms with 20 processors. Using Potential 2 for the same set of 1000 MD steps with 0.04 fs, the results predict a linear increase in the simu lation time with increase in the number of processors (Figure 7 5 ). This is expected, since the increase in the number of processors increases the total communication necessary in between the neighboring nodes. Figure 7 5 Total simulation C PU time for running 1000 MD steps with number of processors for a 2D spatial decomposition with 2160 atoms in each processor. Both shell model and rigid ion model show a linear increase in time with increase in processors. In addition to the above observa tions, it is interesting to point out that the 16 processor calculation is special. This is because there are two different possible ways the 2D spatial decomposition can be achie ved, 2 x 8 and 4 x 4. Figure 7 6 shows the schematic of both the decompositions. Even though the total number of interfaces per node is the same in both the cases, the 2 x 8 decomposition shows ~4.5% less time for the 1000 MD steps simulated in PAGE 173 173 comparison to the 4 x 4 decomposition. This can be explained by the difference in the co mmunication time. Figure 7 6 2D spatial decomposition of the 6 x 24 x 48 system with 16 processors A) 2 x 8, and B) 4 x 4 processors along yand z direction respectively. Further analysis on even larger number of processors will be helpful to benchmark the point beyond which the code is not efficient for the PT system. The effect of temperature on the phase transformation is also very important to substantiate the shell model implementation which is not considered here. With a suitable P T polycrystal structure, this code will be very helpful to analyze the effect of grain boundaries and surfaces along with domains for PT and related systems. B A PAGE 174 174 CHAPTER 8 GENERAL CONCLUSIONS Atomistic simulation methods have been successfully used to establ ish the effect of surfaces, domain walls and grain boundaries on the ferroelectric properties in PbTiO3. These effects are characterized based on energies, cell by cell polarization variations and structural analyses. 8.1. Surfaces In this thesis, it was determined that under open circuit electrical boundary conditions, t he (100) surface characterizations in PbTiO3 determined the Out polarization surfaces are more stable than the In polarization surfaces for both PbO and TiO2terminations. For all cases, the polarization was reduced at or near the surface layers and attains the bulk value away from the surface. Comparing the deviation in polarization from the bulk, the Out terminated surfaces showed a larger reduction in polarization compared to the In ter minations. The interplanar spacing was reduced in the first layer as expected for all the terminations with atoms moving towards the bulk. For free standing films, polarization parallel to the surface was predicted to be the most stable surface polarizatio n combination. The current status of atomic level description of substrate effect is limited to only very thin films (less than 10 unit cells) [35, 36, 101, 103] using first principles and atomic level simulations. The fre e standing films considered in this study are a very special case. In most physical situations, th in films are grown on a substrate. The presence of a substrate substantially modifies both the electrostatics of the system and the strain state, and consequently the stability of the system. Fong et al. [32] and Streiffer et al. [118] experimentally investigated the presence of stripe domains in ultrathin PT thin films grown on SrTiO3 (insulating) s ubstrates. A minimum thickness of three unit cells of PT thin film was established to exhibit any polarization. The PAGE 175 175 effect of SrRuO3 substrate on polarization of PbZr0.2Ti0.8O3 (PZT) film was calculated from experimentally observed atomic positions [156, 157] The profile showed a reduction in polarization due to the presence of the surface, similar to Figure 48. In addition to the surface effect, polarization was reduced due to the presence of the interface between SrRuO3:PZT. With the knowledge acquired from the surface studies, the effect of substrates can be investigated systematically with film thickness. The important aspects to address will be the influence of the interface, the effect of different types of substrates (metallic or in sulating), and the effects of strain due to the planar mismatch between the substrate and film on polarization. In addition, experimental results illustrated a decrease in ferroelectric transition temperature ( Tc) with decreasing film thickness for PT thin films on a SrTiO3 substrate [32] The LGD formalism underestimates the Tc for thinner films; for example, for a film of ~ 2.5 nm thickness, LGD theory predicts the Tc ~ 330 K compared to ~ 800 K measured experimentally [158, 159] There is currently very little under standing on the mechanism(s) responsible for this enhancement of ferroelectricity for very thin films (< 10 nm). Due to the inability of incorporating the effect of temperature in first principles calculations, there is little understanding of the influenc e of temperature on ultrathin ferroelectric films. Insight into this issue can be obtained with the molecular dynamics approach using interatomic potentials. 8.2. Domain Wall s In this thesis, the structure and energetics of 180 domain walls were analyzed with electronic level simulations for different wall terminations. The analysis of the energetics of the domain walls showed that Pb centered walls are more favorable for (100) terminations, in agreement with the literature. The study was extended for (110 ) domain walls, and the OO centered walls were predicted to be more stable. Also, the energies of (100) Pb centered (124 PAGE 176 176 mJ/m2) and (110) OO centered domain walls (121 mJ/m2) were found to be comparable to each other. These similar energies indicate that the energetically favorable path for the wall motion might be a combination of (100) Pbcentered and (110) OO centered domain walls for 180 domains. This suggests that the domain wall may prefer to avoid the high energy (100) Ti centers and (110) PbTiO c enters, and nanofacet during the wall motion. However, it is important to note that the domain wall areas for (100) and (110) terminations are different. Thus, the energy gain due to faceting along a (100) (110) combination may or may not necessarily high er than the (100) terminated domain walls. The nucleation and growth model proposed by Shin et al. [123] presents a similar combination of (100) (110) domain walls in PT. It will be definitely interesting to verify this mechanism of predicted wall motion using atomistic simulations. In addition to the energy comparisons, the structural aspects of the influence of domain wall were addressed by comparing the influence o f different atomic relaxation schemes in this work. In the literature, only tetragonal z direction relaxation of the atoms was allowed during structural optimization. In this study the atoms were allowed to relax in all the three directions. These additional degrees of freedom resulted in a key new result: the optimized domain wall structures possessed additional polarization components normal to the domain wall ( Pn) along with the primary ferroelectric polarization ( Pz) for both (100) and (110) domain wall s. The maximum magnitude of Pn was calculated to be (~2.4 C/cm2) within one unit cell around the domain wall and goes to zero away from the domain wall. The direction of these in plane polarizations points away from the domain wall. However, no in plane t ransverse polarization (Pt) is observed for any of the domain walls considered in this study. PAGE 177 177 The presence of Pn along with Pz required a simultaneous fitting to calculate the domain wall width. Two different analyses performed by expressing the polarizat ions as a polynomial expansion in the odd powers of the primary order parameter Pz yielded domain walls differing in width by a factor of two. Similar effect is observed in recent theoretical calculation by Lee et al. [134] for Y walls in LiNbO3. None of the calculations performed for domain walls considered defects in the simulation systems. However, defects play an important role in bulk ferroelectric materials. The presence of oxygen vacan cies and the diffusion of ions have reported as the possible origins for polarization fatigue. First principles calculations are widely used to address the influence of defects in bulk PT [160164] While most of the studies discuss the effect of these defects in single crystal bulk polarization, there are a few studies on the interaction of these defects at the domain wall [161] with the formation energies of oxygen vacancies reported to be much lower at the domain wall than in the bulk. It will be interesting to characterize the effect of Pb and O vacancies on the polarization variation across the domai n wall. The symmetry of this domain wall is maintained around the Pb atoms present at the 4mm sites. Therefore, the effect of breaking the symmetry by Pbvacancies and O vacancies on the polarization could provide further structure polarization interaction around the domain wall. Analysis of this issue could be extended to the (110) OO centered domain wall, where the effect of O vacancies will be characterized. 8.3. Grain Boundar ies The influence of 5 (310) symmetrical tilt boundaries on polarization wer e characterized for SrTiO3 and PT using interatomic potentials. For both systems, the reverse anion twinned PAGE 178 178 structures was calculated to be energetically favorable. The structural investigation was extended to calculate polarization in each system. For cub ic SrTiO3, a very small amount of polarization parallel ( Pt ~ 0.9 C/cm2) and ( Pn~ 0.9 C/cm2) to the GB were determined for the reverse anion twin. The Pt were found to orient along the same direction across each domain wall, while maintaining inversion s ymmetry for the whole bicrystal system around the center of each grain. The Pn oriented in a head tail fashion at the GB, showing similar polarization behavior at each GB. For AFD SrTiO3, two different interatomic descriptions resulted in different polariz ation characteristics. On one hand, the Akhtar potential predicted the presence of Pz along with Pn and Pt. A magnitude of ~20 C/cm2 polarization was predicted for both Pn and Pt, while Pz values were approximately one fourth of this value. On the other hand, the Tinte potential predicted only the presence of Pn and Pt, with Pz being zero magnitude of ~10 and 7 C/cm2 polarization is predicted for both Pn and Pt respectively. The polarization characteristics were observed to be the same as the cubic SrTiO3 GBs. Thus while the empirical potentials strongly indicate that the GB does have an effect on the polarization, a DFT study will be required to provide a definitive characterization of the effect, Similar to SrTiO3 GB, for the PbTiO3 GB with Up Up polari zation configuration polarizations of magnitude ~11 and 4 C/cm2 were developed for Pn and Pt respectively. The primary polarization ( Pz) was reduced at the GB, but was difficult to quantify. In particular, the relaxation of atoms in all the three directio ns (A XYZ scheme) resulted in ~ 0.45 % increase in lattice parameter along the primary polarization direction, which added ~2.3% to the bulk polarization Pz. However, optimization of atoms only perpendicular to the GB (R Z scheme), developed a Pn value of ~ 6 C/cm2. For the Up Dn polarization condition across the GB, a non PAGE 179 179 zero (~ 3 C/cm2) normal polarization ( Pn) was developed through the entire GB system. Additional calculations on all other primary polarization orientation are necessary for a complete understanding of the 5 (310) symmetrical tilt boundaries. The calculations performed in this work focused on fully dense GB structures. Similar to the issue of defects at the domain wall, previous work on GBs addressed the influence of Schottky defects i n perovskites [144, 145] Since this study is the first for CSL boundaries in PT, it will also be necessary to characterize the structure and energetics of defects at the GBs. In addition, this study can be extended to include other CSL boundaries for both PT and SrTiO3. The effect of the AFD order parameter in SrTiO3 on the structure of different CSL boundaries will also provide more structural understanding of these special boundar ies. The results discussed in this dissertation provide us with great deal of confidence in the ability to use atomistic simulations to address the effects on polarization in PT. In order to simulate complex and large systems, the current first generation parallel shell model code will be very useful. The optimization of this code will help characterize complex ferroelectric systems efficiently. Finally, the energetics of the surfaces, domains and GBs calculated from the first principles calculations in th is work can be utilized as a database to generate the next generation interatomic interactions to describe the PT system. The charged optimized many body (COMB) [165] formalism developed by the Computational Materials Science Focus Group (CMSFG) at the University of Florida (UF) will allow the incorporation of metal oxide semiconductor systems, thereby allowing the predict ion of characteristics of materials structures such as life size ferroelectric nanodevices. PAGE 180 180 LIST OF REFERENCES [1] J. Valasek, "Piezoelectric and allied phenomena in Rochelle salt," Phys. Rev. vol. 15, pp. 537538, 1920. [2] J. Va lasek, "Piezo electric and allied phenomena in Rochelle salt," Phys. Rev. vol. 17, pp. 475481, 1921. [3] D. Dimos and C. H. Mueller, "Perovskite thin films for high frequency capacitor applications," Annu. Rev. Mater. Sci. vol. 28, pp. 397419, 1998. [4] D. L. Polla and L. F. Francis, "Processing and characterization of piezoelectric materials and integration into microelectromechanical systems," Annu. Rev. Mater. Sci. vol. 28, pp. 563597, 1998. [5] J. F. Scott, "High dielectric constant thin films for dynamic random access memories (DRAM)," Annu. Rev. Mater. Sci. vol. 28, pp. 79 100, 1998. [6] O. Auciello, J. F. Scott, and R. Ramesh, "The physics of ferroelectric memories," Phys. Today vol. 51, pp. 22 27, 1998. [7] M. E. Lines and A. M. Glass, P rinciples and Applications of Ferroelectrics and Related Materials Oxford: Clarendon Press, 1979. [8] J. F. Scott, Ferroelectric Memories : Springer, 2000. [9] K. Uchino, Ferroelectric Devices New York: Marcel Dekker, Inc., 2000. [10] S. O. Kasap, Prin ciples of electronic materials and devices 3rd ed: McGraw Hill International, 2006. [11] C. Kittel, Introduction to Solid State Physics 6th ed. New York: John Wiley & Sons, 1986. [12] M. T. Dove, Introduction to Lattice Dynamics Great Britain: Cambrid ge University Press, 1993. [13] T. Kolodiaynyi, "Phonon Dispersion," http://www.chembio.uoguelph.ca/educmat/chm729/Phonons/cont.htm April 16, 1997. [14] T. Shimada, K. Wakah ara, Y. Umeno, and T. Kitamura, "Shell model potential for PbTiO3 and its applicability to surfaces and domain walls," J. Phys.: Condens. Matter vol. 20, pp. 325225, 2008. [15] S. Tinte, M. G. Stachiotti, S. R. Phillpot, M. Sepliarsky, D. Wolf, and R. L. Migoni, "Ferroelectric properties of BaxSr1xTiO3 solid solutions obtained by molecular dynamics simulation," J. Phys.: Condens. Matter vol. 16, pp. 34953506, 2004. PAGE 181 181 [16] C. A. Randall, R. E. Newnham, and L. E. Cross, "History of the first ferroelectric oxide, BaTiO3," http://ceramics.org/wp content/uploads/2009/05/first_ferroelectric_oxide_ba_tio3.pdf 2004. [17] A. J. Moulson and J. M. Herbert, Electr oceramics: Materials, Properties, Applications New York: J. Wiley, 2003. [18] Y. Xu, Ferroelectric materials and their applications Amsterdam: North Holland, 1991. [19] G. Shirane, R. Pepinsky, and B. C. Frazer, "X Ray And Neutron Diffraction Study Of Ferroelectric PbTiO3," Acta Crystallogr. vol. 9, pp. 131140, 1956. [20] A. S. Sidorkin, Domain Structure in Ferroelectrics and Related Materials Cambridge, UK: Cambridge International Science, 2006. [21] R. J. D. Tilley, Understanding solids: the scie nce of materials West Sussex, England: John Wiely & Sons Ltd, 2004. [22] M. J. Haun, E. Furman, S. J. Jang, H. A. McKinstry, and L. E. Cross, "Thermodynamic Theory Of PbTiO3," J. Appl. Phys. vol. 62, pp. 33313338, 1987. [23] C. L. Wang and S. R. P. Sm ith, "Landau Theory Of The Size Driven Phase Transition In Ferroelectrics," J. Phys.: Condens. Matter vol. 7, pp. 71637171, 1995. [24] M. D. Glinchuk, E. A. Eliseev, V. A. Stephanovich, and R. Farhi, "Ferroelectric thin film properties: Depolarization f ield and renormalization of a "bulk" free energy coefficients," J. Appl. Phys. vol. 93, pp. 11501159, 2003. [25] B. Jiang and L. A. Bursill, "Phenomenological theory of size effects in ultrafine ferroelectric particles of lead titanate," Phys. Rev. B v ol. 60, pp. 99789982, 1999. [26] J. W. Hong and D. N. Fang, "Size dependent ferroelectric behaviors of BaTiO3 nanowires," Appl. Phys. Lett. vol. 92, pp. 012906, 2008. [27] H. Ishiwara, M. Okuyama, and Y. Arimoto (Eds.), Ferroelectric Random Access Memo ries: Fundamentals and Applications Berlin: Springer, 2004. [28] K. M. Rabe, C. H. Ahn, and J. M. Triscone (Eds.), Physics of Ferroeletrics: A Modern Perspective. Berlin / Heidelberg: Springer, 2007. [29] "International technology roadmap for semiconduc tors (ITRS)," Semiconductor Industry Association, San Jose 2003. [30] Y. Arimoto and H. Ishiwara, "Current status of ferroelectric random access memory," Mater. Res. Soc. Bull. vol. 29, pp. 823 828, 2004. PAGE 182 182 [31] S. Hong (Ed.), Nanoscale Phenomena in Ferroelectric Thin Films : Kluwer Academic, 2003. [32] D. D. Fong, G. B. Stephenson, S. K. Streiffer, J. A. Eastman, O. Auciello, P. H. Fuoss, and C. Thompson, "Ferroelectricity in ultrathin perovskite films," Science, vol. 304, pp. 16501653, 2004. [33] C. Li chtensteiger, J. M. Triscone, J. Junquera, and P. Ghosez, "Ferroelectricity and tetragonality in ultrathin PbTiO3 films," Phys. Rev. Lett. vol. 94, pp. 047603, 2005. [34] P. Ghosez and K. M. Rabe, "Microscopic model of ferroelectricity in stress free PbT iO3 ultrathin films," Appl. Phys. Lett. vol. 76, pp. 27672769, 2000. [35] J. Junquera and P. Ghosez, "Critical thickness for ferroelectricity in perovskite ultrathin films," Nature vol. 422, pp. 506509, 2003. [36] M. Sepliarsky, M. G. Stachiotti, and R. L. Migoni, "Surface and substrate effects on the ferroelectric properties of PbTiO3 ultrathin films," Ferroelectrics vol. 335, pp. 3 12, 2006. [37] R. H. Mitchell, Perovskites: Modern and Ancient Thunder Bay, Ontario: Almaz Press, 2002. [38] W. Won g Ng, T. Holesinger, G. Riley, and R. Guo, Perovskite Oxides for Electronic, Energy Conversion, and Energy Efficiency Applications Westerville, Ohio: The American Ceramic Society, 2000. [39] P. Goudochnikov and A. J. Bell, "Correlations between transition temperature, tolerance factor and cohesive energy in 2+: 4+perovskites," J. Phys.: Condens. Matter vol. 19, pp. 176201, 2007. [40] M. R. Levy, R. W. Grimes, and K. E. Sickafus, "Disorder processes in A3+B3+O3 compounds: implications for radiation toler ance," Philos. Mag. vol. 84, pp. 533545, 2004. [41] M. D. Graef and M. E. McHenry, Structure of Materials: An Introduction to Crystallography, Diffraction and Symmetry Cambridge: University Press, 2007. [42] S. Piskunov, E. Heifets, R. I. Eglitis, and G. Borstel, "Bulk properties and electronic structure of SrTiO3, BaTiO3, PbTiO3 perovskites: an ab initio HF/DFT study," Comput. Mater. Sci. vol. 29, pp. 165178, 2004. [43] S. Piskunov, E. A. Kotomin, E. Heifets, J. Maier, R. I. Eglitis, and G. Borstel "Hybrid DFT calculations of the atomic and electronic structure for ABO3 perovskite (001) surfaces," Surf. Sci. vol. 575, pp. 75 88, 2005. PAGE 183 183 [44] M. Dawber, K. M. Rabe, and J. F. Scott, "Physics of thin film ferroelectric oxides," Reviews of Modern Physi cs vol. 77, pp. 10831130, 2005. [45] R. E. Cohen, and Krakauer, H., "Lattice dynamics and origin of ferroelectricity in BaTiO3: Linearized augmentedplane wave total energy calculations," Phys. Rev. B vol. 42, pp. 64166423, 1990. [46] R. E. Cohen, an d Krakauer, H., "Electronic structure studies of the differences in ferroelectric behavior of BaTiO3 and PbTiO3," Ferroelectrics vol. 136, pp. 6583, 1992. [47] R. E. Cohen, "Theory of ferroelectrics: a vision for the next decade and beyond," J. Phys. Ch em. Solids vol. 61, pp. 139146, 2000. [48] S. Tinte and M. G. Stachiotti, "Surface effects and ferroelectric phase transitions in BaTiO3 ultrathin films," Phys. Rev. B vol. 64, pp. 235403, 2001. [49] S. Tinte, M. G. Stachiotti, M. Sepliarsky, R. L. Mi goni, and C. O. Rodriguez, "Atomistic modelling of BaTiO3 based on first principles calculations," J. Phys.: Condens. Matter vol. 11, pp. 96799690, 1999. [50] M. Sepliarsky, Z. Wu, A. Asthagiri, and R. E. Cohen, "Atomistic model potential for PbTiO3 and PMN by fitting first principles results," Ferroelectrics vol. 301, pp. 55 59, 2004. [51] M. Sepliarsky, S. R. Phillpot, D. Wolf, M. G. Stachiotti, and R. L. Migoni, "Atomic level simulation of ferroelectricity in perovskite solid solutions," Appl. Phys. Lett. vol. 76, pp. 39863988, 2000. [52] M. Sepliarsky, S. R. Phillpot, D. Wolf, M. G. Stachiotti, and R. L. Migoni, "Ferroelectric properties of KNbO3/KTaO3 superlattices by atomiclevel simulation," J. Appl. Phys. vol. 90, pp. 45094519, 2001. [53] M. Sepliarsky, M. G. Stachiotti, and R. L. Migoni, "Interface effects in ferroelectric PbTiO3 ultrathin films on a paraelectric substrate," Phys. Rev. Lett. vol. 96, pp. 137603, 2006. [54] M. Sepliarsky, M. G. Stachiotti, and R. L. Migoni, "Surface recon struction and ferroelectricity in PbTiO3 thin films," Phys. Rev. B vol. 72, pp. 014110, 2005. [55] M. G. Stachiotti, "Ferroelectricity in BaTiO3 nanoscopic structures," Appl. Phys. Lett. vol. 84, pp. 251253, 2004. [56] M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids Oxford: Oxford Science, 2004. PAGE 184 184 [57] R. M. Martin, Electronic Structure: Basic Theory and Practical Methods New York: Cambridge University Press, 2004. [58] D. S. Scholl and J. A. Steckel, Density Functional Theory: A Prac tical Introduction Hoboken, New Jersey: John Wiley & Sons, Inc., 2009. [59] P. Hohenberg and W. Kohn, "Inhomogeneous Electron Gas," Phys. Rev. vol. 136, pp. B864B871, 1964. [60] W. Kohn and L. J. Sham, "Self Consistent Equations Including Exchange And Correlation Effects," Phys. Rev. vol. 140, pp. A1133A1138, 1965. [61] J. P. Perdew, "Accurate Density Functional for the Energy: Real Space Cutoff of the Gradient Expansion for the Exchange Hole," Phys. Rev. Lett. vol. 55, pp. 16651668, 1985. [62] N W. Ashcroft and N. D. Mermin, Solid State Physics New York: Holt, Rinehart and Winston, 1976. [63] H. J. Monkhorst and J. D. Pack, "Special points for Brillouin zone integrations," Phys. Rev. B vol. 13, pp. 51885192, 1976. [64] D. Vanderbilt, "Soft self consistent pseudopotentials in a generalized eigenvalue formalism," Phys. Rev. B vol. 41, pp. 78927895, 1990. [65] P. E. Blochl, "Projector Augmented Wave Method," Phys. Rev. B vol. 50, pp. 1795317979, 1994. [66] G. Kresse and J. Joubert, "From Ultrasoft Pseudopotentials to the Projector AugmentedWave Method," Phys. Rev. B vol. 59, pp. 17581775, 1999. [67] G. Henkelman, A. Arnaldsson, and H. Jonsson, "A fast and robust algorithm for Bader decomposition of charge density," Comput. Mater. Sci. vol. 36, pp. 354360, 2006. [68] E. Sanville, S. D. Kenny, R. Smith, and G. Henkelman, "Improved gridbased algorithm for Bader charge allocation," Journal Of Computational Chemistry vol. 28, pp. 899908, 2007. [69] W. Tang, E. Sanville, and G. Henkelm an, "A grid based Bader analysis algorithm without lattice bias," J. Phys.: Condens. Matter vol. 21, pp. 084204, 2009. [70] D. Wolf, P. Keblinski, S. R. Phillpot, and J. Eggebrecht, "Exact method for the simulation of Coulombic systems by spherically truncated, pairwise r1 summation," J. Chem. Phys. vol. 110, pp. 82548282, 1999. PAGE 185 185 [71] D. Fincham, "Optimisation of the Ewald Sum for Large Systems," Molecular Simulation vol. 13, pp. 1 9, 1994. [72] B. G. Dick Jr. and A. W. Overhauser, "Theory of the Die lectric Constants of Alkali Halide Crystals," Phys. Rev. vol. 112, pp. 90103, 1958. [73] A. Asthagiri, Z. Wu, N. Choudhury, and R. E. Cohen, "Advances in first principles studies of transducer materials," Ferroelectrics vol. 333, pp. 6978, 2006. [74] M. Sepliarsky, A. Asthagiri, S. R. Phillpot, M. G. Stachiotti, and R. L. Migoni, "Atomic level simulation of ferroelectricity in oxide mateirals," Curr. Opin. Solid State Mater. Sci. vol. 9, pp. 107113, 2005. [75] L. Verlet, "Computer "Experiments" on Classical Fluids. I. Thermodynamical Properties of Lennard Jones Molecules," Phys. Rev. vol. 159, pp. 98103, 1967. [76] C. W. Gear, Numerical Initial Value Problems in Ordinary Differential Equations Englewood Cliffs, NJ: Prentice Hall, 1971. [77] D. G. Luenberger, Linear and Nonlinear Programming: Kluwer Academic Publishers, 2004. [78] H. C. Andersen, "Molecular dynamics simulations at constant pressuure and/or temperature," J. Chem. Phys. vol. 72, pp. 23842393, 1980. [79] M. Parrinello and A. Rah man, "Crystal Structure and Pair Potentials: A Molecular Dynamics Study," Phys. Rev. Lett. vol. 45, pp. 11961199, 1980. [80] M. Parrinello and A. Rahman, "Strain fluctuations and elastic constants," J. Chem. Phys., vol. 76, pp. 26622666, 1982. [81] H. J. C. Berendsen, J. P. M. Postma, W. F. v. Gunsteren, A. DiNola, and J. R. Haak, "Molecular dynamics with coupling to an external bath," J. Chem. Phys. vol. 81, pp. 36843690, 1984. [82] S. A. Adelman and J. D. Doll, "Generalized Langevin Equation Approach for Atom/Solid Surface Scattering: Collinear Atom/Harmonic Chain Model," J. Chem. Phys. vol. 61, pp. 42424245, 1974. [83] J. D. Doll, L. E. Myers, and S. A. Adelman, "Generalized Langevin Equation Approach for Atom/Solid Surface Scattering: Inelasti c Studies," J. Chem. Phys. vol. 63, pp. 49084914, 1975. [84] S. A. Adelman and J. D. Doll, "Generalized Langevin Equation Approach for Atom Solid Surface Scattering General Formulation for Classical Scattering Off Harmonic Solids," J. Chem. Phys. vol 64, pp. 23752388, 1976. PAGE 186 186 [85] S. Nos, "A unified formulation of the constant temperature molecular dynamics methods," J. Chem. Phys. vol. 81, pp. 511519, 1984. [86] W. G. Hoover, "Canonical Dynamics Equilibrium Phase Space Distributions," Physical Review A vol. 31, pp. 16951697, 1985. [87] G. Shirane, S. Hoshino, and K. Suzuki, "X Ray study of the phase transition in lead titanate," Phys. Rev. vol. 80, pp. 11051106, 1956. [88] G. Shirane, J. D. Axe, J. Harada, and J. P. Remeika, "Soft Ferroel ectric Modes In Lead Titanate," Phys. Rev. B vol. 2, pp. 155159, 1970. [89] S. A. Mabud and A. M. Glazer, "Lattice Parameters And Birefringence In PbTiO3 Single Crystals," J. Appl. Crystallogr. vol. 12, pp. 4953, 1979. [90] G. Burns and B. A. Scott, "Lattice Modes In Ferroelectric Perovskites PbTiO3," Phys. Rev. B vol. 7, pp. 30883101, 1973. [91] V. G. Gavrilyachenko, R. I. Spinko, M. A. Martynenko, and E. G. Fesenko, "Spontaneous Polarization And Coercive Field Of Lead Titanate," Sov. Phys. Soli d State vol. 12, pp. 12031204, 1970. [92] M. Ahart, A. Asthagiri, Z. G. Ye, P. Dera, H. K. Mao, R. E. Cohen, and R. J. Hemley, "Brillouin scattering and molecular dynamics study of the elastic properties of Pb(Mg1/3Nb2/3)O3," Phys. Rev. B vol. 75, pp. 144410, 2007. [93] G. Kresse and J. Furthmller, "Efficiency of ab initio total energy calculations for metals and semiconductors using a plane wave basis set," Comput. Mater. Sci. vol. 6, pp. 15 50, 1996. [94] G. Kresse and J. Furthmller, "Efficient i terative schemes for ab initio total energy calculations using a plane wave basis set," Phys. Rev. B vol. 54, pp. 1116911186, 1996. [95] G. Kresse and J. Hafner, "Ab initio molecular dynamics for liquid metals," Phys. Rev. B vol. 47, pp. 558561, 1993. [96] R. D. King Smith and D. Vanderbilt, "First principles investigation of ferroelectricity in perovskite compounds," Phys. Rev. B vol. 49, pp. 58285844, 1994. [97] J. Padilla and D. Vanderbilt, "Ab initio study of BaTiO3 surfaces," Phys. Rev. B vol 56, pp. 16251631, 1997. [98] J. Padilla and D. Vanderbilt, "Ab initio study of SrTiO3 surfaces," Surf. Sci. vol. 418, pp. 6470, 1998. PAGE 187 187 [99] N. Bickel, G. Schmidt, K. Heinz, and K. Mller, "Ferroelectric relaxation of the SrTiO3 (100) surface," Phys. Rev. Lett. vol. 62, pp. 20092011, 1989. [100] B. Meyer, J. Padilla, and D. Vanderbilt, "Theory of PbTiO3, BaTiO3, and SrTiO3 surfaces," Faraday Discuss. vol. 114, pp. 395405, 1999. [101] A. M. Kolpak, N. Sai, and A. M. Rappe, "Short circuit boundary conditions in ferroeletric PbTiO3 thin films," Phys. Rev. B vol. 74, pp. 054112, 2006. [102] B. Meyer and D. Vanderbilt, "Ab initio study of BaTiO3 and PbTiO3 surfaces in external electric fields," Phys. Rev. B vol. 63, pp. 205426, 2001. [103] N. Sai, A. M. Kolpak, and A. M. Rappe, "Ferroeletricity in ultrathin perovskite films," Phys. Rev. B vol. 72, pp. 020101(R), 2005. [104] L. Bengtsson, "Dipole correction for surface supercell calculations," Phys. Rev. B vol. 59, pp. 1230112304, 1999. [105] J. Neugebauer and M. Scheffler, "Adsorbate Substrate And Adsorbate Adsorbate Interactions Of Na And K Adlayers On Al(111)," Phys. Rev. B vol. 46, pp. 1606716080, 1992. [106] B. Meyer and D. Vanderbilt, "Ab initio study of ferroelectric domain walls in PbT iO3," Phys. Rev. B vol. 65, pp. 104111, 2002. [107] W. Zhong, R. D. KingSmith, and D. Vanderbilt, "Giant LO TO splittings in perovskite ferroelectrics," Phys. Rev. Lett. vol. 72, pp. 36183621, 1994. [108] G. C. Benson, P. Balk, and P. White, "Contrib ution of surface distortion to the surface energy of alkali halide crystals," J. Chem. Phys. vol. 31, pp. 109115, 1959. [109] H. Tang, X. Bouju, C. Joachim, C. Girard, and J. Devillers, "Theoretical study of the atomic force microscopy imaging process on the NaCl(001) surface," J. Chem. Phys. vol. 108, pp. 359367, 1998. [110] J. Vogt and H. Weiss, "The structure of NaCl(100) and KCl(100) single crystal surfaces: a tensor low energy electron diffraction analysis," Surf. Sci. vol. 491, pp. 155168, 2001. [111] E. Heifets, R. I. Eglitis, E. A. Kotomin, J. Maier, and G. Borstel, "Ab initio modeling of surface structure for SrTiO3 perovskite crystals," Phys. Rev. B vol. 64, pp. 235417, 2001. [112] V. Gopalan, K. L. Schepler, V. Dierolf, and I. Biaggio, "Ferroelectric Materials," in Handbook of Photonics M. C. Gupta and J. Ballato, Eds., 2nd ed. FL.: CRC Press LLC, 2006, pp. 6 1/6 66. PAGE 188 188 [113] R. C. Bradt and G. S. Ansell, "Transmission Electron Microscopy Of Ferroelectric Domain Boundaries In Barium Titan ate," IEEE Trans. Electron Devices vol. ED16, pp. 594597, 1969. [114] M. Foeth, A. Sfera, P. Stadelmann, and P. A. Buffat, "A comparison of HREM and weak beam transmission electron microscopy for the quantitative measurement of the thickness of ferroelectric domain walls," J. Electron. Microsc. vol. 48, pp. 717 723, 1999. [115] C. J. Lu, S. B. Ren, H. M. Shen, J. S. Liu, and Y. N. Wang, "The effect of grain size on domain structure in unsupported PbTiO3 thin films," J. Phys.: Condens. Matter vol. 8, p p. 80118016, 1996. [116] C. J. Lu, H. M. Shen, S. B. Ren, and Y. N. Wang, "X ray evidence for the absence of 90 degrees domains in polycrystalline PbTiO3 thin films with grains below a critical size," Appl. Phys. A: Mater. Sci. Process. vol. 65, pp. 395401, 1997. [117] S. Stemmer, S. K. Streiffer, F. Ernst, and M. Ruhle, "Atomistic Structure Of 90 Degrees Domain Walls In Ferroelectric PbTiO3 Thin Films," Philos. Mag. A vol. 71, pp. 713724, 1995. [118] S. K. Streiffer, J. A. Eastman, D. D. Fong, C. T hompson, A. Munkholm, M. V. R. Murty, O. Auciello, G. R. Bai, and G. B. Stephenson, "Observation of nanoscale 180 degrees stripe domains in ferroelectric PbTiO3 thin films," Phys. Rev. Lett. vol. 89, pp. 067601, 2002. [119] S. Venkatesan, B. J. Kooi, J. T. M. De Hosson, A. H. G. Vlooswijk, and B. Noheda, "Substrate influence on the shape of domains in epitaxial PbTiO3 thin films," J. Appl. Phys. vol. 102, pp. 104105, 2007. [120] T. J. Yang, V. Gopalan, P. J. Swart, and U. Mohideen, "Direct observation o f pinning and bowing of a single ferroelectric domain wall," Phys. Rev. Lett. vol. 82, pp. 41064109, 1999. [121] W. Cao, G. R. Barsch, and J. A. Krumhansl, "Quasi One Dimensional Solutions For Domain Walls And Their Constraints In Improper Ferroelastics ," Phys. Rev. B vol. 42, pp. 63966401, 1990. [122] X. R. Huang, X. B. Hu, S. S. Jiang, and D. Feng, "Theoretical model of 180 degrees domain wall structures and their transformation in ferroelectric perovskites," Phys. Rev. B vol. 55, pp. 55345537, 1997. [123] Y. H. Shin, I. Grinberg, I. W. Chen, and A. M. Rappe, "Nucleation and growth mechanism of ferroelectric domain wall motion," Nature vol. 449, pp. 881884, 2007. [124] J. Padilla, W. Zhong, and D. Vanderbilt, "First principles investigation of 180 degrees domain walls in BaTiO3," Phys. Rev. B vol. 53, pp. R5969R5973, 1996. PAGE 189 189 [125] S. Poykko and D. J. Chadi, "Ab initio study of 180 degrees domain wall energy and structure in PbTiO3," Appl. Phys. Lett. vol. 75, pp. 28302832, 1999. [126] G. Sag hi Szabo, R. E. Cohen, and H. Krakauer, "First principles study of piezoelectricity in PbTiO3," Phys. Rev. Lett. vol. 80, pp. 43214324, 1998. [127] T. Shimada, Y. Umeno, and T. Kitamura, "Ab initio study of stress induced domain switching in PbTiO3," Ph ys. Rev. B vol. 77, pp. 094105, 2008. [128] Y. H. Shin, V. R. Cooper, I. Grinberg, and A. M. Rappe, "Development of a bond valence molecular dynamics model for complex oxides," Phys. Rev. B vol. 71, pp. 054104, 2005. [129] R. K. Behera, B. B. Hinojosa, S. B. Sinnott, A. Asthagiri, and S. R. Phillpot, "Coupling of surface relaxation and polarization in PbTiO3 from atomistic simulation," J. Phys.: Condens. Matter vol. 20, pp. 395004, 2008. [130] S. Choudhury, Y. L. Li, N. Odagawa, A. Vasudevarao, L. Tia n, P. Capek, V. Dierolf, A. N. Morozovska, E. A. Eliseev, S. Kalinin, Y. S. Cho, L. Q. Chen, and V. Gopalan, "The influence of 180 degrees ferroelectric domain wall width on the threshold field for wall motion," J. Appl. Phys. vol. 104, pp. 084107, 2008. [131] R. D. King Smith and D. Vanderbilt, "Theory Of Polarization Of Crystalline Solids," Phys. Rev. B vol. 47, pp. 16511654, 1993. [132] Y. F. Duan, H. L. Shi, and L. X. Qin, "Studies of tetragonal PbTiO3 subjected to uniaxial stress along the caxis, J. Phys.: Condens. Matter vol. 20, pp. 175210, 2008. [133] D. A. Scrymgeour, V. Gopalan, A. Itagi, A. Saxena, and P. J. Swart, "Phenomenological theory of a single domain wall in uniaxial trigonal ferroelectrics: Lithium niobate and lithium tantalate," Phys. Rev. B vol. 71, pp. 184110, 2005. [134] D. Lee, R. K. Behera, P. Wu, H. Xu, S. B. Sinnott, S. R. Phillpot, L. Q. Chen, and V. Gopalan, "Mixed Bloch Neel Ising Character of 180 degree Ferroelectric Domain Walls," Phys. Rev. B vol. 80, pp. 060102(R ), 2009. [135] B. A. Strukov, S. T. Davitadze, S. G. Shulman, B. V. Goltzman, and V. V. Lemanov, "Clarification of size effects in polycrystalline BaTiO3 thin films by means of the specific heat measurements: Grain size or film thickness?" Ferroelectrics vol. 301, pp. 157162, 2004. [136] L. Sun, Y. F. Chen, W. H. Ma, L. W. Wang, T. Yu, M. S. Zhang, and N. B. Ming, "Evidence of ferroelectricity weakening in the polycrystalline PbTiO3 thin films," Appl. Phys. Lett. vol. 68, pp. 37283730, 1996. PAGE 190 190 [137] B. A. Strukov, S. T. Davitadze, S. N. Kravchun, S. A. Taraskin, M. Goltzman, V. V. Lemanov, and S. G. Shulman, "Specific heat and heat conductivity of BaTiO3 polycrystalline films in the thickness range 20 1100 nm," Journal of Physics Condensed Matter vol 15, pp. 43314340, 2003. [138] G. Burns and B. A. Scott, "Raman studies of underdamped soft modes in PbTiO3," Phys. Rev. Lett. vol. 25, pp. 167170, 1970. [139] G. Gottstein and L. S. Shvindlerman, Grain Boundary Migration in Metals: Thermodynamics, K inetics, Applications Boca Raton, FL: CRC Press, 1999. [140] G. Hasson, Herbeuva.I, J. Y. Boos, M. Biscondi, and C. Goux, "Theoretical And Experimental Determinations Of Grain Boundary Structures And Energies Correlation With Various Experimental Resul ts," Surf. Sci. vol. 31, pp. 115137, 1972. [141] D. Wolf and S. Yip, Materials Interfaces: Atomic level Structure and Properties : Springer, 1992. [142] N. A. Benedek, A. L. S. Chua, C. Elsasser, A. P. Sutton, and M. W. Finnis, "Interatomic potentials f or strontium titanate: An assessment of their transferability and comparison with density functional theory," Phys. Rev. B vol. 78, pp. 064110, 2008. [143] V. P. Dravid and V. Ravikumar, "Atomic structure and properties of the (310) symmetrical tilt grai n boundary (STGB) in SrTiO3 Part II: Comparison with experimental studies," Interface Sci. vol. 8, pp. 177187, 2000. [144] M. Imaeda, T. Mizoguchi, Y. Sato, H. S. Lee, S. D. Findlay, N. Shibata, T. Yamamoto, and Y. Ikuhara, "Atomic structure, electron ic structure, and defect energetics in [001](310) Sigma 5 grain boundaries of SrTiO3 and BaTiO3," Phys. Rev. B vol. 78, pp. 245320, 2008. [145] V. Ravikumar, V. P. Dravid, and D. Wolf, "Atomic structure and properties of the (310) symmetrical tilt grain boundary (STGB) in SrTiO3. Part I: Atomistic simulations," Interface Sci. vol. 8, pp. 157175, 2000. [146] M. J. Akhtar, Z. U. N. Akhtar, R. A. Jackson, and C. R. A. Catlow, "Computer Simulation Studies Of Strontium Titanate," J. Am. Ceram. Soc., vol. 78, pp. 421428, 1995. [147] J. Crawford and P. Jacobs, "Point defect energies for strontium titanate: A pair potentials study," J. Solid State Chem. vol. 144, pp. 423 429, 1999. [148] M. A. McCoy, R. W. Grimes, and W. E. Lee, "Phase stability and interfa cial structures in the SrO SrTiO3 system," Philos. Mag. A vol. 75, pp. 833846, 1997. PAGE 191 191 [149] S. Sekiguchi, M. Fujimoto, M. Nomura, S. B. Cho, J. Tanaka, T. Nishihara, M. G. Kang, and H. H. Park, "Atomic force microscopic observation of SrTiO3 polar surfac e," Solid State Ionics vol. 108, pp. 7379, 1998. [150] B. S. Thomas, N. A. Marks, and B. D. Begg, "Developing pair potentials for simulating radiation damage in complex oxides," Nucl. Instrum. Methods Phys. Res., Sect. B vol. 228, pp. 288292, 2005. [151] K. R. Udayakumar and A. N. Cormack, "NonStoichiometry In Alkaline Earth Excess Alkaline Earth Titanates," J. Phys. Chem. Solids vol. 50, pp. 55 60, 1989. [152] J. L. Wohlwend, R. K. Behera, I. Jang, S. R. Phillpot, and S. B. Sinnott, "Morphology an d growth modes of metal oxides deposited on SrTiO3," Surf. Sci. vol. 603, pp. 873880, 2009. [153] R. Wahl, D. Vogtenhuber, and G. Kresse, "SrTiO3 and BaTiO3 revisited using the projector augmented wave method: Performance of hybrid and semilocal functionals," Phys. Rev. B vol. 78, pp. 104116, 2008. [154] A. Vasudevarao, S. Denev, M. D. Biegalski, Y. Li, L. Q. Chen, S. Trolier McKinstry, D. G. Schlom, and V. Gopalan, "Polarization rotation transitions in anisotropically strained SrTiO3 thin films," Appl Phys. Lett. vol. 92, pp. 192902, 2008. [155] S. Denev, A. Kumar, M. D. Biegalski, H. W. Jang, C. M. Folkman, A. Vasudevarao, Y. Han, I. M. Reaney, S. Trolier McKinstry, C. B. Eom, D. G. Schlom, and V. Gopalan, "Magnetic color symmetry of lattice rotati ons in a diamagnetic material," Phys. Rev. Lett. vol. 100, pp. 257601, 2008. [156] C. L. Jia, V. Nagarajan, J. Q. He, L. Houben, T. Zhao, R. Ramesh, K. Urban, and R. Waser, "Unit cell scale mapping of ferroelectricity and tetragonality in epitaxial ultra thin ferroelectric films," Nat. Mater. vol. 6, pp. 6469, 2007. [157] P. Muralt, "Ferroelectric thin films The emancipation of ferroelectricity," Nat. Mater. vol. 6, pp. 89, 2007. [158] D. D. Fong and C. Thompson, "In situ synchrotron X ray studies of ferroelectric thin films," Annu. Rev. Mater. Res. vol. 36, pp. 431465, 2006. [159] G. B. Stephenson and K. R. Elder, "Theory for equilibrium 180 degrees stripe domains in PbTiO3 films," J. Appl. Phys. vol. 100, pp. 051601, 2006. [160] E. Cockayne a nd B. P. Burton, "Dipole moment of a PbO vacancy pair in PbTiO3," Phys. Rev. B vol. 69, pp. 144116, 2004. [161] L. X. He and D. Vanderbilt, "First principles study of oxygenvacancy pinning of domain walls in PbTiO3," Phys. Rev. B vol. 68, pp. 134103, 2003. PAGE 192 192 [162] H. Mestric, R. A. Eichel, T. Kloss, K. P. Dinse, S. Laubach, S. Laubach, P. C. Schmidt, K. A. Schonau, M. Knapp, and H. Ehrenberg, "Iron oxygen vacancy defect centers in PbTiO3: Newman superposition model analysis and density functional calcul ations," Phys. Rev. B vol. 71, pp. 134109, 2005. [163] C. H. Park and D. J. Chadi, "Microscopic study of oxygenvacancy defects in ferroelectric perovskites," Phys. Rev. B vol. 57, pp. R13961R13964, 1998. [164] C. H. Park and D. J. Chadi, "Effect of i nterstitial hydrogen impurities on ferroelectric polarization in PbTiO3," Phys. Rev. Lett. vol. 84, pp. 47174720, 2000. [165] J. G. Yu, S. B. Sinnott, and S. R. Phillpot, "Charge optimized many body potential for the Si/SiO2 system," Phys. Rev. B vol. 75, pp. 085311, 2007. PAGE 193 193 BIOGRAPHICAL SKETCH Rakesh Kumar Behera was born in Cuttack, Orissa, India. He rece ived his Bachelors degree in metallurgical e ngineering from Regional Engineering College (currently known as National Institute of Technology), Durgapur India in 2000. He came to the United States for his higher studies and joined the graduate program in Mechanical Engineering Department at Louisiana State University (LSU), Baton Rouge, Louisiana, in Fall 2002. H e recei ved his Maste r of Science in mechanical engineering from LSU in August 2004. After completion, he joined the doctoral degree under Prof. Simon R. Phillpot at the University of Florida (UF) Gainesville. He received his M.S and Ph.D. in materials science an d engineering from t he University of Florida in spring of 2008 and fall of 2009 respectively After joining the respective universities, he brought a lot of good charm to both LSU and UF. LSU own the 2003 and 2007 BCS national championships in football an d 2009 national championship in baseball. LS U became the first school to win two BCS national championship s UF also did extremely well by winning the 2006 and 2008 BCS national championships in football and en route to get another chance for a back to bac k national title opportunity. UF became the first school in NCAA history to win both football and basketball national championships in the same year (2006) and won the 2007 basketball national champions. He believe s the dominating performance by LSU and UF will continue in future. Finally He is happy to say Geaux Tigers and Go Gators 