<%BANNER%>

Evolution of Intrinsic Point Defects in Fluorite-Based Materials

Permanent Link: http://ufdc.ufl.edu/UFE0024282/00001

Material Information

Title: Evolution of Intrinsic Point Defects in Fluorite-Based Materials Insights from Atomic-Level Simulation
Physical Description: 1 online resource (197 p.)
Language: english
Creator: Aidhy, Dilpuneet
Publisher: University of Florida
Place of Publication: Gainesville, Fla.
Publication Date: 2009

Subjects

Subjects / Keywords: bi2o3, ceramics, dft, fluorite, md, mgo, molybdenum, radiation, simulation, sofc, uo2
Materials Science and Engineering -- Dissertations, Academic -- UF
Genre: Materials Science and Engineering thesis, Ph.D.
bibliography   ( marcgt )
theses   ( marcgt )
government publication (state, provincial, terriorial, dependent)   ( marcgt )
born-digital   ( sobekcm )
Electronic Thesis or Dissertation

Notes

Abstract: Point defects play a significant role in materials properties. On the one hand, point defects are exploited in applications such as fuel cells; on the other hand, they are deleterious to performance of in such applications as nuclear fuels. To tailor materials with desired properties, it is thus important to understand point defect behavior in materials under specific environments/applications. This thesis is focused on understanding the evolution of point defects in materials, primarily fluorite-related, used in solid oxide fuel cell (SOFC) and nuclear applications. Molecular-dynamics simulation (MD) and density functional theory (DFT) are used as materials modeling tools. Cubic bismuth oxide (delta-Bi2O3) is a promising SOFC electrolyte and a fluorite-based model material, in which, it is found that while intrinsic oxygen vacancies contribute in high oxygen diffusion, under certain circumstances, they themselves act as limiting factors by forming a < 110 > - < 111 > vacancy-ordered system. It is found that high cationic polarizability plays a significant role in achieving high oxygen diffusivity. It is also found that the oxygen diffusivity may be limited due to the formation of covalent character bonds by some Bi ions. Fluorite-structured UO2 is the fuel in almost all operating nuclear reactors. Due to irradiation, UO2 undergoes damage resulting in the formation of Frenkel pairs (FPs) on both uranium and oxygen sub-lattices. It is found that, while damage on oxygen sub-lattice alone is not deleterious, in that oxygen FPs annihilate by mutual recombination, presence of FPs on uranium sub-lattice cause long-lasting damage and make UO2 less radiation tolerant. Regardless of FPs on oxygen sub-lattice, FPs on uranium sub-lattice nucleate new O FPs, which form clusters. The oxygen vacancies? sequestration by uranium vacancies is found to be the mechanism of cluster formation. The point-defect evolution is rocksalt MgO is also elucidated. MgO is an important engineering material and has nuclear applications. As in UO2, it is found that when defects are present on both Mg and O sub-lattice, new FPs, in this case, both on Mg and O are formed. Defects present only on either Mg or O sub-lattice annihilate by mutual recombination. While UO2 and MgO studies are done on single crystals to elucidate FPs elimination by mutual recombination only, to understand their evolution in the presence of GBs, similar studies are done on nanocrystalline Mo. This allows elucidation of dominant mechanism of point-defect elimination. It is found that for the grain size and temperature under study, mutual recombination of FPs dominates and GBs have little effect on the point-defect elimination.
General Note: In the series University of Florida Digital Collections.
General Note: Includes vita.
Bibliography: Includes bibliographical references.
Source of Description: Description based on online resource; title from PDF title page.
Source of Description: This bibliographic record is available under the Creative Commons CC0 public domain dedication. The University of Florida Libraries, as creator of this bibliographic record, has waived all rights to it worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law.
Statement of Responsibility: by Dilpuneet Aidhy.
Thesis: Thesis (Ph.D.)--University of Florida, 2009.
Local: Adviser: Phillpot, Simon R.

Record Information

Source Institution: UFRGP
Rights Management: Applicable rights reserved.
Classification: lcc - LD1780 2009
System ID: UFE0024282:00001

Permanent Link: http://ufdc.ufl.edu/UFE0024282/00001

Material Information

Title: Evolution of Intrinsic Point Defects in Fluorite-Based Materials Insights from Atomic-Level Simulation
Physical Description: 1 online resource (197 p.)
Language: english
Creator: Aidhy, Dilpuneet
Publisher: University of Florida
Place of Publication: Gainesville, Fla.
Publication Date: 2009

Subjects

Subjects / Keywords: bi2o3, ceramics, dft, fluorite, md, mgo, molybdenum, radiation, simulation, sofc, uo2
Materials Science and Engineering -- Dissertations, Academic -- UF
Genre: Materials Science and Engineering thesis, Ph.D.
bibliography   ( marcgt )
theses   ( marcgt )
government publication (state, provincial, terriorial, dependent)   ( marcgt )
born-digital   ( sobekcm )
Electronic Thesis or Dissertation

Notes

Abstract: Point defects play a significant role in materials properties. On the one hand, point defects are exploited in applications such as fuel cells; on the other hand, they are deleterious to performance of in such applications as nuclear fuels. To tailor materials with desired properties, it is thus important to understand point defect behavior in materials under specific environments/applications. This thesis is focused on understanding the evolution of point defects in materials, primarily fluorite-related, used in solid oxide fuel cell (SOFC) and nuclear applications. Molecular-dynamics simulation (MD) and density functional theory (DFT) are used as materials modeling tools. Cubic bismuth oxide (delta-Bi2O3) is a promising SOFC electrolyte and a fluorite-based model material, in which, it is found that while intrinsic oxygen vacancies contribute in high oxygen diffusion, under certain circumstances, they themselves act as limiting factors by forming a < 110 > - < 111 > vacancy-ordered system. It is found that high cationic polarizability plays a significant role in achieving high oxygen diffusivity. It is also found that the oxygen diffusivity may be limited due to the formation of covalent character bonds by some Bi ions. Fluorite-structured UO2 is the fuel in almost all operating nuclear reactors. Due to irradiation, UO2 undergoes damage resulting in the formation of Frenkel pairs (FPs) on both uranium and oxygen sub-lattices. It is found that, while damage on oxygen sub-lattice alone is not deleterious, in that oxygen FPs annihilate by mutual recombination, presence of FPs on uranium sub-lattice cause long-lasting damage and make UO2 less radiation tolerant. Regardless of FPs on oxygen sub-lattice, FPs on uranium sub-lattice nucleate new O FPs, which form clusters. The oxygen vacancies? sequestration by uranium vacancies is found to be the mechanism of cluster formation. The point-defect evolution is rocksalt MgO is also elucidated. MgO is an important engineering material and has nuclear applications. As in UO2, it is found that when defects are present on both Mg and O sub-lattice, new FPs, in this case, both on Mg and O are formed. Defects present only on either Mg or O sub-lattice annihilate by mutual recombination. While UO2 and MgO studies are done on single crystals to elucidate FPs elimination by mutual recombination only, to understand their evolution in the presence of GBs, similar studies are done on nanocrystalline Mo. This allows elucidation of dominant mechanism of point-defect elimination. It is found that for the grain size and temperature under study, mutual recombination of FPs dominates and GBs have little effect on the point-defect elimination.
General Note: In the series University of Florida Digital Collections.
General Note: Includes vita.
Bibliography: Includes bibliographical references.
Source of Description: Description based on online resource; title from PDF title page.
Source of Description: This bibliographic record is available under the Creative Commons CC0 public domain dedication. The University of Florida Libraries, as creator of this bibliographic record, has waived all rights to it worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law.
Statement of Responsibility: by Dilpuneet Aidhy.
Thesis: Thesis (Ph.D.)--University of Florida, 2009.
Local: Adviser: Phillpot, Simon R.

Record Information

Source Institution: UFRGP
Rights Management: Applicable rights reserved.
Classification: lcc - LD1780 2009
System ID: UFE0024282:00001


This item has the following downloads:


Full Text

PAGE 1

1 EVOLUTION OF INTRINSIC POINT DEFECTS IN FLUORITE BASED MATERIALS: INSIGHTS FROM ATOMIC -LEVEL SIMULATION By DILPUNEET SINGH AIDHY A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FUL FILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2009

PAGE 2

2 2009 Dilpuneet S. Aidhy

PAGE 3

3 To my parents with love and gratitude

PAGE 4

4 ACKNOWLEDGMENTS At this very juncture of my life, where I now bridge my educational and professional career, I want to take this opportunity to thank some very important people who have made my PhD a wonderful journey. I express my deepest gratitude towards my mentor Prof. Simon Phillpot, who se advising skills, besides scientific knowledge, I find to be exceptional. His enthusiasm towards science, trust and patience with students blended with a great sense of humor has made my PhD thoroughly enjoyable. He has always provided liberty to explore things even not directly related to research and has shared his wisdom. I feel fortunate to have him as an advisor, responsibility of which, he fulfilled perfectly. I also want to thank Dr. Dieter Wolf with whom I have had many thought -provoking discussi ons from which this thesis has benefited immensely. I have never met a more humble person than him. I admire his love for nature and thank him for stimulating me to explore western US. I also thank Dr. Paul Millet, Dr. Tapan Desai and Srujan, with whom I e njoyed our collaborative work on nuclear materials. Thanks to the friends I made in Idaho, their generosity and warmth softened the blow of Idaho winters. Many thanks to my committee members, Profs. Susan Sinnott, Juan Nino, Eric Wachsman and Jason Weav er for serving on my committee and providing valuable suggestions to my research work. Special thanks to Prof. Sinnott for co advising me, Prof. Nino for teaching me crystallography and Prof. Wachsman for keeping me honest from experimental standpoint. I want to mention some of the very dear friendships that I feel fortunate to have developed during the years. Special thanks to Tabassum for giving me moral support especially during the semester -ending result times, Rakesh for his always -welcome attitude an d

PAGE 5

5 exceptionally delicious food, and Sankara for being a perfect roommate. Besides, all current and past members of the Sinpot group, thanks also to Priyank, Chan Woo, Taku, Peter and Tao. All of you did bear me very kindly. In the end, I want to mention m y dad and my mom, to thank whom I do not have enough words. They have been the supporting pillars of my life and have never let the distance of seven seas interfere in keeping their warmth of love around me. I am indebted to you, dad, especially for lettin g me to pursue a degree in the US, exactly a month after you underwent heart surgery, and you mom, for being especially strong during that time. I do not want to forget my little sister who throughout these years looked forward more to souvenirs than me. It would be inappropriate to not to thank the Almighty for making this journey easy and fruitful.

PAGE 6

6 TABLE OF CONTENTS Page ACKNOWLEDGMENTS .................................................................................................................... 4 LIST OF TABLES .............................................................................................................................. 10 LIST OF FIGURES ............................................................................................................................ 11 ABSTRACT ........................................................................................................................................ 15 CHAPTER 1 GENERAL INTRODUCTION .................................................................................................. 17 1.1 Sources of Energy ................................................................................................................. 17 1.2 Energy and Environment ...................................................................................................... 19 1.3 Cl ean Technology ................................................................................................................. 20 1.4 Challenges in Materials Science .......................................................................................... 22 1.5 Motivation for This Work .................................................................................................... 23 1.6 Significance of Simulation ................................................................................................... 24 1.7 Objective of This Work ........................................................................................................ 25 1.8 Organization of Dissertation ................................................................................................ 26 2 DEFECTS IN SOLIDS ............................................................................................................... 27 2.1 Introduction ........................................................................................................................... 27 2.2 Intrinsic Point Defe cts .......................................................................................................... 27 2.2.1 Frenkel Defect ............................................................................................................ 28 2.2.2 Schottky Defect .......................................................................................................... 29 2.2.3 Th ermodynamics of Point Defects: Equilibrium Concentration ............................. 29 2.3 Transport ................................................................................................................................ 32 2.3.1 Atomistic Theory of Diffusion .................................................................................. 32 2.3.2 Diffusion Mechanisms ............................................................................................... 34 2.3.2.1 The vacancy mechanism ................................................................................. 34 2.3.2.2 The i nterstitial mechanism .............................................................................. 34 2.3.2.3 The interstitialcy mechanism .......................................................................... 35 2.4 Crystallography of Materials ................................................................................................ 36 2.4.1 Fluorite Crystal Structure: Bi2O3 and UO2 ............................................................ 36 2.4.2 Rocksalt Crystal Structure: MgO .............................................................................. 37 2.4.3 Space Group Fm 3 m .................................................................................................. 37 2.4.4 Body-Centered Cubic: Space Group Im 3 m ............................................................ 39 3 SIMULATION M ETHODOLOGY ........................................................................................... 40 3.1 Simulation Methodologies ................................................................................................... 40 3.2 Molecular Dynamics Simulation ......................................................................................... 41

PAGE 7

7 3.3 General MD Algorithm ........................................................................................................ 41 3.4 Periodic Boundary Condition ............................................................................................... 44 3.5 Interatomic Interactions ........................................................................................................ 45 3.5.1 Long-range Interactions ............................................................................................. 47 3.5.2 Short -Range Interactions ........................................................................................... 53 3.5.3 Elect ronic Polarizability Shell Model .................................................................... 55 3.5.4 Thermodynamic Conditions (ensemble) for Materials Simulation ......................... 56 4 IONIC CONDUCTIVITY OF SOLID OXIDE FUEL CELL ELECTROLYTE................... 58 4.1 Fuel Cell ................................................................................................................................ 58 4.2 Solid Oxide Fuel Cell ........................................................................................................... 59 4.3 Fluorite Based Cubic Bismuth Oxide .................................................................................. 61 4.3.1 A Model Material ....................................................................................................... 61 4.3.2 Bi2O3 Limitations .................................................................................................... 63 4.3.3 Ionic Radii or Polarizability? ..................................................................................... 64 4.4 Simulation Methodology ...................................................................................................... 66 4. 5 Results .................................................................................................................................... 67 4.5.1 Polarizability Indeed .................................................................................................. 67 4.5.1.1 Polarizable and nonpolarizable ions .............................................................. 69 4.5.1.2 Onlyoxygen ions polarizable ......................................................................... 71 4.5.1.3 Onlybismuth ions polarizable ........................................................................ 72 4.5.2 Oxygen Diffus ion Mechanism ................................................................................... 73 4.6 Discussion .............................................................................................................................. 75 5 VACANCY ORDERED STRUCTURE OF CUBIC BISMUTH OXIDE ............................. 77 5.1 Introduction ........................................................................................................................... 77 5.2 Ordering of the Vacancies: MD Simulation ........................................................................ 80 5.2.1 Neighboring Vacancies of a Va cancy ....................................................................... 81 5.2.2 Continuous Ordered Structure of Vacancies ............................................................ 83 5.2.3 Expansion and Contraction of the Tetrahedron ........................................................ 85 5.3 Density Functional Theory (DFT) Calculations of the Structure of Bi2O3 ....................... 86 5.3.1 DFT Methodology ...................................................................................................... 86 5.3.2 Vacancy Ordering....................................................................................................... 87 5.3.3 Bonding and Charge Distribution.............................................................................. 88 5.3.4 Electronic Structure of Vacanc y -Ordered -Bi2O3 .................................................. 91 5.4 Discussion .............................................................................................................................. 93 6 CRYSTALLOGRAPHIC ANALYSIS OF THE STRUCTURE OF CUBIC BISMUTH OXIDE ......................................................................................................................................... 94 6.1 Introduction ........................................................................................................................... 94 Bi2O3 ....................................................... 94 6. 2.1 Space Group Analysis ................................................................................................ 97 6.2.2 Bond Lengths .............................................................................................................. 99 6.2.2.1 Bond lengths from DFT calculations ........................................................... 100

PAGE 8

8 6.2.2.2 Cation displacements in the <111> vacancy ordered oxygen sublattices 100 6.2.2.3 Cation displacements in the <110> vacancy ordered oxygen sublattice .. 101 6.2.2.4 Anion displacement ....................................................................................... 102 6.2.2.5 Bond lengths from Shannon radii ................................................................. 103 6.2.2.6 Subunit Cell Distortion ........................................................................................ 103 6.2.2.7 Comparison with Bixbyite Structure .................................................................... 106 6.3 Relation to Experiment ....................................................................................................... 107 6.4 Discussion ............................................................................................................................ 109 6.5 Conclusion and Outlook ..................................................................................................... 110 7 RADIATION DAMAGE IN MgO .......................................................................................... 112 7.1 Nuclear Energy.................................................................................................................... 112 7.1.1 Nuclear Power .......................................................................................................... 112 7.1.2 Challenges for Materials Science ............................................................................ 112 7.2 Radioactivity ....................................................................................................................... 113 7.2.1 Alpha ( ) Decay ....................................................................................................... 114 7.2.2 Beta ( ) Decay .......................................................................................................... 115 7.2.3 Gamma ( ) Decay ..................................................................................................... 116 7.2.4 Neutron Capture ....................................................................................................... 116 7.3 Nuclear Fission.................................................................................................................... 116 7.4 Radiation Damage ............................................................................................................... 117 7.4.1 Standard Radiation Damage Method ology Using MD Simulation ....................... 118 7.4.2 Kinetically-evolving Irradiation induced Defects Method .................................... 120 7.4.3 Materials Under Study ............................................................................................. 122 7.5 Interatomic Potential and MD Simulation ......................................................................... 123 7.6 Results .................................................................................................................................. 124 7.6.1 FPs Only on One Sublattice ..................................................................................... 124 7.6.2 FPs on Both Sublattices ........................................................................................... 125 7.7 Cluster -Formation Mechanism ........................................................................................... 128 7.7.1 Crystallography of the Cluster ................................................................................. 130 7.8 Discussion ............................................................................................................................ 132 8 RADIATION DAMAGE IN UO2 ............................................................................................ 135 8.1 Introduction ......................................................................................................................... 135 8.2 Simulation Methodology .................................................................................................... 137 8.2.1 I nteratomic Potential and Molecular Dynamics Method ....................................... 137 8.2.2 Diffusion-controlled Kinetic Evolution of Defects ................................................ 138 8.2.3 Frenkel -pair Defects on a Single Sub -lattice .......................................................... 139 8.2.4 FPs on Both U and O Sub-lattices ........................................................................... 142 8.2.5 Vacancy Clustering: Schottky Defects ................................................................... 143 8.2.6 Interstitial Clustering: Cuboctahedral Clusters ...................................................... 145 8.2.7 Evolution of Disproportionately Defected Systems ............................................... 150 8.3 Discussion ............................................................................................................................ 153 9 EVOLUTION OF POINT DEFECTS IN BCC MOLYBDENUM1 ...................................... 157

PAGE 9

9 9.1 Introduction ......................................................................................................................... 157 9.2 Simulation Methodology .................................................................................................... 158 9.3 Vacancy and Interstitial Diffusivities ................................................................................ 159 9.4 Grain Boundaries as Sources and Sinks ............................................................................ 161 9.4.1 Grain Boundaries as Vacancy Sources ................................................................... 161 9 .4.2 Source -strength Analysis ......................................................................................... 164 9.4.3 Grain Boundaries as Vacancy Sinks ....................................................................... 167 9.4.4 Grain Boundaries as Interstitial Sinks ..................................................................... 169 9.4.5 Discussion on GB Source/Sink Strength ................................................................ 171 9.5 Frenkel -pair Annihilation ................................................................................................... 172 9.5.1 FP Annihilation by Mutual Recombination ............................................................ 172 9.5.2 Frenkel Pair Annihilation in a Nanocrystalline Structure ...................................... 174 9.6 Conclusion ........................................................................................................................... 177 10 SUMMARY AND FUTURE WORK ..................................................................................... 178 10.1 Bi2O3............................................................................................................................... 178 10.1.1 Cation Polarizability ............................................................................................... 178 10.1.2 Oxygen Diffusion Mechanism .............................................................................. 178 10.1.3 Bonding ................................................................................................................... 179 10.1.4 Crystallographic Analysis of Structure ................................................................. 179 10.2 MgO, UO2 and Mo ............................................................................................................ 180 10.2.1 Methodology ........................................................................................................... 180 10.2.2 MgO Interstitial Clustering .................................................................................... 180 10.2.3 UO2 Clustering ....................................................................................................... 181 10.2.4 Defect Annihilation in Nanocrystalline Mo ......................................................... 182 10.3 Conclusions and Future Work .......................................................................................... 182 10.3.1 Bi2O3 .................................................................................................................... 182 10.3.2 Radiation Damage .................................................................................................. 185 LIST OF REFERENCES ................................................................................................................. 188 BIOGRAPHICAL SKETCH ........................................................................................................... 197

PAGE 10

10 LIST OF TABLES Table page 4 1 Characteristics of different fuel cells. ................................................................................... 59 4 2 Parameters for the short range potential for -Bi2O3 from Jacobs and Mac Donaill ........ 67 5 1 Relative stability of the ordered structures in 2x2x2 system. .............................................. 88 5 2 Bi O bond length in Bi -centered <110> and <111> vacancy-ordered oxygen sub lattice ....................................................................................................................................... 89 5 3 Electronic charge on the ionic species. ................................................................................. 91 6 1 Bi2O3 ......................................................... 98 6 2 Bi O bond lengths in the <110> and <111> vacancy ordered s ublattice ........................ 105

PAGE 11

11 LIST OF FIGURES Figure page 1 1 Worlds primary energy consumption by fuel -type in year 1973 and year 2004 .............. 17 1 2 Concentration of CO2 in the atmosphere over the last millennium. ................................... 19 1 3 Increase in the earth temperature over the last 150 years. ................................................... 20 2 1 Point defects. .......................................................................................................................... 28 2 2 Point defects in ionic materials. ............................................................................................ 28 2 3 Diffusion by vacancy mechanism. ........................................................................................ 34 2 4 Diffusion by interstitial mechanism ...................................................................................... 35 2 5 Diffusion by interstitialcy mechanism. ................................................................................. 35 2 6 Fluorite crystal structure. ....................................................................................................... 36 2 7 Rocksalt crystal structure. ...................................................................................................... 37 2 8 T hree different symmetry operations .................................................................................... 38 2 9 Body -centered crystal struc ture ............................................................................................. 39 3 1 Two -dimensional periodic system of a simulation cell ....................................................... 44 3 2 Inter ionic potential between Zr4+ and O2 -............................................................................ 46 3 3 Total Coulombic energy ........................................................................................................ 49 3 4 Charg e on the truncat ed sphere ............................................................................................. 50 3 5 E nergies of the charge -neutralized and charged systems .................................................... 51 3 6 E xact Madelung energy and approxi mate Madelung energy .............................................. 51 3 7 D amped charge-neutralized and undamped charge -neutralized energy ............................. 52 4 1 Working of a SOFC ............................................................................................................... 59 4 2 Conductivity versus temperature relationship ...................................................................... 60 4 3 Fluorite crystal structure. ....................................................................................................... 62 4 4 Relative conductivity of Bi2O3 doped with different dopants .......................................... 65

PAGE 12

12 4 5 Comparison of time constant for different doped Bi2O3 systems. ................................... 65 4 6 Ionic radius and polarizability of different dopants added in Bi2O3................................ 68 4 7 Oxygen ion diffusion thr ough two Bi ions ........................................................................... 68 4 8 Gradually increasing Bi polarizability .................................................................................. 69 4 9 MSD of oxygen in polarizable and non-polarizable Bi2O3 .............................................. 70 4 10 MSD of the oxygen ions with gradually increasing oxygen polarizability ........................ 71 4 11 MSD of the oxygen ions with gradually increasing bismuth polarizability ....................... 73 4 12 Possible oxygen diffusion directions .................................................................................... 74 4 13 Oxygen diffusion mechanism ................................................................................................ 75 5 1 Three different models of -Bi2O3 ........................................................................................ 79 5 2 Vacancy vacancy pair -distribution function ........................................................................ 81 5 3 Ordering of vacancies ............................................................................................................ 82 5 4 A continuous network of vacancies in the oxygen sublattice ............................................ 83 5 5 2x2x2 fluorite superstructure with vacancy ordering in <110> and <111> dir ections ...... 84 5 6 Bi Bi pair distribution function ............................................................................................. 8 5 5 7 Tetrahedron expansion and contraction ................................................................................ 86 5 8 Possible vacancy Bi2O3 ............................................................. 87 5 9 Two kinds of Bi centered oxygen sub lattice cubes ............................................................ 89 5 10 Electron lone pair charge around two different Bi -centered cubes .................................... 90 5 11 Bi O bonds in a tetrahedra in the unit cell ............................................................................ 92 6 1 The 2x2x2 fluorite related superstructure viewed from three different planes .................. 95 6 2 Bi2O3 along [001] direction .................................. 95 6 3 Layered sequence of the oxygen sub-Bi2O3. ...................... 96 6 4 Crystallographic sites of the space group Fm 3 ................................................................... 99 6 5 <111> vacancy ordered oxygen sub-lattice. ....................................................................... 100

PAGE 13

13 6 6 <110> vacancy -ordered oxygen sublattice. ....................................................................... 102 6 7 Distortion in the superstructure ........................................................................................... 104 7 1 Two phases of radiation damage. ........................................................................................ 119 7 2 Typical Frenkel pair (FP) defects during radiation damage using PKA. .......................... 119 7 3 Evolution of FPs on Mg sublattice. ..................................................................................... 125 7 3 Evolution of FPs on O sublattice. ........................................................................................ 125 7 5 Evolution of FPs on both sublattices ................................................................................... 126 7 6 Total number of Mg FPs present in the simulation cell ..................................................... 127 7 7 Close up of the all interstitial cluster .................................................................................. 128 7 8 Two distinct mechanisms control ling the growth of an interstitial cluster ....................... 129 7 9 A (MgO)4 octopolar.. ........................................................................................................... 131 7 10 Crystallography of MgO rocksalt lattice.. .......................................................................... 132 8 1 -U4O9 crystal structure. ...................................................................................................... 136 8 2 Evolution of the point defects prese nt only on the oxygen sublattice ............................. 139 8 3 N umber of O FPs present when the system is initialized with three different defect conditions .............................................................................................................................. 140 8 4 Evolution of the point defects present only on the uranium sublattice ........................... 141 8 5 Evolution of the point defects simultaneously present on both oxygen and uranium sub -lattices ............................................................................................................................ 142 8 6 Schottky defect formation ................................................................................................... 144 8 7 T hree different kinds of cuboctahedral (COT) clusters ..................................................... 145 8 8 P rogressive formation cuboctahedral (COT) cluster ......................................................... 147 8 9 Clustering of oxygen interstitials in the absence of vacancies .......................................... 148 8 10 Clustering of U and O vacancies ......................................................................................... 149 8 11 O FP evolution in a system initially containing 200 O FPs and fixed U FPs ................... 150 8 12 N umber of O FPs as a function of number of initially created U FPs .............................. 152

PAGE 14

14 9 1 MSD for single -crystal structures ....................................................................................... 160 9 2 Arrhenius plot for Di and Dv ................................................................................................ 161 9 3 V acancy population within a single d = 20 nm grain at T = 2900 K ................................ 162 9 4 N ucleation of two lattice vacancies ..................................................................................... 163 9 5 Evolution of the vacancy concentration in the under and over -saturated nanocrystalline structures .................................................................................................... 166 9 6 V acancy population in a single grain of the over -saturated sample.................................. 168 9 7 O versaturated interstitial distribution .................................................................................. 169 9 8 Quantification of the time evolution of the self interstitial concentration ....................... 170 9 9 Frenkel pair distribution in a single crystal. ....................................................................... 173 9 10 Recombination rate in a single crystal ................................................................................ 173 9 11 Evolving nanocrystalline structures over -saturated with the FPs ..................................... 175 9 12 V acancy and interstitial concentrations for the nanocrystalline structure ........................ 176

PAGE 15

15 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fu lfillment of the Requirements for the Degree of Doctor of Philosophy EVOLUTION OF INTRINSIC POINT DEFECTS IN FLUORITE BASED MATERIALS: INSIGHTS FROM ATOMIC -LEVEL SIMULATION By Dilpuneet Singh Aidhy May 2009 Chair: Simon R. Phillpot Major: Materials Scie nce and Engineering Point defects play a significant role in materials properties. On the one hand, point defects are exploited in applications such as fuel cells; on the other hand, they are deleterious to performance of in such applications as nuclear f uels. To tailor materials with desired properties, it is thus important to understand point defect behavior in materials under specific environments/applications. This thesis is focused on understanding the evolution of point defects in materials, primaril y fluorite related, used in solid oxide fuel cell (SOFC) and nuclear applications. Molecular dynamics simulation (MD) and density functional theory (DFT) are used as materials modeling tools. Cubic bismuth oxide ( Bi2O3) is a promising SOFC electrolyte and a fluorite -based model material, in which, it is found that while intrinsic oxygen vacancies contribute in high oxygen diffusion, under certain circumstances, they themselves act as limiting factors by forming a <110> <111> vacancy -ordered system. It is found that high cationic polarizability plays a significant role in achieving high oxygen diffusivity. It is also found that the oxygen diffusivity may be limited due to the formation of covalent character bonds by some Bi ions. Fluorite -structured UO2 is the fuel in almost all operating nuclear reactors. Due to irradiation, UO2 undergoes damage resulting in the formation of Frenkel pairs (FPs) on both

PAGE 16

16 uranium and oxygen sub lattices. It is found that, while damage on oxygen sub -lattice alone is not deleter ious, in that oxygen FPs annihilate by mutual recombination, presence of FPs on uranium sublattice cause long lasting damage and make UO2 less radiation tolerant. Regardless of FPs on oxygen sublattice, FPs on uranium sublattice nucleate new O FPs, whic h form clusters. The oxygen vacancies sequestration by uranium vacancies is found to be the mechanism of cluster formation. The point defect evolution is rocksalt MgO is also elucidated. MgO is an important engineering material and has nuclear applicati ons. As in UO2, it is found that when defects are present on both Mg and O sub lattice, new FPs, in this case, both on Mg and O are formed. Defects present only on either Mg or O sub lattice annihilate by mutual recombination. While UO2 and MgO studies are done on single crystals to elucidate FPs elimination by mutual recombination only, to understand their evolution in the presence of GBs, similar studies are done on nanocrystalline Mo. This allows elucidation of dominant mechanism of point -defect eliminat ion. It is found that for the grain size and temperature under study, mutual recombination of FPs dominates and GBs have little effect on the point -defect elimination.

PAGE 17

17 CHAPTER 1 GENERAL INTRODUCTION 1.1 Sources of Energy Availability of different forms of energy is essential for the development of an economy. Affordable energy has enabled industrialization and human development. Since industrialization, man has used different forms of energy and explored more. On the one hand, the pursuit of energy has enab led the discovery of new materials, on the other, materials have been used to produce energy and convert it from one form to another; nuclear materials and photovoltaic silicon cells are respectively two good examples. Energy and materials therefore hold a very close relationship. There are different fuels for energy production depending on their availability and affordability, e.g., coal, oil, hydro, nuclear. Figure 1 1 gives a graphical view of the energy derived from different sources/materials over the last 30 years.1 Figure 1 1. Worlds primary energy consumption by fuel type in year 1973 and year 2004 (adapted from reference 1 ). It shows the worlds significant dependence on the fossil fuels. A comparison shows that there is only modest change of the energy consumed from one source to another. Coal and oil have been used as the major sources of energy for a

PAGE 18

18 long period of time. Even now, coal continues to contribute the same share, primarily due to its widespread availability and developed infrastructure. In contrast, the dependence on oil has decreased by 10 % and the nuclear has risen by 6 %. Nevertheless, overall the graph shows that there is still significant dependence on the fossil fuels With fossil fuels being so indispensable, their depletion is worrisome as more countries industrialize. The energy consumption of the world has increased 20 times from year 1850 to 2000. In 2005, total energy consumption was estimated to be around 447 qu ads (quadrillion BTUs, or quads). If we continue to use resources at the same rate, the world is expected to run out of fossil fuels before end of the century. While coal is still expected to last for another 164 years,2 other sources of energy are not so abundantly available. Uranium deposits are only sufficient for another 70 80 years.3 Similarly, with a 2% increase in production per year, the oil resources are likely to peak in the year 2016.4 The energy consumption is further compounded by the highly populous and fast growing economies, in particular China and India. The demand for energy will increase tremendously as they make strides towards development. To put the energy situation into perspective, if people in these two countries were to enjoy the same transport -lifestyle as in the United States, the demand for oil would increase by 100 million barrels per day (BPD) above the present 83 million BPD.5 This exemplifies the enormous energy challenge that the world faces in the near future. To fulfill t he energy needs of future generations, it is important to develop new methods and sources of energy production. In addition to the future energy problem, the world also faces an imminent problem of even bigger magnitude the challenge to maintain Earths environment clean.

PAGE 19

19 1.2 Energy and Environment The reliance on fossil fuels is causing long lasting damage to the environment. All combustible fossil fuels result in the production of greenhouse gases such as CO2, CH4 and N2O. A major environmental concer n with the emission of the greenhouse gases is entrapment of the additional energy in the environment at levels far above the pre industrial levels. It is now widely accepted that the heat entrapment directly contributes to climate change; the most obvious effect is global warming. Figure 1 2 shows the concentration of the CO2 in the atmosphere over the last millennium.6 Over the period 10001800, the CO2 concentration remained virtually constant. In contrast, with the industrialization since then, the CO2 concentration has gradually increased and is now reaching alarming concentrations. Figure 1 2. Concentration of CO2 in the atmosphere over the last millennium. The concentration of CO2 has increased exponentially over the last 50 100 years, matchin g the industrialization around the world. Red line indicates sharp increase in the last 150 years. (Adapted from reference 6 ) This increase in the CO2 concentration has led to global warming. It has been show n from various scientific studies that the earth temperature increased by 0.6 0.2 C in the 20th century.3 The plot in Figure 1 3 illustrates the increase in the temperature of

PAGE 20

20 the Earth over a period of 150 years. There h as been 0.2 C increase per decade in the earths temperature over the last 30 years. The Intergovernmental Panel for Climate Change estimates an increase of 1.8 4 C in the Earths temperature within the next century.7 Another simulation study predicts tha t even with the greenhouse gases levels stabilized to year 2000 level, by the end of 21st century, the earths temperature will increase by 0.5 C.8 Clearly, the unchecked emissions of the greenhouse gases will severely raise sea levels and pose serious th reats to people living on the coastal cities such as London, New York, Mumbai, Shanghai and other low -lying areas. Figure 1 3. Increase in the earth temperature over the last 150 years. The increase in the temperature corresponds with the increa se in the concentration of the CO2. Red line indicates alarming increase in the last 30 years. (Adapted from reference 6 ) 1.3 Clean Technology The projected impacts of CO2 on the environment indicate that sig nificant and immediate measures are required to save the planet from further environmental degradation. While there are measures for developing new technologies and improving upon current ones, reducing CO2 concentration by alternative means is also among the prospective solutions. One of the proposed ways of reducing CO2 from the environment

PAGE 21

21 is its sequestration. Capturing CO2 from industrial sites and burying it deep under the ground is being approached as a possibility. It can also be transported and sto red deep under the oceans. These methods are still under consideration and some ongoing projects are in their initial stages. However, there are concerns for their success: besides the largescale infrastructure required for such an effort, one of the key concerns is the long-term stability of CO2 underground. Quite possibly, if the sequestered CO2 reservoir breaks and CO2 makes its way to the surface of earth, being heavier than air, it can displace oxygen from the surface.9 There are also concerns over th e impacts to the marine life if CO2 were to be stored deep inside the ocean.10 Research is still going on to make CO2 sequestration a realizable prospect. Current environmental demands, on the one hand, require a limited or diminishing use of the fossil f uels; ever -increasing energy demands on the other, entail development of environmental -friendly new feasible methods of energy production. A swift transition from the current energy system to an efficient and environment friendly one will require multiple technological advances. Some of the possible technologies are nuclear fission, solar, hydro, wind, natural gas, fuel cells. Some of these technologies have been used for a long period of time, while others are relatively new. However, by no measure are any of these capable of carrying the weight of the economy as a whole, some due to lack of availability while others due to their drawbacks. Nuclear fission provides high power density and small waste volume but the potential for groundwater contamination and nuclear accidents is worrisome. Solar photovoltaics have no significant air pollution but can lead to land competition. Hydropower has a sustainable energy supply but it can impact fish habitat and lead to loss of terrestrial habitat. Geothermal energy has no

PAGE 22

22 significant air pollution; however, there can be potential impacts from brine disposal. Only a judicious use of the each will ensure that in pursuit of more energy and reducing greenhouse emission, they do not harm the environment significantly. We ca n say that driven by the huge energy demands, all of the technologies will have to advance simultaneously. 1.4 Challenges in Materials Science Technological advancement in all these areas poses formidable challenges for materials science. If fossil fuel s remain the most important sources of energy before other sources take the forefront, CO2 levels have to be controlled. In industry, the CO2 can be absorbed by materials and released during regeneration so that it can be sequestered underground. Amines ar e one class of promising materials; development of other low cost, more efficient materials is required.11 To take advantage of wind energy, development of stronger and lighter materials for wind turbine blades is required. Currently, the blades are made of fiber composite materials. Off -shore installation requires much stronger materials that can withstand severe environmental conditions. Some of the other challenges are, providing sufficient stiffness to blades to prevent them from deflection, adequate fat igue life under variable wind load and prevention from bucking fracture.12 Similarly, challenges lie ahead for the photovoltaics, which are very promising renewable sources of energy. Increased cell efficiency and reduced manufacturing costs are the foremos t; application of thin films instead of silicon to reduce overall cost. Devices that can work for much larger wavelengths of sunlight spectrum are also needed for the realizing solar power as a dependable source of energy.13 Hydrogen is being looked as anot her major energy carrier in the future. However, to realize a hydrogen economy, better materials for transport through pipes,

PAGE 23

23 storage without loss, easy accessibility, lightweight materials for car -fuel tanks, and efficient conversion of hydrogen into elec tricity are required. Improvement in the storage capacity of the batteries and faster recharging capability are also materials challenges in the battery industry. There are also significant challenges for the fuel -cell industry including high efficiency, l onger lifetime comparable to that of the vehicle, reducing losses. Nuclear fission processes are very good energy producers but challenges lie for the safe and longtime storage of the nuclear wastes, the reuse of the spent fuel, and lower installation cost s. 1.5 Motivation for This Work There is extensive ongoing research in all the above -mentioned fields. (For a review please see reference 14). T he research in this thesis is centered on two energy related fields: solid oxide fuel cells (SOFCs) and nuclear materials. In particular, the focus is on understanding fundamental mechanisms in the evolution of the intrinsic point defects that play key role in materials performance. Defects in materials play a very significant role in determining the properties of the materials. On the one hand, they can be deleterious to the strength of materials, whereas on the other, they can be useful in tailoring the materials properties to desired performance. The advantages of the defects can be observed in the functioning o f the SOFCs. The high diffusivity of the oxygen vacancies (point defects) in the electrolyte material is one of the keys in achieving high ionic conductivity, which is important for good power density. In contrast, due to certain structural (atomistic and electronic) changes, the oxygen diffusivity can be limited. To avoid such structural changes, it is therefore important to understand the fundamental mechanisms that lead to them.

PAGE 24

24 In nuclear reactors, the point defects are highly deleterious as they damag e the materials and its performance. UO2 is used as a fuel in nuclear reactors. It is damaged under extreme environments by the irradiating particles, thus creating point defects. Severe damage due to high concentration of the point defects can limit the t ime of the fuel in the nuclear reactor. In order to enhance the radiation tolerance of UO2, the fundamental mechanism of damage, i.e., point -defect evolution under different conditions needs to be understood. The challenges in materials science lie over v arious length scales; in this thesis, they are attacked from the atomistic scale (using computer simulation). 1.6 Significance of Simulation Atomistic simulations play a significant role in complementing experiments by revealing mechanisms that are eith er interconnected and not easy to disentangle, or are not accessible on the experimental length and time scale. For example, dopants added to stabilize cubic phase of Bi2O3 lower its ionic conductivity, which ultimately reduces power density. The dopants have different atomic radii and ionic polarizability than the host cation (Bi). Experimentally, it is not clear if the conductivity decrease is due to the difference in the radii or the ionic polarizability of the dopant. Using atomistic simulations, arti ficial dopants can be created, in which the ionic polarizability can be varied while keeping the atomic radii constant and vice versa. The dependence of ionic conductivity on polarizability or radii can then be disentangled. In a nuclear reactor, irradiat ion events take place on a very short time scale (~ nanoseconds); capturing them on the experimental time scale is not possible. Atomistic simulations are very helpful in not only capturing those events but are also capable in revealing the fundamental mec hanisms of irradiation -induced processes.

PAGE 25

25 Another exemplification of atomic simulation can be seen in understanding the role of the grain boundaries (GBs). It has been suggested that nanocrystalline materials are more radiation tolerant than their bulk c ounterparts.15 The point defects formed due to irradiation are eliminated by two mechanisms (1) by vacancy -interstitial recombination, (2) elimination at the GBs. Not only do these mechanisms occur a very shot time scale but also take place simultaneously. Using simulations, it is possible to disentangle them and reveal the dominant mechanism between the two. In the following chapters the above examples are discussed in more detail. 1.7 Objective of This Work This work is focused on understanding the evoluti on of the intrinsic point defects in three ceramic materials Bi2O3, UO2 and MgO, and one metallic material Mo. While to some extent, point defects are advantageous in fuel -cell electrolyte materials, they are certainly disadvantageous in nuclear materials. The objective of this work is to understand the fundamental materials responses to point defects, thus revealing the conditions under which materials lose their desired properties. In Bi2O3, while desirable oxygen vacancies act as charge carries and enha nce power density; under certain atomistic environments they themselves act as degradation agents. The environments that have not been understood from experiments until now are unraveled here. In nuclear materials, while the point defects are deleterious, they are not always damaging to materials. Only under specific atomistic environments do materials get damaged, understanding of such environments can help in tailoring better materials. During the course of this research, it was found that conventional ra diation -damage MD simulation is not capable enough to reveal complete evolution of point defects. Hence, a new methodology is developed to perform such simulations. In additions to MD, density

PAGE 26

26 functional theory (DFT) is also used to validate MD simulation results and to reveal certain electronic phenomena that are beyond MD. To characterize materials from structural viewpoint, crystallography is also used to illustrate the crystal structures. 1.8 Organization of Dissertation The dissertation is organized o n the basis of materials. It begins by giving some fundamentals of materials science relevant to this work in Chapter 2. Important aspects of MD simulation are provided in Chapter 3. Chapters 4, 5 and 6 are focused on Bi2O3, where the point -defect evolutio n, oxygen diffusion and their effect on the structure of Bi2O3 are discussed in detail. Chapters 7 and 8 are focused on the fundamentals of radiation damage and point defect evolution in MgO and UO2. Finally, the effect of microstructure on the evolution of point defects in BCC Mo is elucidated in Chapter 9. The conclusions of this work and envisioned future work are given in Chapter 10.

PAGE 27

27 CHAPTER 2 DEFECTS IN SOLIDS 2.1 Introduction On the atomic level, a crystalline material is a collection of atoms tha t are arranged in a periodic fashion providing three -dimensional long -range structural and chemical order. There are a large number of different crystal structures, varying from simple structures of pure metals to more complex ones found in ceramic materia ls. During processing of materials, it is very difficult to maintain an infinite three dimensional long range order. More often, this longrange order is disrupted by defects. These defects can be of different kinds classified by the spatial dimensions: point defects (zero dimensional), dislocations and stacking faults (one dimensional), grain boundaries and surfaces (two dimensional), pores and cracks (three -dimensional). In this work, the effect of point defects on the materials behavior is studied. 2.2 Intrinsic Point Defects Intrinsic Point defect involves one atom, which is either missing from a crystalline site or is additionally present in an open space in a crystal. An atom missing at its crystalline site is called a vacancy ; an additionally present atom is known as an interstitial. Figure 2 1 (a) shows a two dimensional long range lattice of atoms. A missing atom in Figure 2 1 (b) shows a vacancy and an additional atom in the open space represents an interstitial ( Figure 2 1 (c)). In this simple per fect lattice, all nearest neighbor atoms are equidistant from each other, whereas, the lattice is locally strained in the presence of the vacancy or interstitial. The vacancy creates an open space into which the neighboring atoms relax, as shown in Figure 2 1 (b). In contrast, the interstitial

PAGE 28

28 constrains the neighboring atoms because the open space between the crystalline atoms is usually smaller than the size of the interstitial ( Figure 2 1 (c)). Figure 2 1. Point defects. A) 2 -D lattice of atoms. B) a vacancy. C) an interstitial 2.2.1 Frenkel Defect A Frenkel defect is a vacancy interstitial pair formed when an atom jumps from its crystal lattice site to an interstitial site leaving behind a vacancy. A Frenkel defect is shown schematically in Figur e 2 2 (a). It was first observed by a Russian scientist Yakov Frenkel in 1926, and hence the name. Figure 2 2. Point defects in ionic materials. A) Frenkel pair B) Schottky defect Using Krger -Vink notation, the Frenkel defect on the cation (Mg) lat tice can be expressed as: MgMg X Mgi V Mg (2 1) A B C A B

PAGE 29

29 2.2.2 Schottky Defect A Schottky defect occurs when a stoichiometric units of cation and an anion leave the crystal site and segregate at an external surface (grain boundary, surface, etc.) leaving behind a vacancy each on cation and anion site. A Schottky defect is shown in Figure 2 2 (b). This defect is unique to ionic materials and stoichiometric numbers of cations and anions have to be missing to maintain electroneutrality in the m aterial. This defect was first discovered by German scientist Walter H. Schottky. The Schottky defect in MgO crystal can be expressed as: nil V Mg VO (2 2) There are other defects, both intrinsic and extrinsic, which are prevalent in the materials. However, we limit ourselves to the above two point defects that are most central to this thesis. 2.2.3 Thermodynamics of Point Defects: Equilibrium Concentration Point defects are a part of the equilibrium state of a chemically pure ide al crystal. Any material under equilibrium has a specific concentration of point defects. It appears unintuitive because it would cost energy (enthalpy change, H ) to create a point defect, which will increase from the lattice energy of a defect -free crystal. However, the equilibrium state of the crystal is determined not by the lowest lattice energy but by the lowest total free energy. With the creation of a def ect (say, a vacancy), there are entropy changes. The increase in the entropy of the crystal lowers the total free energy of the crystal, thereby, stabilizing the crystal by forming an equilibrium concentration of vacancies. For a given temperature, the equ ilibrium concentration of vacancies can be analytically calculated as shown below.16

PAGE 30

30 Free energy of a perfect crystal The free energy is given as: Gperf Hperf TSperf (2 3) where H is the enthalpy, S is the entropy and T is the absolute temp erature of crystal. The total entropy consists of contributions from two types of entropies, the configurational entropy and the vibrational entropy. S Sconfig ST (2 4) Configurational entropy is due to the different arrangements th at the atoms can have. For a perfect crystal, Sconfig = 0 as there is only one way of arranging N atoms. In Equation 2 4, the vibrational entropy, ST, is the contribution to the entropy due to the vibration of the atoms. In simple terms, it can be given by ST Nk (ln kT h 1 ) (2 5) where N is the number of atoms involved, k is the Boltzmans constant and is the vibrational frequency of the atoms in the perfect crystal. Therefore, putting back to Eq. 2 3, Gperf Hperf NkT (ln kT h 1 ) (2 6) Equation. 2 6 is used later to compare against the Gdef to calculate the equilibrium concentration of vacancies for a given temperature. Now the attention is turned to the crystal with defects to formulate Gdef. Free energy of the crystal with defect s If it takes hd joules to create one defect, and if the number of defects created are nv then, Hdef Hperf hdnv (2 7)

PAGE 31

31 In the defected crystal, since there are nv new crystal lattice sites, therefore, N atoms can now be distributed on N + nv sites. There is a change in the configurational entropy, which is represented as: Sconfig k ( N ln N N nv nvln nvnv N ) (2 8) Due to the presence of the vacancies, the nearest neighbor atoms will essentially vibrate with a different frequency, If the coo rdination number of each vacancy is the total number of atoms affected by the creation of new vacancies is nv. The vibrational entropy of term is given by: S k ( N nv)(ln kT h 1 ) nvk (ln kT h 1 ) (2 9) By combining Equations. 2 7 to 2 9, the free energy of a c rystal containing defects is: Gdef Hperf nvhd T ( Sconfig S ) (2 10) The change in the free energy by creating the defects is given as: G Gdef Gperf (2 11) i.e., G nvhd kTnvln kT ( N ln N N nv nvln nvnv N ) (2 12) At a constant T a G versus nv plot g oes through a minimum which illustrates that the creation of defects lowers the free energy of the system thus making it energetically favorable. Beyond the equilibrium concentration of the defects, further increase in the defects is no longer energeticall y favorable and the free energy increases again. The equilibrium number of vacancies, neq, can be obtained at the minimum of the curve, i.e., when G / nv = 0.16

PAGE 32

32 neqN exp( hd T SvibkT ) (2 13) Hence, the system is thermodynamically driven towards neq under equilibrium. This feature of point defects will be witnessed in Chapter 9 in which the creation of the vacancies from the grain boundaries driving the system towards equilibrium concentration of vacancies will be illustrated. 2.3 Transport 2.3.1 Atomistic Theory of Diffusion Atoms constantly vibrate due the energy provided by the surroundings. At sufficiently high temperatures, the energy is high enough for a vigorous vibration, which can overcome the migration energies (and attempt frequencies) leading to diffusion of atoms. There are various materials phenomenon that occur due to diffusion: carbonization, ionic conduction, creep to name a few. Under a concentration gradient, the randomly diffusing atoms are dir ected in specific direction resulting in a flux. Ficks law gives the relationship between the flux of atoms under a concentration gradient. The flux of atoms, J and the concentration gradient, dc/dx are related by: J Dcx (2 14) where D is the self -diffusion coefficient. The self -diffusion coefficient, D of a diffusing ion or atom is given by:16 D 2 (2 15)

PAGE 33

33 where is the frequency of the successful jumps, is the e lementary jump distance and is a geometric constant. If the diffusion is by a vacancy mechanism, then, =1/ ( is the coordination number of the vacancy). The jump frequency, can be written in terms of the probability of having an energy to make a jump, and the probability that the adjacent site is vacant. That is (2 16) oexp( HmkT ) (2 17) o is the vibration of the atoms. Since a vacancy can make a jump to any of the available sites surr ounding it, Putting all the variables, Dvac 2oexp( HmkT ) (2 18) In atomistic simulations, the self -diffusivity, DA, of the atoms can be directly calculated from the mean square displacement (MSD) of the atoms. DA 1 6 t | ri ro|2 (2 1 9) The vacancy diffusivity Dv and the interstitial diffusivity Di are given by: Dv DAcvfv (2 20) Di DAcifi (2 21) where fv and fi are the correlation factors. Correlation factor depends upon the crystal st ructure and diffusion mechanism.

PAGE 34

34 2.3.2 Diffusion Mechanisms To understand the role of different atomistic species in diffusion, it is important to understand their diffusion mechanisms. Atoms, vacancies and interstitials have different diffusion mechanism s.17 2.3.2.1 The vacancy mechanism An atom diffuses by a vacancy mechanism when it moves into a vacancy occupying a nearest neighbor lattice site. In this mechanism, the atom and the vacancy swap lattice sites; they move in the opposite direction to each other. Figure 2 3 shows the diffusion of an atom through the vacancy mechanism. Figure 2 3. Diffusion by vacancy mechanism. An atom in (A) diffuses in the direction shown by arrow. (B) Atom and vacancy swap positions. By this mechanism, the vacanc y has a higher diffusivity than ion because the ion has to wait for a vacancy to come in its vicinity, whereas, a vacancy can exchange sites with any of the neighboring atoms. 2.3.2.2 The interstitial mechanism In this mechanism, the defect (interstitial) occupies an interstitial site before and after the diffusion jump. Since the diffusing species and the defect are the same, they move in the same direction. Figure 2 4 shows the diffusion through the interstitial A B

PAGE 35

35 mechanism. The interstitial can occupy the lattice site only if it comes across a vacancy (vacancy site). Figure 2 4. Diffusion by interstitial mechanism. An interstitial in (A) diffuses in the direction shown by arrow in (B). 2.3.2.3 The interstitialcy mechanism This is a displacement dif fusion mechanism in which an interstitial atom displaces a lattice atom from the site and occupies it. Each interstitial makes only one jump and has to wait for another interstitial to displace it from it lattice site. Figure 2 5. Diffusion by i nterstitialcy mechanism. An interstitial in (A) knocks a lattice atom and occupies the lattice site in (B). Figure 2 5 shows the interstitialcy mechanism of diffusion. A B A B

PAGE 36

36 2.4 Crystallography of Materials In this section, the crystal structure of the material s under investigation is introduced. In this thesis, while the primary focus is on the evolution of point defects in the fluorite -based materials ( Bi2O3 and UO2), some key point -defect mechanisms are elucidated from other materials (MgO and Mo), which a re used as general model materials. -Bi2O3 and UO2 have fluorite crystal structures, MgO has a rocksalt crystal structure, and Mo has body centered cubic (BCC) structure. 2.4.1 Fluorite Crystal Structure: -Bi2O3 and UO2 In the fluorite crystal structure, the cations form an FCC network with the anions occupying all the eight tetrahedral sites thereby forming a simple cubic lattice (Figure 2 6). Figure 2 6. Fluorite crystal structure. The cations (shown in blue) form a FCC network and anions (in red) occupy the tetrahedral sites forming a simple cubic lattice. There are therefore four cations and eight anions in a fluorite unit cell. The coordination number of the cations is eight and that of anions is four. The space group of fluorite crystal st ructure is Fm 3 m (number 225). The cations occupy the Wyckoff site 4a and the anions occupy 8b.

PAGE 37

37 While UO2 has fluorite crystal structure, -Bi2O3 has a fluorite related crystal structure. In Bi2O3, there are six anions instead of eight and remaining two anion sites are vacant. Due to the vacant sites, the crystal structure is distorted and has the lower symmetry space group, Fm 3 The structural implications will be discussed later in detail in Chapter 6 For now, Bi2O3 is analyzed within the Fm 3 m space group. 2.4.2 Rocksalt Crystal Structure: MgO In the rocksalt crystal structure, both cations and anions form interpenetrating FCC networks ( Figure 2 7). In a rocksalt unit cell, there are equal num ber of cations and anions, therefore, the structure has a stoichiometric ratio of 1:1. The space group of rocksalt crystal structure is same as that of the fluorite, i.e., Fm 3 m Figure 2 7. Rocksalt crystal structure. Both catio ns and anions have a multiplicity of 4, therefore they occupy Wyckoff sites 4a and 4b respectively. 2.4.3 Space Group Fm 3 m All the three materials mentioned above have a space group of Fm 3 m The number of this space gr oup is 225 and the point group symmetry is m 3 m .18 Fm 3 m is more fully

PAGE 38

38 written as F 4 m 3 2 m The first letter, the centering type, represents the translational symmetry in three dimensions. Here, F represen ts a Face Centered Cubic symmetry operation. The three entries following the centering type are three kinds of symmetry directions of the lattice belonging to the space group. According to their position in the sequence, the symmetry directions are referre d to as primary, secondary and tertiary directions. For the cubic lattice, the primary symmetry directions are 100 i.e., [100], [010] and [001]; the secondary symmetry directions are 111 ; and the tertiary symmetry d irections are 110 Figure 2 8. Schematics showing three different symmetry operations represented by the space group F 4 m 3 2 m A) A four -fold rotation symmetry along [100]. B) A three -fold inversion symmetry along [111]. C) A two -fold rotation symmetry along [110]. In the F 4 m 3 2 m symmetry, 4/m represents four -fold rotational symmetry along the primary direction and two perpendicular mirror planes in perpendicular direction to rotation. 3 represents the roto -inversion along the body diagonal of the unit cell which is also the secondary direction. The symmetry operation 2/m represents the two-fold rotation along the edge of the unit cell (tertiary direction) along with two perpe ndicular A B C

PAGE 39

39 mirror planes in perpendicular direction to the rotation. The symmetry operations for four -fold rotation, three -fold roto-inversion and two-fold rotation are given schematically in Figures 2 8 (a), (b) and (c) respectively using a fluorite crystal structure. 2.4.4 Body -Centered Cubic: Space Group Im 3 m Mo has a BCC crystal structure with one atom occupying the center of a cubic unit cell. The schematic representation is given in Figure 2 9 The space group of BCC crystal struct ure is Im 3 m The space group number is 229, the Wyckoff position is 2a and the point group symmetry is also m 3 m .18 In this space group, I represents bodycentered translational symmetry The rest of the symmetry directions are same as the fluorite crystal structure. Figure 2 9. Body-centered crystal structure.

PAGE 40

40 CHAPTER 3 SIMULATION METHODOLOGY 3.1 Simulation Methodologies In principle, all properties of materials are describabl e by using quantum mechanics (QM).19 However, unraveling all materials properties from QM is unrealistic due to small system sizes (102 103 atoms) that QM can handle; materials engineering, in contrast, involve larger length scales upto of the order of 1022 atoms. To elucidate various materials behaviors over wide length scale, different simulation methodologies are used to bridge this gap. At the engineering scale (103 m and above), macroscopic stresses, large temperature gradients, etc., that drive the materials behavior are captured by continuum simulation methods.20 One step down the length scale is the mesoscale (106 103 m); grain boundaries (GBs), voids, etc. that govern the properties of materials at this scale are modeled by kinetic Monte Carlo21 and phase -field methods.22 Further down the length scale lie atomistic scale, modeling of which is done by using molecular dynamics simulation (1010 106 m) using interatomic potentials23, 24 and force fields.25 Atomic diffusion, dislocations and other atomistic phenomenon are accurately captured by this method. At the lowest level lies the electronic scale (1012 1010 m), modeling of which is done by QM methods to elucidate various electronic phenomenon such as atomic bonding and surface adhesion.26 Whil e electronically and atomistically informed multiscale modeling is an area of intense study, individual methods are widely used to understand specific properties in more detail. The work in this thesis focuses primarily on the atomistic scale, hence, mole cular dynamics (MD) simulation is used to understand materials behavior. However, to develop a profound understanding and to validate the hypothesis derived at various

PAGE 41

41 junctures from MD, the electronic -level simulation using density functional theory (DFT ) is also used. 3.2 Molecular Dynamics Simulation MD simulation involves modeling of atoms as point particles, the interactions among which are defined by an empirical potential (plus Coulombic interactions if particles are charged). The dynamics is obtai ned from the Newtons equations of motion. The system of atoms is set up with basic initial conditions like temperature, pressure and number of particles. Each atom is assigned a mass, charge and position, and the evolution of the system is obtained by adv ancing time using an integration scheme. With progressing time, the evolution of the system is followed and various properties including radial distribution functions, dynamics structure factor, phase -diagrams, thermal expansion coefficient, diffusion, pla stic deformation can be calculated. The accuracy of the calculated properties, however, depends entirely upon the fidelity of the interatomic potential used to model the material. Before discussing the potential and its implications, general MD simulation algorithm is discussed. 3.3 General MD Algorithm In MD simulation of a crystalline solid, all atoms are initially assigned with the positions that define the crystal structure of the material. The atoms are also assigned random velocities based on the tem perature given by the expression: 1 2 mv2 3 2 kBT (3 1) where, m and v are the mass and velocity of the atoms, kB is the Boltzmans constant and T is the desired temperature. Before beginning the simulation, the angular momentum is se t to zero.

PAGE 42

42 The forces on the atoms are calculated from the interatomic potential (described later) as: F i U r i (3 2 ) where, U is the potential energy that an atom i experiences due to atom j ri, is the distance between the tw o atoms, and F i is the force. The motion of the atom i, is expressed by the second law of motion where the force, F i is related to the acceleration as F i mi a i (3 3 ) where, mi and ai is t he mass and acceleration of atom i. Equation 3 2 may also be written as mi2 r it2 F i (3 4 ) The velocity and position of atoms at a time t can be calculated from the acceleration of atoms. To advance time by a small increment dt a n integrator is used to obtain the position, velocity and acceleration of the atoms. There are various types of integrating algorithms available including Verlet,27 leap -frog,28, 29 velocity -V erlet, Beemans,30 p redictor -c orrector31 and s ympletic integrators;32, 33 they differ by accuracy and computational load. In this thesis, a 5th-order Gear predictor -corrector method34 is used as an integrating method built in our in -house MD code. This method uses a higher -order Taylor expansion of the atomic position about the current position, x(t) up to the 5th derivative. The positions of the atoms at the next time increment are predicted by using the current position, velocity, acceleration, third, fourth and fifth order derivatives of x(t) given as:

PAGE 43

43 xp( t dt ) x ( t ) v ( t) dt 1 2 a ( t ) dt2 1 3 d3x dt3 ( t ) dt3 1 4 d4 x dt4 ( t ) dt4 1 5 d5 x dt5 ( t ) dt5 (3 5 ) Similarly, velocity, acceleration and higher order derivatives are also predicted. This prediction is not very accurate as it is mere a Taylor expansion. The second step in the predictor -corrector method is the force evaluation step. Here, the force s on the atoms are evaluated using the predicted positions. Using this force, the correct accelerations ac(t+dt) may then be calculated. The error between the predicted acceleration ap(t+dt) and the correct acceleration ac(t+dt) can be given by: a ( t dt ) a c( t dt ) a p( t dt ) (3 6 ) The position, velocity and higher order derivatives are then corrected based on the corrected acceleration. The corrected position, velocity and acceleration is be given by: x c( t dt ) x p( t dt ) co a ( t dt ) (3 7 ) v c( t dt ) v p( t dt ) c1 a ( t dt ) (3 8 ) a c( t dt ) a p( t dt ) c2 a ( t dt ) (3 9 ) where, co, c1, c2 are empirically -determined numerical constants 3/20, 251/360 and 1 respectively. The higher order derivatives are corrected similarly. The predictor step provides an initial guess to which the successive corrective step should converge to correct trajectories. Iterative corrective steps can be employed to obtain a highly accurate answer; in MD however, only one or two corrective iterations are performed to save the computation time. This process is repeatedly followed for a desired amount of time or until the system researches equilibrium or the dynamical behavior of interest has taken place. The time step used for advancing the simulation time is one of the most crucial enti ties. To capture the atomic motion correctly, time step, dt is typically taken to be 1 fs.

PAGE 44

44 It is however necessary to calculate the correct time step from test simulations by checking the conservation of the energy. 3.4 Periodic Boundary Condition There a re two major types of boundary conditions:35 isolated boundary condition (IBC) and periodic boundary condition (PBC). In IBC, the atoms are surrounded by vacuum and are assumed to be far away from any kind of external forces. IBC is generally suited for stu dying clusters of atoms. In contrast, PBC is more prevalent boundary condition for studying bulk materials. In PBC, the cubic super -cell is replicated in three dimensions thus filling the space to form an infinite lattice (Figure 3 1). Figure 3 1. Two dimensional periodic system of a simulation cell. The shaded box is the actual simulation box; surrounding boxes are its periodic images. Also shown is the displacement of atom (marked as 1) by an arrow going out of the box. Due to PBC, its images als o displace in same fashion and another atom enters the simulation box from bottom. (Reproduced from reference 34) By applying PBC, the problem of surfaces is overcome. In a PBC simulation, if an atom moves in the super -cell, its image also moves exactly in the same way in the neighboring super cells. Further, if an atom leaves the super -cell, its image enters from the other side of the super -cell. PBC also allows traveling waves in the system thus allowing the system to behave like a large medium. In this manner, PBCs allows to

PAGE 45

45 model a three -dimensional bulk system. However, if applied incorrectly, an atom can interact with its own image and cause artifact in the results. It is particularly true in the case of defects. To prevent such interaction, a large enough system of atoms is therefore necessary to model materials such that the images are fairly separated and do not interact with each other. (For a detailed discussion on PBC, reader is directed to reference 34 ). 3.5 Interatomic Interactions The origin of interatomic -forces theory can be traced back to the work of Maxwell, Clausius and others (before 1900). It was around that time that an understanding of the at om was being fully accepted, though it was still thought as indivisible. Certainly before that time, the understanding on the origin of the interatomic forces was not prevalent. Clausius introduced the virial theorem in which he showed that the virial coef ficients depend ( B, C, etc.)on the forces that atoms exert on each other.36 pV RT ( 1 B V C V2 ...) (3 10) B is due to collision between pair of atoms, C is due to tertiary collision, etc. The nature of forces was however not known. Van der Waals first showed that the interaction energy between the atoms in a gas consisting of attractive and repulsive potential can be approximated by: ( p a V2 )( V b ) nkT (3 1 1 ) He suggested that atoms could not come close beyond certain distance without overcoming the repulsive forces. It became clearer that at smaller distances, repulsive forces were the dominant forces and an inverse power law could be used to represent the interatomic potential.

PAGE 46

46 U( r ) C1rn C2rm (3 1 2 ) where, C1 and C2 are constants, n and m are integers (n>m) and r is the distance between two atoms. The first term is the repulsive term and second term is the attractive term. (The Lennard-Jones R -6 potential is based on the same form). More broadly though, the origin of the attractive and repulsive forces was still not clear. The fundamental understanding was later developed by the theories from Pauli (1925), London (1930) and others when the quantum nature of particles was recognized. They will be discussed later in some detail. Figure 3 2. Inter -ionic potential between Zr4+ and O2as a function of distance between the two ions. The total potential (ETotal) is a sum of short -range (ES R) and Coulombic (ECoul) potential. The short range potential is mainly repulsive and the Coulombic is attractive. The parameters for the potential are taken from Khan et al.37 The above discussion is based only on the short range interatomic forces, which were then being understood in small molecules, primarily gases like H2, N2 and noble gases. During the same period, theory behind ionicity of materials was also being developed

PAGE 47

47 with the first works done by Madelung and Born.16 It was realized that in the ionic crystals, that atoms were charged and therefore also interacted through Coulombic forces. These interactions are long range and go beyond nearest neighbors. Hence, ionic materials not only have the short range interactions but also long-range Coulombic interactions. To illustrate both types of potentials and their effective sum, the inter -ionic potential energy as a function of distance between Zr4+ and O2ions is shown in Figure 3 2. Since this thesis primarily deals with ionic materials, the origin of both kinds of potentials will be delv ed into. The long range, difficult to calculate, Coulombic interactions are discussed first. 3.5.1 Long-range Interactions The long -range forces arise from the Coulombic interaction of the charged ions, here shown are Zr4+ and O2 -. In ionic system, ions ar e periodically arranged with all cations being at a specific distance from anions. To arrange the atoms in the periodic fashion, the work done to bring ions from infinity to a distance r apart is given by the Coulombic law: ECoul q1q2r (3 13) where, and q1 and q2 are the charges of the ions (here, +4 and 2) and ECoul is the Coulombic energy in eV and r is the inter -ionic spacing in If two ions are oppositely charged, there is an attractive force between them and the Coulombic energ y is negative. If the ions are of the same charge, then they repel each other and the Coulombic energy is positive. The Coulombic interactions are not limited to the closest neighbors but span large distances, as illustrated by the long ECoul tail in Figur e 3 2. To accurately calculate Ecoul experienced by every ion, one has to essentially extend to very long (virtually

PAGE 48

48 infinite) distances, which under is not computationally feasible. As a result, one has to truncate the Coulombic interactions at a certain distance, Rc, which can lead to non convergence (discussed later). To overcome non-convergence, there are methods to accurately calculate ECoul term. Two prevalent methods are the Ewald method34 and direct summation method.38 Ewald summation method is more precise but is also more computationally cumbersome than the direct -summation method. The computational load, in the Ewald summation method, increases as an order of N2 in comparison to N in the direct -summation method. Sinc e throughout this thesis, the focus is on understanding the equilibration of point defects, which occurs on longer MD time scale, a relatively less expensive, direct -summation method is chosen. The method has been shown to work well not only for crystallin e but also disordered systems.38 Direct -summation method The foremost problem that causes non -convergence is that with choice of any truncation radius ( Rc), the system summed over is practically never charge neutral. To overcome this problem, in the direct -summation method,38 charge neutrality is achieved by viewing an ionic crystal system as a charge -neutral molecule consisting of Bravias lattice sites such that the molecules are never broken so as to preserve charge neutrality. For example, instead of considering NaCl as two interpenetrating FCC networks of Na and Cl, it is thought of as (NaCl)4 molecule. When it is truncated at radius Rc, the molecule as a whole is retained. To achieve this in the direct -summation method, the net charge on the truncation sphere is neutralized by adding extra charge on the surface of the truncation sphere of every ion, thus enabling good energy convergence. Lets discuss by taking examples, both, for the problem and its remedy from reference 38.

PAGE 49

49 Without applying any convergence criteria, Figure 3 3 shows a comparison between the energy per ion as a function of cut -off radius, Rc, and a fully -converged value ( Rc ). The energy for cut -off Rc is shown as open circles and the fully converged result is shown as a horizontal broken line. It is clear that the energy deviates considerably from the correct value when the system is left charged after truncation. Figure 3 3. Tota l Coulombic energy per ion obtained by summing the Coulomb pair potential over a cut -off radius Rc (shown in open circles). The straight line indicates the fully converged value ( Rc ). The arrows mark those truncations that have exactly or nearly neutr al systems. (Reproduced from reference 38) Some of the data points that lie on the line (shown by arrows) are those that are summed over neutral or nearly neutral units. The solution proposed in the direct s ummation method is to place neutralizing charges on the surface of the truncation sphere, such that all those ions that are incomplete molecules can form charge -neutral molecules. This is shown in Figure 3 4. The correct energy (EMad: Madelung energy) can then be calculated by subtracting the neutralizing energy from the calculated energy before neutralizing.

PAGE 50

50 Figure 3 4. For every charge qi, an opposite charge qj is placed on the surface of the truncation sphere of radius Rc. (Reproduced from refer ence 38) It can be represented as: Etot Mad( Rc) Etot( Rc) Etot neutr( Rc) (3 1 4 ) 1 2 i 1 Nj i ( rij Rc)qiqjrij 1 2 i 1 Nj 1 ( rij Rc)qiqjRc (3 15) With the neutralizing charges, there is a dramatic effect on the Madelung energy of the sphe rically truncated sphere. A comparison between the charged neutralized and charged truncated sphere is shown in Figure 3 5. The plot shows a significant decrease in the energy fluctuation. However, even with the charge neutralization scheme, the approximat ed Madelung (from equation 3 15) is still not comparable with the exact Madelung energy. A comparison between the two is shown in Figure 3 6, in which the solid circles are same as in Figure 3 5 but on much finer scale. The approximate Madelung energy osci llates about the exact Madelung energy in a slightly damped manner.

PAGE 51

51 Figure 3 5. A comparison between the energies of the charge-neutralized and charged systems. The solid spheres represent charge-neutralized system and the open spheres represen t charge system. (Reproduced from reference 38) Figure 3 6. A comparison of the exact Madelung energy with the approximate Madelung energy. The approximate Madelung energy oscillates about the exact Madelung energy in a slightly damped manner. (Reproduced from reference 38) In order to decrease the fluctuation, a damping function is introduced via the complementary error function. In simple terms, the damped, charge -neutralized coulomb pair potential can be expressed as:

PAGE 52

52 Etot Mad( Rc) 1 2 i 1 Nj i ( rij Rc)( qiqjerfc (rij) rij limrij R c{ qiqjerfc (rij) rij }) ( erfc (Rc) 2 Rc 1 / 2 ) qi 2 i 1 N (3 16) where is the damping parameter. The damped charge neutralized Madelung energy can now be compared with the undamped ch arge -neutralized (Figure 3 6). Such comparison is shown in Figure 3 7. The damped charge -neutralized energy oscillates much as in comparison to the undamped charge -neutralized energy. Further, the energy converges well to the exact Madelung energy within a distance of two lattice parameters. Figure 3 7. A comparison of the damped charge -neutralized Madelung energy with the undamped charge -neutralized energy. The plot also shows the exact Madelung energy represented by broken horizontal line. (Reproduced from reference 38 ) The direct -summation method is built into our in house MD code and has been used throughout the thesis.

PAGE 53

53 3.5.2 Short -Range Interactions Previously, it is discussed that the shor t range interactions are mainly repulsive. However, there are also small attractive forces that act between the atoms. (The short range potential between Zr4+ and O2is shown in Figure 3-2, represented by diamonds in blue) In this section, the origin of both the short range interaction types, including the classic works of Pauli, Debye, Lennard Jones,39, 40 Buckingham,41 Born and London42 will be discussed. From equat ion 3 13, the longrange attractive potential between the oppositely charged ions should bring the two atoms closer to each other such that as the distance between them goes to zero, the atoms should collapse and fuse! However, this does not happen. It follows that there must be some repulsive potential that keeps the two atoms certain distance ap art. To understand the repulsive forces, one has to go back to Paulis exclusion principle,26 which states that no two electrons can have the same quantum numbers. From the orbital viewpoint, there can be only two electrons in one orbital and they have to be of opposite spin. It follows that when two atoms are brought closer to each other, the electron clouds (orbitals) of both atoms penetrate into each other, which creates high density of electrons in a small region, thus, increasing the energy. The charged electron clouds repel each other giving rise to repulsive potential. This is the repulsive potential that is expressed as first term in the Lennard -Jones expression, C/rn. The repulsive energy increases rapidly as the two atoms are brought close to each other. Hence, the interpenetration of the electronic clouds is the reasoning behind the short range repulsive forces.

PAGE 54

54 By using quantum methods to analyze noble gases, Bleick and Mayer36 later gave a more refined form that shows that the repulsive energy follows as exponential form rather than power law. With this in mind, equation 3 11 can now be given as: U( r ) C1er C2r 6 (3 17) where, C1, C2 and are the adjustable param eters. This is the form given in the Buckingham potential widely used in the ionic systems. While the understanding on the repulsive forces was primarily achieved after the advent of the quantum mechanics, some understanding on the attractive part of the Van der Waals forces was established before the development of quantum theory. In 1912 and later, Keesom suggested that the dipole moment, was the most important constant of the forces between two molecules.42 The orient ation of the two dipoles was regarded as the interaction force between the two molecules, and was given as the interpretation for the attractive part in the Van der Waals forces. Later, Debye and Falckenhagen (1920 1922) suggested that apart from orientati on effect, there is also an induction effect; that the charge distribution around a molecule changes in the presence of an electric field (due to polarizability ).42 At this point, London (1930) questioned both concepts and argued from quantum mechanics that rare gases were spherically symmetrical and none showed these interactions. He showed that interaction between the rare molecules could only be explained by the dispersion effect.42 ( For d etailed discussion, please refer reference 42). The sole reasoning was based on the uncertainty principle that a particle cannot be at absolute rest. He considered two spherically symmetrical isotropic harmonic oscillators with no permanent moment, and each with polarizability Using elastic and dipole

PAGE 55

55 interactions energies, and the harmonic oscillator with frequency he showed that the attractive energy between two rare gas molecules is given by: U 3 4 ho2r6 (3 18) According to London, this is the force that illustrates the Van der Waals attractive force between the two uncharged rare atoms, the existence of which is without any permanent dipole or higher multiple preconditioning. The at tractive force here is proportional to r6, which is the second (attractive) term in the Lennard Jones expression, or for that matter, in the Buckingham potential form. Hence, by combining all forces acting on an atom due to the surrounding atoms, the en ergy expression can be given as: U ( r ) q1q2r A exp( r ) C r6 (3 19) where, the first term is the long range Coulombic term, second term is the shortrange repulsive term and third term is the Van der Waals attractive term. A, and C are adjustable par ameters. In this thesis, this potential form has been used to model all ionic materials. 3.5.3 Electronic Polarizability Shell Model The electronic polarizability, is a measure of the effectiveness of an electric field E in polarizing a material. It can be expressed as p E (3 20) where, p is the dipole moment. To account for the electronic polarizability of the atoms in MD simulation, the shell model devised by Dick and Overhauser43 is often used. In the shell model, an atom is described as consisting of a core and a shell. It is analogous to a

PAGE 56

56 real atom, which has a nucleus surrounded by a cloud of electrons. In the shell model, core and shell represent nucleus and inner core electrons, and valence electrons of the atom respectively. The core and shell carry charge X e and Y e respectively; the total charge of the atoms is given by (X+Y) e The shell is attached to the core by a spring of force constant k such that the polarizability of the atom can be given by 1 4 o ( Y2k ) (3 21) or, 14.3994 ( Y2k ) (3 22) where, o is the permittivity of the free space. The units of polarizability are 3. 3.5.4 Thermodynamic Conditions (ensemble) for Materials Simulation In an MD simulation, certain pre -conditions are set to model materials, which replicate the environmental conditions for a real material. That is, the number of atoms, their position and chemical identity, temperature, pressure, volume are pre -defined in simulation to obtain certain physical properties. A system under a fixed number of atoms (N), constant volume (V) and temperature (T), is a simulation under canonical ensemble (NVT). Similarly, fixed number of atoms, constant pressure (P) and constant temperature constitute a n isothermal -isobaric ensemble. There are other system conditions, such as micro -canonical ensemble (NVE), grand -canonical ensemble ( VT) and isobaric isoenthalpic ensemble (NPH) under which systems can be simulated.34 Here, and H represent chemical potential and enthalpy respectively. In this work, all the simulations have been performed using the isothermal isobaric ensemble (NPT), the methodologies for which are discussed in brief.

PAGE 57

57 There are various algorithms that can b e applied to maintain constant pressure.44, 45 In the extended Lagrangian method given by Andersen,44 one ca imagine a set of pistons of certain mass is attached to the system of atoms. The force on the pistons is varied which changes the volume of the simulation box and maintains a constant pressure. The time scale for the box -volume fluctuations is roughly the same as the time for the sound waves to cross the simulation box.44 This method applie d pressure hydrostatically to the simulation cell. Parrinello and Rahman46, 47 later extended this method to allow changes in shape and size of the box. In this thesis, the P -R method is used to simulate the system under NPT conditions. To control the tempera ture of the atoms in the simulation box, a thermostat is used. The working of the thermostat is based on equation (3 1). To achieve a desired temperature, the velocities of the atoms are rescaled such that the total kinetic energy matches the desired tempe rature. One of the simplest of the many thermostats48, 49, 50 is the velocity -rescaling thermostat. In this method, to achieve a desired temperature Td, the rescaled velocity vs of the atoms is given by: vs TdT vi (3 23) where, vi is the initial velocity of the atoms and Ti is the initial temperature. Other widely used methods to control temperature of the system are Nos Hoover,48 51 Berendsen52 and generalized Langevin.50

PAGE 58

58 CH APTER 4 IONIC CONDUCTIVITY O F SOLID OXIDE FUEL CELL ELEC TROLYTE 4.1 Fuel Cell A fuel cell is a device that converts electrochemical energy into electricity. Generally, the sources of electrochemical energy are oxygen/air and hydrogen, or oxygen/air and hydrocarbons. The reaction between the hydrogen and oxygen ions releases electrons, flow of which generates electricity. The main byproducts of the reaction are water and heat. Compared to conventional (non renewable) sources of energy, byproducts of the fuel cell are relatively more environmental friendly. Also to some extent, fuel cells can decrease our dependability on the nonrenewable sources of energy, particularly, coal and oil. There are various kinds of fuel cells, such as phosphoric acid fuel cell ( PAFC), polymer electrolyte membrane fuel cell (PEMFC), molten carbonate fuel cell (MCFC), and solid oxide fuel cell (SOFC). The fuel cells are characterized by the electrolyte used in them. For example, in PAFCs, the electrolyte used is phosphoric acid, in PEMFCs, the electrolyte used in a polymer. Similarly, in MCFCs, the electrolyte used is molten carbonate and in SOFCs, a solid ceramic electrolyte is used. Fuel cells can also be distinguished by the ion that carries the charge. PAFCs and PEMFCs rely on hydrogen ion whereas MCFCs and SOFCs use oxygen transport. In the working of a fuel cell, the transporting ion diffuses through the electrolyte that is sandwiched between a cathode and an anode. If hydrogen is the charge carrier ion, it reacts with oxygen o n the cathode. Otherwise, if the transporting ion is oxygen, it reacts with hydrogen on the anode. In both cases, the reaction between hydrogen and oxygen produces electrons, which travel through the outer circuit producing electricity. The general charact eristics of these fuel

PAGE 59

59 cells are listed in Table 4 1. Since these fuel cells have widely ranging power productions, they have different applications. Broadly speaking, PEMFCs are used for transportation energy generation and small communication devices; MC FCs can be used in marine equipments while PAFCs and SOFCs have anticipated applications in distributed electricity generation. This thesis is focused on SOFCs. Table 4 1. Characteristics of different fuel cells. Fuel Cell Electrolyte Charge Carrier Opera ting Temperature (C) Fuel PAFC Phosphoric Acid H + 100 200 H 2 PEMFC Polymer H + 90 H 2 MCFC Molten Carbonate CO 3 2650 CH 4 SOFC Ceramic O 2 600 800 CH 4 4.2 Solid Oxide Fuel Cell A SOFC is shown schematically in Figure 4 1 At the cathode, oxygen molecules (O2) ionize to form oxygen ions (2O2 -). The oxygen ion transports through the solid electrolyte to the anode where it reacts with hydrogen to produce water and electrons as shown in the schematic. The electrons travel through the outer circuit th us producing electricity. Figure 4 1. Working of a SOFC. (Reproduced from reference 53) Fast transport of the oxygen ion through the high conducting electrolyte is important for achieving high power density in the SOFC. Efficient conductors are Electrolyte

PAGE 60

60 therefo re required to achieve high oxygen conductivity. Doped zirconia is one such high oxygen conducting electrolytes and is widely used in current SOFCs. A comparison of the conductivity temperature plot of yttria -stabilized zirconia (YSZ) with other electrolyt es is shown in Figure 4 2. In addition to high conductivity, it also has good chemical stability in both the oxidizing and reducing environments that commonly occur in SOFCs.54, 55 Hence, for working temperatures around 1000 C, YSZ is a very good SOFC electr olyte. Figure 4 2. Conductivity versus temperature relationship of various electrolytes used in SOFCs. (Adapted from Reference 56, courtesy S. Omar57) However, such high temperatures cause problems including degradation of electrode due to demixin g.58 More importantly, it requires an all -ceramic construction, as metals cannot withstand such high temperatures. This in turn raises cost of the fuel cell. To avoid these problems, intermediate temperatures, i.e., 500 700 C, are therefore the desirable working temperatures of SOFCs. While intermediate temperatures prevent materials degradation, on the downside, the conductivity of the electrolyte also

PAGE 61

61 decreases. Therefore, one of the foremost challenges in SOFC design is to develop materials that can ma intain high conductivity at intermediate temperatures. Doped ceria (gadolinium -doped ceria (GDC) shown in Figure 42) is one such material. However, it has associated problem of low transference number. Transference number, ti = i/(i+e), i and e are i onic and electronic conductivities respectively, i.e., it shows appreciable electronic conductivity, thus decreasing the open circuit potential (OCP) of the cell. Bi2O3 is a material that does not suffer from either of these problems. It has higher conductivity than both YSZ and stabilized ceria (plot for Bi2O3 in Figure 4 2) and also has a high transference number.59 Because of these properties, it is a candidate material for electrolyte in SOFCs. However, it has other associated problems that are dis cussed later. 4.3 Fluorite -Based Cubic Bismuth Oxide 4.3.1 A Model Material Cubic bismuth oxide ( Bi2O3) serves as a model material to understand oxygen diffusion and other atomistic phenomenon in the fluorite -based materials. Fluorite -structure material s, e.g., ZrO2 and CeO2 cannot be directly used as electrolytes in SOFCs. In the operating oxygen partial pressure conditions, since the oxygen conduction takes place by vacancy mechanism, to facilitate conduction, oxygen vacancies are required. In a non -pr imitive fluorite unit cell, for example ZrO2, there are four Zr4+ cations and eight O2anions (Figure 4 3 (a)). As such there are no inherently present oxygen vacancies. Oxygen vacancies are therefore extrinsically created. A standard procedure to create stoichiometric vacancies is by doping with trivalent oxides

PAGE 62

62 (e.g. lanthanides: Y2O3, Er2O3, Gd2O3). An electrochemical reaction for doping of ZrO2 by Y2O3 using Krger -Vink notation can be given as: Y2O3 2 ZrO22 Y Zr 3 OO x VO (4 1) When a trivalent dopant oxide (Y2O3) is introduced into tetravalent oxide (ZrO2), the fewer number of oxygen atoms present in dopant oxide are balanced by oxygen vacancies to maintain stoichiometry. The positively -charged oxygen vacancies mediate in oxygen conduction. However, incorporation of the dopant creates charge and size imbalance on the cation sites. The doped cation sites (negatively charged) and vacant anion sites (positively charged) tend to associate with each other through electrostatic Figure 4 3. Fluorite c rystal structure. A) ZrO2. B) Bi2O3. Ideal fluorite crystal structure in (a) has all eight oxygen sites filled, in contrast to (b) which has two vacant oxygen sites (vacancies). Color scheme in blue are cations, in brown are anions and in green are oxy gen vacancies. Coulombic interactions. Consequently, it leads to formation of associated defect structures such as Y Zr VO Y Zr .60 In contrast, it has been shown recently that oxygen vacancies tend to associate with Zr4+ ions instead of Y3+ ions and is due to larger size of Y3+ ion.61 Regardless of the underlying mechanism, oxygen vacancy association limits oxygen diffusion and hence, decreases conductivity. Hence, although dopants provide vacancies for oxygen conduction, the oxygen conductivity gets limited due to associated B A

PAGE 63

63 phenomenona. Furthermore, dopant cation has different polarizability (usually lower) than host cation that also contributes in decreasing conductivity. The delta phase of bismuth oxide, Bi2O3, is a very convenient model mate rial for doped fluorite electrolytes in that it allows these individual effects to be understood separately. -Bi2O3 has a fluorite based crystal structure but it has six oxygen ions instead of eight (Figure 43(b)). Since Bi is a trivalent cation, for eve ry two Bi cations, only three oxygen ions are required to balance it stoichiometrically. Therefore, the structure can be considered as a fluorite structure with an intrinsic deficiency of two oxygen atoms per unit cell. Although not strictly crystallographically correct since they are an intrinsic part of the structure, it is convenient to refer to the two unoccupied crystallographic sites that would be occupied in fluorite as vacant sites or vacancies. Because these vacancies are present in a stoichiom etric structure, they are charge neutral. The presence of these intrinsic oxygen vacancies in -Bi2O3 means that it has the oxygen vacancy density of a highly doped fluorite structure, however without the dopants themselves. As a result, in pure -Bi2O3 th e effects of dopant cation size and polarizability are absent, allowing the effects of the oxygen vacancies on the structure to be determined in isolation. In addition to being a model for highly-doped fluorite electrolytes, doped Bi2O3 itself has consider able potential as SOFC material. In particular, DyWSB (Dy2O3 and WO3 co doped Bi2O3) has the highest ionic conductivity62 (0.53 S/cm at 500C) among all the fluorite -based materials thus far characterized. 4.3.2 -Bi2O3 Limitations Although Bi2O3 is a promising material, it has still not found wide applications. One of the noticeable limitations from the conductivity plot (Figure 4 2) is its limited

PAGE 64

64 phase stability. Pure Bi2O3 fluorite phase is stable only over the very short temperature range between 825 C and 730 C; it melts beyond 825 C63 and transforms to a monoclinic phase below 730 C.64, 65 For its use in intermediate -temperature SOFC (IT SOFC), melting is not as much a problem as is the phase transformation, upon which, the oxygen c onductivity decreases considerably. Hence, the high conducting fluorite phase has to be stabilized at intermediate temperatures. Lanthanide dopants are added to Bi2O3 to stabilize its high temperature fluorite phase at lower temperatures. However, doping concomitantly decreases its conductivity. A conductivity profile of yttria stabilized Bi2O3 in Figure 4 2 shows lower conductivity than pure Bi2O3. Nevertheless, it is still higher than other materials justifying its usage. 4.3.3 Ionic Radii or Polar izability? Doping decreases the ionic conductivity of Bi2O3. A relative conductivity plot as a function of time (in hours) for different doped Bi2O3 systems is shown in Figure 4 4. The conductivity profile for different dopants illustrates that time ta ken for decrease in conductivity varies for different dopants. While in Yb doped system, the relative conductivity decreases quickly, in Dydoped system, the conductivity decay is less marked. As discussed above, two possible factors for decreased conducti vity are different dopant radius and polarizability from host cation. To elucidate their individual effect, Figure 4 5 shows plots of the conductivity decay time constant as a function of dopant radius (Figure 4 5(a)) and polarizability (Figure 4 5(b)). Ti me constant is defined as the time taken for significant conductivity decrease. Both plots show that the time constant follows a linear dependence on the cation radius and polarizability. Collectively, it means that as the ionic radii and polarizability ar e increased, the conductivity decrease is less

PAGE 65

65 severe. However, the individual effect of either radii or polarizability cannot be elucidated separately as both increase linearly and simultaneously for any given dopant. Figure 4 4. Relative conducti vity of -Bi2O3 doped with different dopants. The time taken for conductivity decrease varies for different dopants. (Reproduced from reference 66) Figure 4 5. Comparison of time constant for different doped -Bi2O3 systems. The dopants are char acterized based on (A) cation radii and (B) polarizability. (Reproduced from references 67 and 66) B A

PAGE 66

66 To engineer better materials for longer -lasting conductivity, understanding the individual effect of either is n ecessary; while experiments are limited on this front, computer simulations prove very useful. 4.4 Simulation Methodology Conventional MD simulation methods as discussed in chapter 3 are used with the interactions between ions described by a combination of electrostatic interactions between the ions, and empirical potentials to describe the short ranged, mainly repulsive, interactions. The form of Jacobs and MacDnaill,68 is adopted in which the interactions are: V ( r ) Z1Z2e2r A exp( r ) B exp( r ) C r6 (4 3) As di scussed in chapter 2, the first term in the potential describes the Coulombic interactions; r is the separation between two ions of charges Z1 and Z2 respectively; e is the charge on an electron. The second and third terms are short ranged repulsive intera ctions, parameterized by A B (energy scales for the repulsive interactions), and and (the range of the repulsive interactions) The last term is the attractive Van der Waals interactions characterized by C The parameters A B and C were fitted to experimental values of pertinent physical properties, including the lattice pa rameter, high frequency dielectric constant for and phases, and the cohesive energies. The details of the potential can be found in Jacobs and Mac Dnaill.68 The potential parameters for the Bi Bi, Bi O, and O O are given in Table 4 2. The polarizability of the ions can be accounted for via the shell model,43 as discussed in Chapter 3; the atom consists of a shell of charge Y|e| which moves with respect to th e massive core of charge X|e|

PAGE 67

67 Table 4 2. Parameters for the short range potential for Bi2O3 from Jacobs and Mac Donaill68 A (eV) () B (eV) ) C (eV/ 6 ) O 2 O 2 1290.0 0.3011 47683.4 0.14438 71.3 Bi 3+ Bi 3+ 159 83.1 0.25643 3087441.4 0.10406 40.36 Bi 3+ O 2 323.0 0.41631 5745.0 0.27279 0.00 The shell and the core are attached by a harmonic spring with force constant k In the presence of an internal or external electric field, the displacement of the shell wit h respect to the core describes the polarizability of the ion, given by equation 322. The values of Y and k are also taken from Jacobs and Mac Dnaill.68 For Bi3+, Y = 30.0 and k =6600.0 eV 2 and for O2-, Y = 3.0 and k =33.0 2. The simulation cell consists of a cubic arrangement of 4x4x4 non-primitive cubic unit cells, each cell containing four bismuth and six oxygen atoms for a total of 640 atoms. Periodic boundary conditions are applied in all three spatia l directions. Before starting the simulation, the oxygen vacancies are randomly distributed on the tetrahedral sites. The time step of 0.05 fs is smaller than is typically used, so as to ensure good energy conservation in microcanonical ensemble simulation s. 4.5 Results 4.5.1 Polarizability Indeed Using shell model43 MD simulation, it is possible to create artificial dopants in which the ionic polarizability can be varied while keeping the ionic radii constant, and vice versa Keeping one constant and varying the other can disentangle the individual effect of either. Figure 4 6 shows a plot of the ionic radii and calculated polarizabilites of corresponding ions. The radii of the ions are taken from Shannon69 and the polarizabil ites

PAGE 68

68 are derived from the relationship between polarizability and ionic radius r (to account for Coulombic screening of outer electrons by 4f electrons) as given by Shirao et al. :70 1.375 3.176 r3 (4 2) Using this relationship, the plot illustrates that the dopants added to Bi2O3 have radius and polarizability considerably less than Bi. Figure 4 6. Ionic radius and polarizability of different dopants added in -Bi2O3. Both radius and polarizability of all dopants are conside rably less than Bi. Figure 4 7. Oxygen ion diffusion through two Bi ions (shown in blue) The displacement of electronic cloud by Bi ions facilitates oxygen diffusion. To understand a physical effect of polarizability, a schematic representation of oxygen ion diffusion through two Bi ions is shown in Figure 4 7.

PAGE 69

69 If Bi ions are polarizable, one can imagine that diffusion of a negatively -charged oxygen ions would be facilitated if the negatively -charged electronic cloud of Bi ions can displace itself away from the oxygen ion. With away electronic displacement, oxygen ion faces less repelling force and hence can easily diffuse through two Bi ions. Therefore, higher ionic polarizability may facilitate higher oxygen diffusion. This is characterized in th e next section. 4.5.1.1 Polarizable and non polarizable ions Using MD simulation, two different systems distinguished by ionic polarizability are modeled. In one system, the ions (Bi and O in Bi2O3) are kept fully polarizable, while in the other, both a re treated as completely non -polarizable. A varying degree of polarizability can by understood from the schematic representation Bi ion in Figure 48. A rigid ion model in Figure 4 8 (a) represents a non -polarizable Bi ion with Bi = 0. Figures 4 8 (b) and (c) represent weakly and strongly polarizable models of Bi ions respectively. Figure 4 8. Gradually increasing Bi polarizability. A) non -polarizable, B) weakly polarizable and C) strongly polarizable. In the presence of electric field, the electro nic cloud of a nonpolarizable model remains un -effected. In contrast, a weakly polarizable and a strongly polarizable models alter their electric fields corresponding to their degree of polarizability. A B C

PAGE 70

70 A fully -polarizable -Bi2O3 model represents real -Bi2O3, whereas, a non polarizable system represents an idealized case of lower polarizable dopant system. In the simulations, t he radii of the ions are kept constant in both systems. Mean square displacement (MSD) of the oxygen ions is used as the physi cal property to distinguish the effect of ion radii or polarizability on oxygen diffusion. MSD is related to diffusivity by the relationship given in equation 2 19. During the simulation, oxygen ions are the only ones that diffuse; since bismuth ions have no vacant sites for diffusion, they vibrate at their respective positions. The MSD profiles of oxygen in two systems are given in Figure 4 9. As expected, a fully-polarizable system (Bi = 2.853 and O = 3.923) allows continuous oxygen diffusion, as give n by a linear MSD profile (Figure 4 9(a)). Figure 4 9. Mean -squared displacement (MSD) of oxygen in polarizable and non polarizable Bi2O3. A) Constant increase in MSD is observed in polarizable Bi2O3 compared with non -polarizable Bi2O3where the diffusion st ops after < 150 ps. B) A magni -polarizable Bi2O3. In contrast, oxygen diffusion is very limited in the nonpolarizable system. The oxygen atoms diffuse only for first few pico -seconds (< 150 ps) during which the syste m equilibrates. After equilibration, the oxygen diffusion essentially ceases as shown by the A B

PAGE 71

71 plateau in the MSD plot (Figure 4 9 (b)). As characterized later in chapter 5, this rapid cessation of oxygen diffusion is a result of the development of a vacancy-ordered structure. By contrast, in the polarizable Bi2O3, there is no such vacancy o rdering on the time scale of MD simulations (0.5 nanoseconds). This qualitative difference in the behavior of the polarizable and non -polarizable systems clearly identifies the polarizability of the ions as being a determ inant of the diffusion behavior. 4.5.1.2 Only -oxygen ions polarizable To further distinguish the effect of the cation and anion polarizability, similar simulations are performed in which, one of the ions (either cation or anion) is treated as polarizable and other is treated as non-polarizable. Figure 4 10. MSD of the oxygen ions in different systems with gradually increasing oxygen polarizability. Bismuth ions are kept nonpolarizable. MSD profile matches with the Figure 4 9 (b) suggesting tha t oxygen polarizability has little effect on oxygen diffusion.

PAGE 72

72 Figure 4 10 shows the MSD of the oxygen ions for the system where bismuth ions are treated as non -polarizable (Bi = 0) and the oxygen polarizability is increased gradually in different syst ems. To put into context, Bi = 0 3 means that the Bi ions are treated as rigid atoms and the electron cloud (or the shell in shell -model) does not respond at all to the electric field. In the same way, O = 0 3 treats oxygen as a rigid atom. O = 3.92 3 is the real oxygen polarizability observed experimentally.71 The MSD profile for various systems shows a similar behavior to that of the non-polarizable system in Figure 4 9 (b). Once the vacancy ordered structure is formed, the oxygen diffusion ceases af ter the system equilibrates and MSD reached a plateau. This illustrates that oxygen polarizability does not control oxygen diffusion. In contrast, when polarizability on ions is reversed, i.e., when bismuth ions are gradually polarized and oxygen ions are treated as non -polarizable, it reveals a different behavior. 4.5.1.3 Only bismuth ions polarizable A gradual increase in the bismuth polarizability gradually increases oxygen diffusion. Figure 4 11 shows oxygen MSD for the gradually increasing bismuth po larizability while keeping zero oxygen polarizability (O = 0) The bottom -most line in Figure 4 11 is the MSD profile for the non-polarized bismuth system, which is same as seen above in Figure 4 9 (b). Gradually increasing the Bi polarizability increases the oxygen diffusion by a corresponding amount. The most polarizable (shown in red) bismuth systems MSD matches with the fully polarizable system in Figure 49 (a). The simulations further indicate that there is a certain minimum ionic polarizability req uired to prevent vacancy ordering. In particular, below a critical polarizability, ordering takes place on an MD accessible time scale, while for higher polarizabilities, there is

PAGE 73

73 considerable diffusion, with no sign of vacancy ordering. This elucidates the reasoning behind the vacancy ordering observed experimentally when Bi2O3 is doped. Hence, this suggests that bismuth (cation) polarizability is most significant to oxygen diffusion. Figure 4 11. MSD of the oxygen ions in different systems wit h gradually increasing bismuth polarizability. Here, oxygen ions are kept nonpolarizable. MSD profile increases gradually as bismuth polarizability is increases. The most polarizable bismuth system matches with Figure 4 9 (a) suggesting that bismuth polar izability is important for high oxygen diffusion. Taken together, simulations elucidate that high polarizable dopants should be considered for phase stabilization to achieve high conductivity. 4.5.2 Oxygen Diffusion Mechanism In this section, the oxygen diffusion in a crystallographic direction in Bi2O3 is discussed. In a Bi2O3 unit cell, to occupy a vacancy, oxygen atoms can diffuse in three different directions, (a) <100>, (b) <110> and (c) <111>. They are schematically shown in Figure 4 12. In Figure 4 12, any of the six oxygen atoms having sufficient activation

PAGE 74

74 energy can make a jump to occupy either of the two available vacancies. In a unit cell, <100> is the shortest diffusion path (=0.5ao), <110> has intermediate diffusion path length along the face diagonal of oxygen sublattice ( 2 / 2 ao) and <111> has the longest path length along the body diagonal of oxygen sublattice ( 3 / 2 ao). Figure 4 12. Possible oxygen diffusion directions. A) <100>, B) <110> and C) <111>. Using MD simulations, the path of any migrating atom can be followed and the prominent diffusion direction can be statistically calculated. From multiple simulations, it is found that oxygen ion prefers to diffuse along <100> diffusion direction. A time averaged plot of Bi and O atoms in the unit cell is shown in Figure 4 13. A three dimensional and a two -dimensional view of the unit cell are shown Figures 4 13(a) and 4 13 (b) respectively. The time averaged path of all Bi atoms is shown in black color and all other colors represent oxygen atoms. The paths of all Bi atoms are limited around ideal Bi positions, thus illustrating only small atomic vibration with no diffusion jump. Similar to Bi atoms, the time averaged path of five oxygen atoms, i.e ., two shown in dark blue, one each in light blue and green, and one of the two in red show vibration at their corresponding sites. These five oxygen atoms do not diffuse for the MD time represented in this plot. The sixth oxygen atom in the unit cell unde rgoes two diffusion jumps marked as 1 and 2 and shown in pink and red lines respectively. In both diffusion processes, A B C

PAGE 75

75 oxygen jumps in <100> direction. Twodimensional view of these diffusion processes clearly show <100> jumps. Figure 4 13. Oxygen diffusion mechanism. A) A 3 D representative view of time averaged position of Bi and O atoms in a unit cell. Bi positions are shown in black and oxygen in other colors. Out of six oxygen atoms, five vibrate at their corresponding positions. One oxygen at om diffuses in <100> direction. It makes two jumps in the diffusion path shown by lines in pink and red. B) A 2 D plot of same oxygen atom making two jumps. 4.6 Discussion According to Boyapati et al.,72 doped -Bi2O3 forms an ordered structure below 600C in which the vacancies tend to order in <111> (occupational ordering), with concomitant oxygen displacement (positional ordering). The occupational and the positional ordering increase as smaller and less polarizable dopants (e.g., Yb3+ and Er3+) are adde d, which result in a decrease in the oxygen conductivity. These dopants have been observed to lead to shorter time constants for oxygen ordering than the time constants associated with bigger and more polarizable dopants, such as Dy3+. In addition, for Dy3 + doped systems Boyapati et al. did not observe complete ordering with ageing. This means that even at relatively low temperatures, the formation of vacancy ordered structure can be prevented by the presence of larger and more polarizable dopants. While t he system simulated is different in that there are no dopants, the simulation results are B A 1 2

PAGE 76

76 consistent with these experiments. The simulations show that polarizability is much more important for high conductivity. Further, the simulations for the fully-polar izable system show no indication of either kinds of ordering. In contrast, the non-polarizable system (the extreme limit of low polarizability; low polarizability ions Yb3+ and Er3+) produce both occupational and positional ordering in doped Bi2O3. Hence, MD simulation elucidates the reasoning behind loss of conductivity with time. The occupational ordering will be discussed (vacancy ordering) in Chapter 5. The implication of vacancy ordering on the structure will be discussed in Chapter 6 when the ions di splacement from crystallographic viewpoint is discussed. This constitutes positional ordering. Oxygen diffusion direction was previously observed experimentally by Boyapati et. al .73 Similar to observation made here, they also observed <100> direction as p rominent oxygen diffusion direction. They further showed that exact oxygen path is not along straight <100> direction but is curved, i.e., the oxygen atom diffuses through the empty octahedral position in the center of the unit cell. However, results obser ved here are partly in agreement. The preferential diffusion is along <100> essentially along straight line with no indication of oxygen passing through octahedral site. It is important to mention there that this diffusion illustration is for the non -polar izable system (Figure 4 8(b)). The exact oxygen diffusion path in the simulated polarizable system could differ, i.e., might occur through octahedral site, as was the case in the systems considered by Boyapati et al.

PAGE 77

77 CHAPTER 5 VACANCY ORDERED STRUCTURE O F CUBIC BISMUTH OXID E 5.1 Introduction In the previous chapter, using MD simulation, it is found that the non polarizability of ions leads to oxygen diffusion cessation. In this chapter, the implications of cessation on the structure of Bi2O3 is discussed. Evidently, the cessation is related to the ordering of the vacancies intrinsically present in Bi2O3. On the one hand, diffusing oxygen vacancies contribute to high oxygen diffusivity; however, on the other, their ordering limits it. In this chapter, the focus is on the non-polarized ionic system of Bi2O3 to elucidate the vacancy -ordering mechanism, and the structural changes that occur in Bi2O3. Up to this point, MD is used to analyze the system. It has allowed to observe oxygen diffusion at high temperatures in a large system of atoms. However, it is well known that MD methods are intrinsically limited in their materials fidelity. By contrast, electronic structure methods, particularly density -functional theory (DFT) approaches a re known to give a much better description of materials structure. Hence, in this chapter, both MD and DFT are used to provide further corroboration for the structure determined from MD simulations and experiments. A non -primitive unit cell of -Bi2O3 is shown in Figure 4 3 (b). Similar to fluorite crystal structure, bismuth ions form a FCC network and oxygen ions sit at the tetrahedral sites. As also mentioned previously, Bi2O3 has six oxygen ions instead of eight, leaving two vacant sites (vacancies). These vacancies can be at any of the eight available anion sites. If the vacancies were to form a simple ordered array, it would seem likely that they would order along <100>, <110> or <111> directions, i.e., <100> as first nearest

PAGE 78

78 neighbors, <1 10> as second nearest neighbors, and <111> as third nearest neighbors. The structure of Bi2O3 has therefore been a subject of considerable discussion, and various models have been put forward with particular interest focusing on possible ordering of vacancies.63 68 74, 75, 76 The Gattow model, posits equal possibility of tetrahedral site occupancy by oxygen ions, i.e., the oxygen sublattice is completely disordered and there is no preferential vacancy ordering (Figure 5 1 (a)). In contrast, according to the Sillen model the oxygen sub lattice is completely ordered, i.e., vacancies are preferentially ordered in the <111> directions (Figure 5 1 (b)). These models are in complete disagreement to each other. Furthermore, they focus only on the site occupancy and do not discuss the possible oxygen displacements. The Willis model74 accounts for such displacements that each oxygen ion is displaced from regular tetrahedral 8c site in <111> direction towards the face of the tetrahedron (Figure 5 1 (c)). In the unit cell, there are eight tetrahedral positions and each tetrahedron has four faces. There are therefore 32 equivalent positions (32f) to which oxygen ions can displace. Since, there are only six oxygen ions in the unit cell, it gives occ upancy of 3/16. In this manner, the Willis model accounts for the oxygen displacement but lacks an accounting for the vacancy occupancy. Hence, in Figure 5 1 (c), only the displacement of oxygen in <111> direction towards one of the tetrahedron faces is sh own, and not the positioning of vacancies. Hence, while the Willis model accounts for positional ordering (oxygen displacement), Gattow and Sillen models account for occupational ordering (vacancy ordering). Individually however, none of the three models a re self -contained in describing the structure completely. Furthermore, these models do not discuss about the cation displacements.

PAGE 79

79 Figure 5 1. Three different models of Bi2O3. A) Gattow model all sites have equal occupancy (in the schem atic no vacancies are shown; equal probability is shown by oxygen in lighter color). B) Sillen model vacancies order in <111> direction. C) Willis model Oxygen displaces in <111> direction towards the face of the Bi tetrahedron. In (C), only the displ acement of oxygen atoms is shown. Willis model does not account for the positioning of vacancies. Recently, from analysis of neutron-diffraction experiments, Boyapati et al .72 concluded that the vacancies in doped-Bi2O3 o rder along <111> with additional ordering in <110> direction, thus accounting for occupational ordering. Moreover, neutron scattering showed that the oxygen ions are displaced from the 8c sites, with the data refinement being consistent with displacement t o the 32f sites.72 This accounts for positional ordering. However, in contrast to these experiments, Walsh et al.76 more recently reported from DFT calculations that the vacan cies prefer nearest neighboring, B A C

PAGE 80

80 i.e., <100> ordering; they did not however explore the structure proposed by Boyapati et al These results might not necessarily contradict each other because the experiments were performed on doped systems, while the DFT c alculations were performed on pure Bi2O3. However, by exploring only three possibilities, i.e., <100>, <110> and <111> and not the one observed by Boyapati et al. it appears that although <100> could be energetically the lowest among the three, but perh aps not the lowest structure of all. This inconsistency is explored by comparing the experimentally consistent MD results with DFT calculations. It is found that the vacancy ordering in combined <110> and <111> predicted by experiments and MD simulation i s consistent with DFT calculations (shown later). Moreover, the results in the previous DFT calculations by Walsh et al. were due to small system size (only one unit cell) considered for calculations. It is found that diffusion cessation in non -polarizable ionic system forms combined vacancy ordering and it is this ordering that is observed in experiments consistent with DFT. Therefore, in this chapter, use both MD simulations and DFT Bi2O3. 5.2 Ordering of the Vacancies: MD Simulation Recalling the MD simulation methodology, before starting the simulation, a 4x4x4 unit cell system was generated by randomly placing the vacancies with no site preference. All the atoms were also placed at the crystallographic positions. This structure was then equilibrated at high temperature to allow oxygen diffusion. After only few oxygen jumps, the structure thus obtained after oxygen diffusion (in the tail end of Figure 49(b)) is the one that is ana lyzed here.

PAGE 81

81 5.2.1 Neighboring V acancies of a V acancy For the purpose of this analysis, the vacancies are assumed to sit at the crystallographic site (8c) of the missing oxygen atom. The vacancy sublattice is then defined by the pattern of ordering of the vacancies. To understand the nature of this ordering, the focus is on the vacancyvacancy neighbors. A pair distribution function (PDF), averaged over all vacancies, for the neighboring vacancies of a vacancy is given in Figure 5 2. Because the vacancies are randomly placed in the initial structure (non equilibrated structure) the resulting PDF in Figure 52(a) is exactly what one would expect from a statistical analysis of such a random structure, with peaks at all the cubic neighbor distances of 0.5, 0.707, 0.866, 1.0a0 etc., with heights in accord with statistical expectations. Figure 5 2. Vacancy -vacancy pair -distribution function. A) Non-equilibrated structure. B) equilibrated structure. As is evident from Figure 5 2(b), the equilibrate d structure is very different. The PDF of each individual vacancy is identical to that shown in Figure 5 2(b), demonstrating that all A B

PAGE 82

82 vacancies are crystallographically identical to each other. Particularly important is the absence of peaks at 0.5a0 and 1a0; this indicates that there are no vacancy -vacancy nearest neighbors along <001>. On the other hand, every vacancy has three vacancy neighbors in the <110> directions. These <110> vacancies form an equilateral triangle with the central vacancy (CV) at the center, as shown in Figure 5 3. In addition, there are two vacancies in the <111> directions normal to the plane of <110> vacancies, and in opposite directions relative to the CV (Figure 5 Bi2O3 has a combined ordering of the vacancies in <110> and <111> directions. Figure 5 3. A central vacancy (CV) is surrounded by three <110> vacancies and two <111> vacancies. This is a repeating configuration of the vacancies throughout the system. This vacancy ordering mechanism is consistent with the occupational ordering observed in neutron scattering experiments73 Bi2O3. Recalling the result from Walsh et al. where <100> was suggested to be the prominent vacancyordering mechanism, the high temperature equilibration MD simulations show that in fact <100> does not exist at all. Instead a combination of ot her two ordering directions is prevalent. 120

PAGE 83

83 5.2.2 Continuous Ordered Structure of Vacancies Neighboring vacancies do not form clusters but a continuous network throughout the system both in <111> and in <110> directions. To represent the continuous network of <110> vacancies, the oxygen sublattice is focused. Since, there are three <110> vacancy neighbors of a vacancy, continuous chains of vacancies are present in all three <110> directions. Two such chains are shown in Figure 5 4(a). Two of the three <110> vacancy neighbors of every vacancy are shown here; the third neighbor is in the direction normal to the plane of the figure. The network formed by the <111> ordered vacancies is shown in Figure 5 4(b). Figure 5 4. A continuous network of v acancies in the oxygen sublattice. A) Two parallel chains of <110> network of vacancies (shown in orange) connected by <111> ordered vacancies (shown in blue). B) <111> network of vacancies creating an open channel through the system. A B z x y

PAGE 84

84 This continuous net work of vacancies in <111> is simultaneously present in all the four families of <111> directions. The <111> network of vacancies acts as a connection between two parallel chains of <110>. As a result, one ordering network is directly connected to another. It is important to recognize that although there are open channels of vacancies available for easy oxygen diffusion along the <111> network, the oxygen diffusion is limited by preferential oxygen diffusion in <100> directions. Because a combined ordering in <110> and <111> is not possible in a single unit cell, a 2x2x2 system of unit -cells is required to observe the complete ordering phenomenon. This is the reason that while working with only one unit cell, Walsh et al. did not observe this combined orderi ng mechanism. Figure 5 5. A 2x2x2 fluorite superstructure with a combined vacancy ordering in <110> and <111> directions. A) Simulation. B) Experiment (reproduced from reference 75). The vacancy or dering mechanism from MD simulation agrees with experiment. A B

PAGE 85

85 The 2x2x2 fluorite related superstructure determined from simulations is shown in Figure 5 5; this superstructure is completely consistent with the structure determined from the analysis of electr on75 and neutron diffraction experiments.72 5.2.3 Expansion and Contraction of the Tetrahedron Although there are numerous ways in which the structure of -Bi2O3 can be visua lized, it is particularly instructive to consider the tetrahedron formed by four bismuth atoms. Each of these tetrahedrons encloses either an oxygen atom or a vacancy. Every Bi atom sits at the corner of eight tetrahedron, six of which contain oxygen atoms two of which contain vacancies. The Bi -Bi inter atomic distances are different in the two cases. Figure 5 6 shows the Bi -Bi PDF for the equilibrated system. Figure 5 6. Bi Bi pair distribution function. Two different first nearest neighbors dis tances are represented by the first two peaks of equal height. There are two distinct nearest -neighbor distances, corresponding to Bi Bi separations of 0.67ao and 0.76ao in comparison to one 0.707ao. If there is an oxygen atom inside the tetrahedron (Figure 5 7(a)) then four out of the six Bi Bi inter atomic distances are short (0.67a0) while two are long (0.76a0). The short bonds are shown by full lines

PAGE 86

86 whereas the long bonds are shown by broken lines. On the other hand, if there is a vacancy inside the t etrahedron (Figure 5 7(b)), all the six Bi -Bi inter atomic distances are long. The volume of the empty tetrahedron is thus larger than that of the full tetrahedron; this is physically very reasonable since the oxygen ion has strong attractive Coulombic int eractions with the Bi ions surrounding it. When an oxygen atom and the vacancy are juxtaposed, the two tetrahedrons share a longer Bi -Bi distance between them (Figure 5 7(c)). Figure 5 7. Tetrahedron expansion and contraction due to a vacancy and an oxygen atom, respectively. A) Oxygen present inside the tetrahedron. Oxygen tetrahedron has four short and two long Bi -Bi interatomic distances. B) Vacancy inside the tetrahedron. Vacancy tetrahedron has all six long Bi Bi interatomic distances. C) A vacancy tetrahedron juxtaposed with an oxygen tetrahedron share a long Bi Bi interatomic distance. 5.3 Density Functional Theory (DFT) Calculations of the Structure of Bi2O3 In this section DFT is used to provide further corroboration for the structure det ermined from experiment and from the MD simulations. 5.3.1 DFT Methodology Vienna Ab initio Simulation Package (VASP) is used for these calculations.77, 78, 79 In particular, the GGA exchange -correlation functional in PW9180 form is used. The projector augmented wave method (PAW)81 was applied to the core electrons (Bi A B C

PAGE 87

87 [Xe]4f145d10), O[He]). A plane wave cutoff of 500eV was used with a k -points grid of 6*6*6 for systems containing 1x1x1 non-cubic unit cells and 3*3*3 k -points grid for 2x2x2 systems. The positions of the atoms were optimized such that the force on each atom had converged to less than 0.0005 eV/ Pe riodic boundary conditions were applied in all the three directions. 5.3.2 V acancy Ordering Walsh et al.76 have recently used DFT to examine defect -ordered structures in Bi2O3. As a check on the implementation, these calculations are repeated with vacancies aligned in the <100>, <110> and <111> (see Figure 58). These calculations require a simulation cell containing only a unit -cell of two Bi2O3 stoichiometric units, i. e., a total of 10 atoms. The geometry optimization was performed using the experimentally obtained lattice parameter of 5.644 .82 Figure 5 8. Possible vacancy Bi2O3. A) <100> B) <110> and C) <111>. The lattice parameter, the a/c value, cell volume and the relative energies were found to be in quantitative agreement with the previously published results. In particular, in agreement with the Walsh et al. ,76 it is found that <10 0> and <110> vacancy -ordered structures show a tetragonal distortion whereas the <111> vacancy-ordered structure A B C

PAGE 88

88 maintains cubic symmetry. The relative energy differences are also consistent, predicting <100> ordered as the most stable structure. However, this result is in disagreement with the earlier analyses both from theory68, 83 and experiments.75 84 This work has extended the DFT work by examining a simulation cell containing 2x2x2 unit cells. Simulations, with <100>, <110> or <111> ordering showed the same energy trends as for the 1x1x1 systems. Table 5 1. Relative stability of the ordered structures in 2x2x2 system. Ordering Relative Energy (eV/Bi) <100> 0. 000 <110> 0.066 <111> 0.542 <110> <111> 0.225 The energy of the structure predicted by the MD simulations and obtained from the analysis of the experiments, i.e., the network of <110> and <111> ordering is also calculated. As Table 5 1 shows, this c ombined ordering in <110> and <111> is significantly lower in energy than any of the three previously considered structures. 5.3.3 Bonding and Charge Distribution To characterize the bonding in the system, the electron density in the structure is examin ed. The electron density contour maps for the individually ordered structures in <100>, <110> and <111> were found to be consistent with those shown by Walsh et al.76 It is useful to examine the structure fro m the perspective of the Bi ion, which is in an octahedral environment with respect to the oxygenvacancy sub -lattice. Walsh et al. found that when the vacancies are in a <100> or <110> arrangement, the Bi electron lone -pair charge is directed towards the vacancies. By contrast, when the vacancies are in a <111> arrangement, no lone pair is formed and the charge is essentially evenly distributed around the Bi ion.

PAGE 89

89 Figure 5 9. Two kinds of Bi -centered oxygen sublattice cubes are present in the co mbined vacancy-ordered structure. A) Bi surrounded by <110> vacancy ordering. B) Bi surrounded by <111> vacancy ordering. The distribution of the electron lone -pair charge in the combined vacancy ordered structure, i.e., <110> and <111>, maintains an analogous charge distribution to that observed in the individual vacancy ordered structures. Table 5 2. Bi O bond length in Bi -centered <110> and <111> vacancy-ordered oxygen sub -lattice <110> <111> Bi O () Bi O5, Bi O6 2.303 2.426 Bi O1, Bi O2 2.399 Bi O3, Bi O4 2.559 However, in the combined vacancy -ordered structure, there are two different environments around the Bi atoms, with some Bi atoms surrounded by <110> ordered vacancies and some by <111> ordered vacancies (i.e., in the arra ngements shown in Figures 5 9(a) and 5 9 (b)). As a result, there are varying bond lengths (Table 5 2) and different electron lone -pair charge distributions ( Figure 5 10) The charge distribution is represented in two (110) planes perpendicular to each other both passing through the central bismuth ion. All of the Bi O bonds are of equal length as shown in Table 5 2. A B

PAGE 90

90 Figure 5 10. Electron lone pair charge around two different Bi -centered cubes. One of the perpendicular (110) pla nes in the oxygen sublattice is shown for both <111> and <110> vacancyordered Bi -centered cubes. A) In <111> vacancy ordered oxygen sub lattice, the electron lone -pair charge is equally distributed. B) All four Bi -O bond lengths are same. C) <110> vacanc y -ordered oxygen sublattice; electron lone -pair charge is directed toward the vacancies. D) Shorter and longer pairs of bond lengths predicting the distorted oxygen sub lattice with an evidence of lesser bonding between Bi O3 and Bi O4. (Contour levels ar e between 0.0 e/3 (red) and 0.4 e/3 (pink). Due to the symmetric ordering of vacancies in <111> cubic sub-lattice (Figures 5 10(a) and 5 10(b)), the electron lone -pair charge is equally distributed to form the Bi O A B C D

PAGE 91

91 bonds. In contrast, the Bi centered <110 > ordered oxygen sublattice has different bond lengths, which occur in pairs. The two perpendicular (110) planes for <110> ordered sublattice are shown in Figure 5 10(c) and 510(d). The plane containing the two vacancies here (in Figure 5 10(c)) has a Bi O bond length that is the same as in the <111> ordered oxygen sub lattice. However, in the perpendicular plane (Figure 510(d)), two of the Bi O bonds are shorter and two are longer thus distorting the Bi centered <110> oxygen sublattice. As expected, t here is a directional distribution of the electron lone -pair charge towards the vacancies. In this distorted sub lattice, there is further an evidence of weaker bonding between the Bi O5 and Bi O6. 5.3.4 Electronic Structure of Vacancy -Ordered -Bi2O3 Due to two different vacancy ordered environments around the Bi ions, there are different Bi O bond lengths. The ionicity of the ions and the surrounding atomic environment affects the bond length. The charge distributions associated with the Bi and O ions ar e characterized using Bader charge analysis85, 86 of the electronic structure determined from DFT (Table 5 3). If Bi2O3 were fully ionic, the Bi and O ions would have formal charges of +3e and 2e respectively. Instea d, Bader charge analysis shows that the Bi ions inside <111> ordered vacancies have a DFT atomic charge of +2.1232e while the Bi ions inside <110> ordered vacancies have a charge of +2.8411e. Table 5 3. Electronic charge on the ionic species. Vacancy or dering environment Ion Charge (e) <111> Bi +2.1232 <110> Bi +2.8411 O 1.77

PAGE 92

92 While Bi ions with charge +2.84 e is very similar to the valence charge of +3 and show ionic characterm, the Bi ions with charge +2.12 e is significantly deviated and hence sh ows a degree of covalency. The charge on the O ions is 1.77e. Numerically, the deviation of the DFT atomic charge from the formal charge gives the amount of charge associated with covalent bonding. It is useful to characterize the charge distribution fro m the perspective of the Bi tetrahedra (Figure 5 11(a)). As shown in Figure 5 11(b), each oxygen ion present inside the tetrahedron is bonded to three Bi ions of higher ionic charge +2.84e, and one of more covalent character: charge +2.12e. Figu re 5 11. A) A tetrahedra in the unit cell. B) Each b ismuth tetrahedra is formed by three Bi (surrounded by <110> ordered vacancies of charge 2.84e and one Bi (surrounded by <111> ordered vacancies) of charge 2.12 e. All four Bi O bonds are of different lengths. The four different bond lengths observed previously in Table 5 2 are shown here. All four Bi -O bonds inside the tetrahedra have different lengths. T h e tetrahedral configuration of Bi ions is same even around the vacancies. A B

PAGE 93

93 5.4 Discussion According to Larrif and Theobald,87 oxygen diffusion in Bi2O3 occurs in the direction of the lone pair of electrons. DFT calculations in this thesis predict the existence of a directional lone pair of electrons when Bi has vacancy neighbors in the <110> direction, and an isotropic electron distributions in Bi has vacancy neighbors in the <111> direction. It appears that the electronic charge distribution plays a significant role in limiting the diffusion. Since a covalent bond is directional and difficult to break, it is speculated that the covalent Bi O bond is the limiting factor in the oxygen diffusion in the vacancy ordered system. As such, it is anticipated that variations in the covalent Bi O bond characteristics might be the key crystallochemical variable chi efly responsible for tailoring the ionic conductivity of bismuth oxide based electrolytes. A corresponding argument may well be applicable to ZrO2 and CeO2 based electrolyte materials.

PAGE 94

94 CHAPTER 6 CRYSTALLOGRAPHIC ANA LYSIS OF THE STRUCTU RE OF CUBIC BISMUTH OXIDE 6.1 Introduction In the previous chapter, it is found that a 2x2x2 superstructure is required to completely define the vacancy ordering mechanism (occupational ordering) in Bi2O3. Due to the ordering, it is observed that there are two different kinds of Bi ions surrounded by different vacancy ordering environments for which the bonding environments from charge -distribution viewpoint is established. In this chapter, the stru ctural distortion that occurs due to the vacancy ordering. i.e., the positional ordering (displacement) of the oxygen ions and the cation displacements (which have until now remained subtle in the wide literature on Bi2O3) is established. A systematic cr ystallographic analysis of the vacancy Bi2O3 is also presented. This Chapter begins by presenting a schematic of the 2x2x2 super -structure Bi2O3, which due to the vacancy ordering differs from t he conventional fluorite space group. This will lead to the discussion on the consequent ionic displacements and structural distortion in the sub -structures. 6.2 Schematic Representation of the Structure of -Bi2O3 As a way to begin describing the 2x2x2 s uperstructure which contains 32 Bi atoms and 48 O atoms, the plane views of the structure obtained from DFT simulation are presented along the (100), (010) and (001) in Figure 6 1. All the columns of cations contain cations that are displaced from the colu mn axis. However, there are certain cations, which do not move with respect to the positions of an ideal double fluorite structure. For example, cations #s 1, 5, 9 and the cations behind them in Figure 61(a) do not displace when viewed from either of th e planar directions. On the other hand,

PAGE 95

95 some cations such as, #s 3 and 7 which are not displaced in the (100) plane (Figure 6 1(a)) are in fact displaced along (001), as can be seen from the other perspectives. Figure 6 1. The 2x2x2 fluorite -related superstructure viewed from three different planes. A) [100], B) [010] and C) [001]. Bi ions are shown in gold and oxygen ions are shown in red. Figure 6 -Bi2O3 along [001] direction showing the relative displacement of the bismuth ions with respect to the original fluorite positions and the distortion of the oxygen sublattice. Furthermore, it is found that the cations move in pairs in different directions. For instance, when viewed from one of the planes (e.g., Figure 6 1(a)), a (010) displacement A B C

PAGE 96

96 of the cations (atoms # 2 & 27 or 12 & 37) is followed by a displacement in (001) (e.g., atoms # 4 & 29 or 39 & 14). In Figure 6 2, a layer sequence of the structure along the <001> direction is presented using the same convention as Galasso88 when describing the pyrochlore structure (another fluorite -related 2x2x2 superstructure). Schematic shows relative displac ement of the bismuth ions relative to the exact fluorite positions (shown in grey). The oxygen ions are also displaced from their respective fluorite positions. The layered positioning of the oxygen and vacancies are shown in Figure 6 3. The oxygen displac ements are not shown here. Figure 6 3. Layered sequence of the oxygen sub-lattice along [001] in 2x2x2 Bi2O3. The displacements of the ions are not shown. Oxygen ions are shown in red and vacancies are shown in white. The above schematics, Figures 6 1 and 6 2, show exact positions of the cations and anions in the 2x2x2 fluorite -related structure of -Bi2O3. Figure 6 3 shows the arrangement of vacancies in different planes. These schematics give an overview of the Bi2O3. Due to the vacancies ordering and ionic displacements, the structure differs from the ideal fluorite space group.

PAGE 97

97 6.2.1 S pace Group Analysis In agreement with experiments, the DFT simulations predict a defect -derivative fluorite crystal structure with ordering in the unoccupied anion sites in the <110> and <111> directions. The space group of an ideal fluorite cubic struct ure is Fm 3 -Bi2O3, the overall cubic symmetry is maintained. However, the lattice parameter is very close to twice that of normal fluorite structures: 11.19761 as determined from the DFT calculations. The unit cell contains 16 Bi2O3 stoichiometric units rather than two as in the fluorite unit cell. The experimentally reported lattice parameter of a unit Bi2O3 is 5.644 ,82 which corresponds to 11.288 for the 2 2 2 superstructure. This DFT error of 0.8% may be taken as representative of the degree of quantitative accuracy that can be expected in the analyses below. The superstructure lowers the symmetry to space group of Fm 3 (202). It is important to recall that Fm 3 is a maximal non -isomorphic t subgroup of Fm3m, where all translational symmetry is maintained but the order of the point group is reduced (from 48 to 24) .18 In essence, Fm3 is derived from the Fm3m in that the 4 -fold rotation along <100> is reduced to a 2 -fold rotation, and that the 2 -fold rotation along <110> and the associated perpendicular mirror planes are lost. The reduction in the 4 -fold rotation to 2 -fold rotation is due to the displacement of the Bi ions in <100>. However, this loss in symmetry is neither a consequence of the vacancies nor their ordering. The resulting space g Bi2O3 are presented in Table 6 1. This superstructure contains three crystallographically distinct Bi sites, Bi(1)

PAGE 98

98 at 4a, Bi(2) at 4b and Bi(3) at 48h, and two crystallographically distinct O sites, O(1) and O(2), both at 96i. Ta ble 6 Bi2O3 determined from density functional theory under the GGA -PAW approximation. Space Group Fm 3 202 Unit Cell a o = 11.19761 Volume = 1404.0283 3 Atom Wyckoff Position x y z Occupancy Bi (1) 4a 0.0000 0.0000 0.0000 1.00 Bi (2) 4b 0.5000 0.0000 0.0000 1.00 Bi (3) 48h 0.2500 0.2229 0.0000 0.50 O (1) 96i 0.1279 0.0962 0.1460 0.25 O (1) 96i 0.4038 0.1460 0.1279 0.25 Taking the occupancy into account, the superstructure is comprised of 32 cations and 48 anions, i.e., 16 formula units per unit cell (Z=16). In the case of the cation, there are two different sites each of multiplicity four, Bi(1) and Bi(2), and one site of multiplicity 48 Bi(3). The occupancy of Bi(1) and Bi(2) sites is 1.0 and of Bi(3) is 0.5. There are thus 64 anion sites available in the structure, of which 48 are occupied by oxygen ions equally distributed on the two crystallographically different sites, O(1) and O(2); the remaining 16 site s are unoccupied (vacancies). The crystallographic sites of 1/8 of the superstructure are represented in Figure 6 4. As reference points, the conventional oxygen fluorite positions (8c in ideal fluorite) are shown in black in Figure 6 4; however, these poi nts do not correspond to ion positions in the superstructure. Each of these fluorite oxygen positions is surrounded by three O(1) or O(2) sites, only one of which is occupied at any given time. The oxygen atoms are displaced from the ideal 8c Wyckoff posit ion to 96i. To maintain the stoichiometry, only six oxygen anions are present in each 1x1x1 subunit; there are also two vacancies per subunit.

PAGE 99

99 Figure 6 4. Crystallographic sites of the space group Fm 3 (202). 1/8 of the 2x2x2 superstructure is shown here. Bi(1) and Bi(2) have an occupancy of 1.0. The Bi(3) sites have an occupancy of 0.5 with only one of the two of nearby available Bi(3) sites being filled at any given time. In the oxygen sublattice, for every oxygen fluorit e position (black) there are three possible oxygen positions, only one of which is filled either by O(1) or O(2). Color scheme: Bi(1) are shown in brown, Bi(2) in violet, Bi(3) in orange, O(1) in red, O(2) in blue and exact (8c) oxygen position in black. These displacements arise to accommodate the crystallochemical imbalance resulting from the unoccupied sites (vacancies). Further, since the vacancies are ordered, these displacements are also arranged in an ordered manner. 6.2.2 Bond Lengths Bond lengths are briefly discussed in the previous chapter. In this chapter, they are discussed in detail with the illustration of the implications on the structural distortion. The bond lengths predicted by DFT are presented and then compared with the estimates from the other, more approximate, method. The relationship between the different bond lengths and the vacancy ordering is also presented.

PAGE 100

100 6.2.2.1 Bond lengths from DFT calculations Bi2O3 combines vacancy ordering in <110> and <111> directions. To characterize the bonding in the structure, it is useful to first examine the structure from the perspective of the Bi ions. The Bi ions are in an octahedral environment with respect to the oxygen-vacancy cubic sub lattice, i.e., the Bi ions are at the centers of O ion cubes. In the combined vacancy-ordered structure, the Bi ions sit in two different environments, t he Bi(1) and Bi(2) ions have vacancy neighbors in the <111> directions (Figure 6 5), while the Bi(3) i ons have vacancy neighbors in the <110> directions (Figure 6 6). 6.2.2.2 Cation displacements in the <111> vacancy ordered oxygen sub -lattices For the <111> vacancy ordered oxygen sublattice, each Bi(1) and Bi(2) ion is surrounded by six O(1) and O(2) io ns respectively (Figure 6 5). As shown in Table 6 2, all the Bi(1) O(1) and Bi(2) O(2) bond lengths are identical. Figure 6 5. <111> vacancy ordered oxygen sublattice. Bi(1) present inside the O(1) sub lattice. All Bi -O bond lengths are 2.426 al l O -O distances are 3.138 and all O O O angles are 72.276 Similarly, Bi(2) is present inside O(2) oxygen sub -lattice and the sub lattice has the same geometry. Moreover, as discussed previously, this uniformity in bond lengths is reflected in the elect ronic degrees of freedom, with the charges between Bi(1) O(1) and Bi(2)-O(2)

PAGE 101

101 being equally distributed.89 This equivalence in the positions of the oxygen ions ensures that the Bi(1) and Bi(2) ions do not move from their crystallographic sites at the center of the cube. These are those non-displaced Bi ions previously discussed in Figure 6 1. Furthermore, all nearest neighbor O -O separations are 3.138 However, the equal Bi O bond lengths do not ensure a perfect -cubic (angles = 90 ) oxygen sub-lattice. Rath er, all the O O O angles in the sublattice are 72.276 a very significant deviation from cubic. 6.2.2.3 Cation displacements in the <110> vacancy ordered oxygen sub -lattice The oxygen ions in the <110> vacancy ordered sublattice are not all in crystall ographically equivalent positions. In this sub lattice, of the six oxygen ions surrounding each Bi(3), three are in O(1) sites and three in O(2) sites. As a result, the <110> oxygen sub lattice is heavily distorted and has three different Bi O bond lengths (see Table 6 2). One of the twelve possible <110> vacancy ordered oxygen sub-lattices is shown in Figure 6 6(a). In this sub lattice, there are three O(1) ions, marked as O(1) 1, O(1) 3 and O(1) 6 and three O(2) ions, marked as O(2) 2, O(2) 4 and O(2) 5. There are three different Bi -O bond lengths that occur in pairs; each pair contains two crystallographically different oxygen ions. In Figure 6 6(a), the Bi(3) O(1) 6 and Bi(3) O(2) 5 bond lengths are 2.303 Similarly, the Bi(3) O(1) 1 and Bi(3) O(2) 2 bond lengths are 2.399 while the Bi(3) O(1) 3 and Bi(3) O(2) 4 bond lengths are 2.559 A consequence of the non-equivalence of the oxygen ions and <110> vacancy ordering is that the central Bi(3) ion is displaced from its ideal positions (shown by dott ed arrow in Figure 6 6(a)) by 0.3036 in one of the <100> directions away from the face of the distorted cube containing the vacancies. In contrast to the simpler geometry of the <111>

PAGE 102

102 oxygen sub lattice, the <110> sub-lattice is heavily distorted with th ese differing bond lengths also resulting in dissimilar O O O angles (Figure 6 6(b)). Figure 6 6 <110> vacancyordered oxygen sub-lattice. A) Bi(3) present inside the oxygen sub lattice formed by three O(1), three O(2) and two vacancies. There are three different pairs of Bi O bond lengths: Bi(3) O(1)6 and Bi(3) O(2)5 = 2.303 Bi(3) O(1)1 and Bi(3) O(2)2 = 2.399, Bi(3) O(1)3 and Bi(3)O(2)4 = 2.559 The dotted arrow shows Bi(3) displacement in <100> direction by 0.3036 away from the plane containing the vacancies and solid arrows on O(2)5 and O(1)6 show the oxygen displacement in a direction close to <110>. B) Schematic showing the geometrical values of the distorted <110> oxygen sub -lattice. Only O O distances and O O O angles are shown; those involving a vacancy are not shown. 6.2.2.4 Anion d isplacement The DFT calculations show a preferential -directional displacement of all of the oxygen ions. The structural analysis shows that every oxygen ion displaces in the direction of the two fir st nearest vacancies ordered in <110>, as shown by the solid arrows in Figure 6 6(a). The se two vacancies provide considerable open space in the structure, into which the ions can relax. The result is that both oxygen ions displace from the fluorite site b y a magnitude of x = 0.235, y = 0.322 and z = 0.033, resulting in a total displacement of 0.4 in a direction close to, but not exactly parallel to, <110>. A B

PAGE 103

103 6.2.2.5 Bond lengths from Shannon radii There are several standard ways to estimate the b ond len gths in ionic materials, including analysis of the Shannon radii.69 In this sub -section section the value of Shannon radii method for analyzing Bi2O3 is assessed The bond lengths calculated from the radii given by Shanno n69 are based on the coordination environment. In ideal structure s without defects, the bond lengths predicted by Shannon radii are quite accurate.69 However, in the presence of defect s or vacant lattice site s as in the Bi2O3 structure the bond lengths can either deviate from the expected values or sometimes can provide only incomplete information, as is found to be the case in the following analysis E very cation i n an ideal fluorite structure is coordinated with eigh t anions. However in Bi2O3, each Bi3+ ion is coordinated with only six oxygen ions, the other two anion sites being vacant. A ccording to the Shannon radii,69 rBi 3+(VI) = 1.03 rBi 3+(VIII) = 1.17 and rO 2-(IV) = 1.38 This yields estimates of the Bi(VI) O of 2.41 This estimated bond length is actually very close to the Bi O bond length of 2.424 in the non -distorted <111> ordered oxygen sublattice. Interestingly, although six coordinated, the DFT calculations yie ld Bi O bond lengths in the <110> sub lattice as 2.559 which is very close to the Bi(VIII) O bond length of 2.55 Hence, in this case the Shannon radii predictions are in close agreement with some Bi O bond lengths given by DFT, but not others. 6.2.2 .6 Sub-unit Cell Distortion Thus far, it is characterized that ionic displacements of both Bi and O ions from their ideal fluorite positions Bi2O3. As a consequence,

PAGE 104

104 the structure is distorted locally. In this section, th ese local distortions from the perspective of the FCC network of Bi ions are characterized. Figure 6 7. Distortion in the superstructure. A) Two interpenetrating FCC networks in 2x2x2 superstructure. B) A non-distorted cubic FCC network formed by four Bi(1), four Bi(2) and six Bi(3) cations. C) A distorted FCC network formed by twelve Bi(3); the other two cations can be either both of same Wyckoff number or different (Bi(1), Bi(2)). The corners of the non -distorted FCC network are formed by Bi(1 ) and Bi(2) whereas those of distorted FCC network are formed by Bi(3) cations. F or clarity o xygen ions are not shown. Color scheme Bi(1) in brown, Bi(2) in violet and Bi(3) in orange. As Figure 6 7(a) shows, the Bi network can be considered either as an FCC network with Bi(1) and Bi(2) ions in the eight corners and Bi(3) ions of the six phases, or as a network with Bi(3) ions at the eight corners and on four of the six faces, the remaining two faces being occupied by a Bi(1) and a Bi(2) ion. These two distinct, interpenetrating FCC networks are shown in Figure 6 7(a). The difference between the two networks can be characterized by the shape of the boxes defined by the corner Bi ions. A B C

PAGE 105

105 Table 6 2. Bi O bond lengths in the <110> and <111> vacancy ordered s ub-lattice. The two <111> vacancy ordered sublattices has Bi(1) and Bi(2) ions surrounded by six O(1) and six O(2) ions respectively. The <110> vacancy ordered sublattice has three O(1) and three O(2) ions surrounding each Bi(3) ion. <111> () Bi(1) O O(1) 2.426 Bi(2) O O(2) 2.426 <110> Bi(3) O O(1) 5, O(2) 6 2.303 O(1) 1, O(2) 2 2.399 O(1) 3, O(2) 4 2.599 Recalling that th e Bi(1) and Bi(2) ions sit at the ideal fluorite sites, the box defined by these ions, Figure 6 7(b), is a perfect cube with the Bi Bi second nearest neighbor distance being half the lattice parameter, 5.5988 By contrast, the FCC network defined by the Bi(3) ions in the corners, which are displaced by 0.3036 along <100> as shown in Figure 6 7(c) is distorted. In Figure 6 7(c), the Bi -Bi distance between ions 2 and 20 and between ions 4 and 21 is not affected by the displacement and remains 5.5988 as is in the undistorted Bi(1) Bi(2) FCC sublattice. However, the anti -parallel displacement of the columns of ions has the consequence that the distance between ions 2 and 4, and between ions 20 and 21 ions increases to 5.6316 (Similar anti -parallel di splacements occur between ions 2 and 15 and between ions 20 and 18.) These displacements distort 96.18 Each 2x2x2 superstructure contains an equal number of these di storted FCC networks oriented along crystallographic axis. Thus, despite these distortions, due to the anti -parallel displacements in all three crystallographic axis, the overall system is cubic. In the distorted FCC network, 12 of the 14 cations are Bi(3) ions with the remaining two being either Bi(1) or Bi(2). Depending on the choice of the Bi(3) FCC network, both of the nondisplaced ions could be at their Wyckoff positions or one could

PAGE 106

106 be in a Bi(1) site and the other in Bi(2) site. In the non -distorte d FCC network (Figure 6 7 (b)), 8 out of the 14 Bi FCC ions are non-displaced, i.e., four each from Bi(1) and Bi(2); the remaining six are Bi(3). Battle et al.90 also postulated a similar kind of FCC network to the non-distorted network that has been prese nted here. In the Y2O3-Bi2O3 system, they expected Bi3+ to occupy all Bi(3) positions (as also seen here), whereas, the rest of the positions at the corners of the unit cell (see Figure 67(a)) could be occupied either of Y3+ and Bi3+. In their analysis, t hey did not, however, distinguish between the non-distorted and the distorted networks. 6.2.2.7 Comparison with Bixbyite Structure The bixbyite structure displays Ia 3 space group symmetry with a cell formula M32O48 consisting of 32 cations, 48 anions and 16 vacant anion sites. The cations occupy 8a and 24d sites at 0, 0, 0 and x 0, 0.25 respectively and the anions occupy 48e sites at x, y, z It is interesting to notice that while the vacancy ordered structure of -Bi2O3 very closely resembles that of the bixbyite structure, it is actually distinct. On the one hand, the vacancies and the cations have similar crystallographic positions and the vacancies have same site occupancy in both the structures, thus, result ing in the same ordering of vacancies. Using coordination defect theory, Bevan and Martin91 have recently shown 3D vacancy networks in <110> directions in bixbyite Mn2O3. These vacancy networks are connected by <111> vacancy linkages. This is identical to the <110> <111> vacancy in Bi2O3.75 89 The cations show similar crystallographic characteristics in both the structures. In particular, in each structure there

PAGE 107

107 are 8 non-displac ed and 24 displaced cations. Further, the vacancy environment around two different cations remains the same (Figures 6 5 and 6 6). On the other hand, the anion positions in the two structures are different, arising from different displacements relative to the cation positions. To elucidate these differences, a DFT simulation starting with Bi2O3 arranged in the bixbyite structure is performed. During the relaxation, the O ions move considerably; the final structure is that originally obtained structure for Bi2O3 with space group symmetry Fm 3 It is quite possible that the differences in the anion sublattice, responsible for the difference between the two structures, arises from the lone pair electrons of Bi, that are absent in other cations crystallizing in bixbyite structures. 6.3 Relation to Experiment The issue of the overall nature of the defect ordering is addressed first. MD simulations show that the high polarizability of Bi ions leads to continuous diffusion in the system. As a result the combined <110> <111> ordered defect structure does not become locked in. The simulations further show that the presence of non-polarizable or weakly polarizable cations does result in the ordered defect structure being produced, and that it h as a structure that is consistent with that found in the experimental analysis of doped Bi2O3 by Wachsman et al.75 It is this structure that is analyzed here and it can be considered as the prototype for the structure present in the doped system. It is thus appropriate to place the crystallographic analysis in the context of the experimental results for doped systems. Anion Displacement: Positional Ordering from DFT

PAGE 108

108 Bi2O3 has long been a subject of discussion with several contradictory models being presented.68 74, 90, 92 A detailed discussion of these models can be found in Boyapati et al.73 Here the issue of doped Bi2O3 is addressed, i.e., the anion displacement within the defect ordered structure. Based on neutron-d iffraction data, Willis et al.74 presented an oxygen sub-lattice model (anion cubic sub lattice) in which all the anions were suggested to be displaced in <111> directions, i.e., a displacement from regular t etrahedral site, 8c (, ), to an interstitial site, 32f (0.3, 0.3, 0.3). Battle et al.84 suggested a similar displacement in pure Bi2O3. In Y2O3-Bi2O3 92 however, they found th at a higher concentration of yttria decreased the concentration of <111> displaced anions, with some evidence of <110> anion displacement. By contrast, the experiments by Boyapati et al. have, however, found that in the displacement occurs in <111> direct ions; they did not observe any signature of <110> displacement. An incr eased concentration of less polarizable dopant s (Y3+ compared to Dy3+) increases the number of <111> displaced anions. As Boyapati et al.73 point out, t he difference between their results and those of Battle et al. may be due to the fact that their samples were not aged while those of Battle et al. were aged. DFT calculations also reveal anion displacement that are close to but not exactly along <110>, an d the observations in this thesis are thus are to a large extent consistent with the results of Battle et al .90, 92 The actual defect structure of -Bi2O3, whic h is only stable at high temperature, has not been fully characterized. While it is possible that this structure, evolving through diffusion, is not related to the ordered structure inherent to the doped systems, it also seems reasonable that at any single instant in time for any small region of the crystal, the

PAGE 109

109 structure would strongly resemble one the crystallographically equivalent variants of the defect -ordered structure. Further structural characterization and simulation are required to resolve this is sue. 6.4 Discussion MD simulations and DFT calculations are consistent with experiments in that Bi2O3 has a 2x2x2 fluorite related superstructure with a combined vacancy ordering in <110> and <111> directions. It is shown that vacancies lead to significa nt structural changes, particularly, the ionic displacements. Battle et al.84 reported the comparison between the displacement of the cations in pure Bi2O3 and Y2O3 doped -Bi2O3. In the pure Bi2O3 the cations were found to be present at the exact crystallographic positions.84 By contrast, in the doped -Bi2O3 cations were displaced in directions. Similar cation displacements in directions are a lso found here. Further analysis shows that the displacement of the cations is correlated with the ordering of the vacancies around them. All cations, which are displaced have an ordering of the vacancies in <110> around them. By contrast, all non-displac ed cations have vacancies ordered in <111> directions. All the cations that displace in either of the directions, i.e., [100], [010] or [001] move by 0.3036 from their crystallographic positions. Given the partial occupancies of some anion and cation s ites, and the vacancy ordering previously described, Table 6 1 fully describes the atomic and vacancy arrangement. The previous analysis72 73 of neutron diffraction data within the fluorite space group considered three possible positions for the oxygen ions 8c, 32f and 48i sites within the fluorite space group. The refinement identified the 32f sites as being most

PAGE 110

110 consistent with the data. The analysis of the positions derived from the simul ation is based on the 2x2x2 fluorite superstructure and identifies the oxygen ion positions as occupying the 96i sites. The results of these experiments showed that not only was there occupancy ordering where vacancies aligned along <111> and <110>, but that there was also positional ordering where the oxygen ions were displaced in <111> from the normal tetrahedral 8c site to interstitial 32f sites.72 73 In the current simulation study th e oxygen sites are indexed with respect to the 2x2x2 Fm 3 space group, which has a larger number of equivalent sites. Thus, the Fm 3 m 32f site is analogous to the 96i in the Fm 3 structure. Th e physical displacement from the 8c sites to the 32f sites in the Fm 3 m space group implies an equal possibility of displacements along the faces of the coordination tetrahedron or the <111> directions. However, in the Fm 3 space group, due to the doubling of the unit cell to 2x2x2, displacements from the ideal tetrahedral position (which is 32f in this space group rather than 8c) to the 96i sites represent in -plane displacements within the {111} planes, i.e., the oxygen ions are displaced towards the edges of the coordination tetrahedron. Due to lower symmetry of the superstructure, for a given oxygen ion, displacement is not possible towards all six edges of the tetrahedron but only towards three of them.93 As su ch, by taking into consideration the doubling of the unit cell, the superstructure presented here represents and captures the correlated anion and cation displacements in even more detail. 6.5 Conclusion and Outlook In Chapters 4,5 and 6 a detailed charact erization of the structure of Bi2O3 from molecular -dynamics and electronic level simulations is presented that enabled a complete crystallographic analysis. MD simulations have shown that atomic polarizability is very

PAGE 111

1 11 important in the oxygen diffusion. Low polarizability leads to a c ombined ordering of the vacancies in <110> and <111> directions; high polarizability leads to a disordered lattice and sustained diffusion. The equilibrated vacancy ordered structure determined from MD simulations is consistent with the order of energies d etermined from DFT calculations. Crystallographic analysis shows that a 2x2x2 superstructure is required to fully describe the vacancy ordered structure of Bi2O3, the crystallographic parameters of which are given. From the electronic structure viewpoint, it is shown that the charge distribution plays a very important role. Due to two different kinds of vacancy environments around Bi atoms, the Bi atoms bo nd differently with O atoms. It is found that Bi forms covalent character bond when vacancies are symmetrically ordered in <111> direction. In contrast it maintains high ionic character bind when vacancies are ordered in <110> direction. IT is speculated t hat because covalent bond is directional and difficult to break, it is the limiting factor in oxygen diffusion. The cessation of the oxygen diffusion is a consequence of the bond formation. Finally, although this work has laid the foundations in underst anding the structure of pure Bi2O3, it is important to understand the mechanisms of diffusion and the difference in the doped systems like YSZ and doped CeO2 to fully characterize the fluorite -based materials.

PAGE 112

112 CHAPTER 7 RADIATION DAMAGE IN MGO 7.1 Nucle ar Energy 7.1.1 Nuclear Power A large amount of energy can be released by the disintegration of an atom. In the year 1911, Ernest Rutherford, brought this to attention from the decay of radium and stated:94 This evolution of heat is enormous, compared with that emitted in any known chemical reaction The atoms of matter must consequently be regarded as containing enormous stores of energy which are only released by the disintegration of the atom In the last hundred years, science has come a long way by de veloping a good understanding of the nuclear fission reactions. It is recognized that a wise use of nuclear energy has immense potential to be a part of the solution to the energy crisis; the most significant nuclear energy contribution being in producing electricity. Currently there are 443 nuclear power plants operating throughout the world that account for 15% of the global electric power generation. With nuclear reactors becoming safer95 and with increasing technological understanding, this sustainable e nergy -production source can be relied upon more. Moreover, with a degrading environment and depleting fossil fuel resources, the world is hard -pressed to advance in this technology. 7.1.2 Challenges for Materials Science While the understanding on nuclea r reactions has been put to applications in nuclear reactors for the last 50 years, we still face enormous challenges from a materials science viewpoint. Two most fundamental challenges are, (a) improving radiation tolerance of fissionable materials to imp rove fuel life inside the reactor, (b) safer disposal

PAGE 113

113 or recycling of the nuclear waste. The first problem is relatively less worrisome than second, for which, we still have no scientific longterm solution. The nuclear waste consists of radioactive materials that emit different kinds of harmful radiations. Materials such as 233U and 239Pu have long half -life period of 160,000 and 24,100 years.96 In addition, minor actinides radioactive such as 237Np, 241Am+243Am and 244Cm are produced in the nuclear fissi on reaction and a pose direct threat to environment, not only by radiation but also by mixing with the underground water. Consequently, in pursuit of preventing the environment degradation from carbon-based fuels, all of a sudden one has to worry about the pollution due to their alternatives.97 Two strategies are broadly envisioned for dealing with nuclear waste: recycling of spent nuclear fuel and direct disposal into a repository. The first strategy involves reprocessing of the nuclear fuels to reclaim U and Pu for the fabrication of mixed oxide fuel (MOX). The actinides can be incorporated into an inert matrix fuel (IMF). After a once through burn up, the MOX and IMF would be sent to a geological repository. The second strategy involves direct disposal o f the spent fuel in a geological repository. The U, Pu and actinides would be stored in materials such as zircon98 and fluorite -based materials99 that are resistant to damage due to radiation emitted by these radioactive materials. These can then be deposited in the geological repositories. Both these strategies are still under research and much needs to be understood, particularly from the radioactivity viewpoint, before they can be widely applicable. 7.2 Radioactivity From Coulombs Law, the positively cha rged protons in the nucleus should repel each other and the nucleus should rapidly fall apart. However, at distances smaller than 1015 m, the protons and neutrons are bound together by a powerful attractive force called

PAGE 114

114 the strong nuclear force, which kee ps the nucleus together.100 This force is independent of the charge on the species. As the number of protons increases inside the nucleus, to maintain its stability, the number of neutrons also increases. For low atomic number elements, the number of neutrons is approximately equal to the number of protons. However, for higher atomic number elements even more neutrons are required to stabilize the nucleus (e.g., 27 59Co 28 60Ni ). However, eventually for elements with very high atomic number, even the increased number of neutrons is no longer able to provide such stability. Hence, to gain stability, the nucleus can disintegrate by radiation. This process of disintegration is known as radioactive decay. 92 235U 94 239Pu are common elements that disintegrate through radioactivity. There are four kinds of radiation decay by which a nucleus can disintegrate: alpha ( ) decay, beta ( ) decay, gamma ( ) decay and neutron decay. 7.2.1 Alpha ( ) Decay When a nucleus disintegrates by emitting an alpha particle ( 2 4He ), the process is known as alpha decay. Equation 71 describes the alpha decay of 94 239Pu 94 239Pu 92 235U 2 4He (7 1) In the -decay, the parent nu cleus ( 94 239Pu ) is different from the daughter nucleus ( 92 235U ). Hence, one chemical element is converted into another. The binding energy librated due to the reaction is converted into kinetic energy of the alpha particle and the daughter nucleus. Alpha particles usually have kinetic energy of ~5 MeV and mean free path length of ~32 m. Because of their small path length, they can be stopped by a sheet of paper and hence are not harmful to human tissue, unless inhaled.

PAGE 115

115 7. 2.2 Beta ( ) Decay Beta decay occurs if a nucleus is unstable due to too many protons or neutrons.100 The excess neutrons or protons transform into other nucleon, changing the neutron to proton ration and stabilizing the nuc leus. There are three types of beta decays such as beta -minus decay, beta -plus decay and electron capture. In the beta -minus decay, a neutron transforms into proton, and emits an electron called as beta minus particle ( -) and an electron antineutrino ( e ). A beta -minus decay can be represented as: 0 1n 1 1p 1 0 e (7 2) In the beta plus decay, a proton changes into neutron by emitting a positron ( +) and an electron neutrino ( e ). The positron i s similar to electron but has an opposite charge to electron. The reaction can be given as: 1 1p0 1n 1 0e (7 3) The third type of beta decay is electron capture. In this decay, a proton captures an orbiting electron and transforms into a neutron while emitting a neutrino. The electron capture decay can be represented as: 1 1p 1 0e 0 1n e (7 4) In all types of beta decay, the parent nucleus has a different atomic number than the daughter nucleus. Equation 7 5 gives a c ommon process of all three types of beta decay occurring in 52 126I and 55 126Cs Both these elements are short -lived radioactive nuclides. 52 126Te 53 126I 54 126XeEC 55 126Cs (7 5)

PAGE 116

116 7.2.3 Gamma ( ) Decay In -decay, a nucleus changes from a higher energy state to a lower energy state releasing energy in the form of packets of electromagnetic radiation known as gamma rays. These rays are the most intense and penetrative among all. In gamma decay, the number of protons and neutrons remain same. 7.2.4 Neutron Capture In this radioactive process, a nucleus captures a neutron and forms an excited nucleus.101 The excited nucleus then emits gamma rays and returns to the normal state. A neutron capture reaction can be represented by equation 7 6. 0 1n 92 238U 92 239U 92 239U92 239U (7 6) A target nucleus, 92 238U captures a neutron and forms 92 238U *, which is an excited state of 92 238U It is an important process in the product ion of 94 239Pu in reactors. 7.3 Nuclear Fission In a nuclear fission reaction, a nucleus divides into two fragments (of smaller nuclei) and three neutrons. A typical fission reaction of 235U can be given as: 0 1n 92 235U 56 144Ba 36 89Kr 3 n0 1 (7 7) For every neutron involved in the nuclear reaction, three new neutrons are produced that further cause fission reactions, in process, setting up a chain reaction. Hence, nuclear fission reactions can be self -sustaining reactions. The fission products have very high kinetic energy; their slowing down is the main mechanism of heat production in the reactor. The energy obtained from a nuclear reaction can be roughly estimated by using the equation E=mc2. For equation 7 7, the energy produced can be estimated as:

PAGE 117

117 E [ m (235U) m (144Ba ) m (89Kr ) 2 mn] c2173 MeV (7 8) Considering other reactions such as kinetic energy of beta rays, energy of gamma rays, energy of neutrinos and the energy associated with the delayed neutrons, the total energy of a fission reaction of 235U is approximately 195 MeV per fission event.101 This energy is much higher than energy from normal chemical reaction which involve energy transfer of ~5 10 eV. 7.4 Radiation Damage In a nuclear reactor, with all of the ab ove mentioned processes occurring simultaneously, during the fuels lifetime, large amount of fission products (10%) are accumulated and the fuel undergoes 2000 3000 displacements -per atom (dpa),102 thus causing gross changes in the microstructure of the n uclear materials such as formation of voids and cracks. In addition, the crystalline nature of the materials is also damaged due to the creation of vacancies and interstitials, collectively known as Frenkel Pairs (FPs), thus limiting the lifetime of the fu el. Similar processes also occur in nuclear waste. The by -products (or nuclear waste) undergo radiation through -decay of fission products (137Cs and 90Sr) and -decay of actinides (U, Np, Pu, Am and Cm). While -particles produce very low energy recoil nuclei and rays, -decay produces high -energy particles (4.5 to 5.5 MeV) and energetic recoil nuclei (70 to 100 keV) along with rays. Due to their high energy, -particles undergo elastic collisions and the cascade effect displaces several thousand atoms, thereby creating FPs. Furthermore, due to the longer half lives of actinides, -decay is generally dominant at longer times.103 Understanding the physical and chemical processes due to the consequence of radiation damage, both inside the nuclear reacto r and in the repository conditions, is necessary to improve the

PAGE 118

118 performance of the fuel in the reactor and to produce longlasting storage materials that can sequester radioactive materials/rays from the environment for thousands of years. The radiation damage phenomena take place at very small time scales (< 1012 s), which makes it difficult to understand them from experiments. Molecular dynamics (MD) simulation is an easy and informative method of calculating the atomistic damage, both chemical and phy sical. In this thesis, MD has been extensively used to elucidate the atomistic phenomenon occurring due to radiation damage. 7.4.1 Standard Radiation Damage Methodology Using MD Simulation Atomic level simulation has long provided important insights into the fundamental processes associated with radiation damage .104105, 106, 107, 108 In particular, simulations ha ve shown themselves to be capable of capturing two important phases of radiation damage. First, the ballistic phase involves the formation of the collision c ascade arising from the high -energy primary knock on atom (PKA) (or the recoil of alpha particle causing atomic displacements). This phase typically lasts only a few picoseconds, its net effect being the generation of FPs and small point -defect clusters.109 Second, the kinetic phase captures the diffusion-controlled dynamical evolution of these FPs and clusters, their annihilation, and the formation of larger longlived defect clusters. The evolution of FPs during these two phases is represented in Figure 7 1. During the ballistic phase, a sudden increase in the defects is represented by a rising curve (Figure 71(a)). The number of defects reaches a maximum (shown by peak) during this phase. At the same time, the high kinetic energy (KE) of the PKA is converted into potential energy (PE) of the defects (Figure 7 1(b)). In Figure 7 1(b), the peak time in the number of defects (also same in Figure 7 1(a)) matches with the peak -time in the PE, thus illustrating the conversion of PKA KE to maximum PE of defec ts.

PAGE 119

119 Figure 7 1. Two phases of radiation damage. A) Evolution of defects during the ballistic and kinetic phase. B) Distribution of kinetic energy (KE) and potential energy (PE) during the evolution of defects. Number of defects in (A) is a magni fied view of defects curve in (B). Figure 7 2. Typical Frenkel pair (FP) defects during radiation damage using PKA. A) Defects formation during the ballistic phase. B) Left -over defects during the kinetic phase. This phase lasts only for few pi coseconds after which, the FPs recombine and annihilate. The fall in the number of defects (and PE) represents the recombination of the FPs. During this time, a gradual transition from the ballistic phase to the kinetic phase takes A B A B

PAGE 120

120 place. During the kineti c phase, most of the energy has been dissipated throughout the system (in the simulation, a thermostat is used to dissipate the energy) and the left-over defects diffuse to either recombine or form cluster. Schematics of the defects during the ballistic ph ase and kinetic phase are shown in Figure 7 2. Large number of defects in Figure 72(a) represents the peak during ballistic phase. Figure 7 2(b) represents the leftover defects during the kinetic phase. Most of the defects created during the ballistic ph ase recombine and annihilate. Thus, it is this kinetic phase that determines the long -time, experimentally accessible behavior of the material: while the recombination events promote radiation tolerance, the cluster -formation process produces the irreversi ble damage that degrades materials performance and ultimately limits lifetime. Previous simulations of collision cascades followed by rather long temperature accelerated dynamics (TAD) simulations107 110 have shown that the defect evolution during the kinetic phase is largely independent of the detailed nature of the initial damage created during the ballistic phase, suggesting that on a longer time scale, the role of the ballistic phase in the cascade simulati on is mostly to introduce non -equilibrium point defects into the system. 7.4.2 K inetically -evolving Irradiation -induced Defects Method Standard collision cascade radiation-damage simulations involving a PKA introduce only rather a small number of such poi nt defects into the system. Thus a key question is whether these simulations are able to produce the wide range of defect environment present in experiments. Moreover, because it has been shown for some ionic systems, such as UO2, that the initial propagat ion direction of the PKA has little effect on the ultimate defect evolution,111 it appears that, at least in some cases, a rather complete

PAGE 121

121 understanding of defect recombination, annihilation at sinks and clustering can be obtained by focusing on the kinetic phase alone. While standard radiation -damage MD simulations using the PKA approach capture the formation of the defects they are not ideal for capturing the longer time evolution of the systems. Most importantly, typically they are not long enough in dur ation to allow the formation of complex defect structures. Here an approach is described that, by circumventing the ballistic phase entirely, allows a much wider range of different defect environments to be explored. In this method, at the beginning of eac h simulation, a specified number of point defects (FPs) are created in a high -temperature equilibrated system To avoid uninteresting and spontaneous interstitial -vacancy recombination events, the vacancies and interstitials are separated by a distance gre ater than the recombination radius. The system is then let to equilibrate at a high temperature and the diffusion -controlled kinetic evolution of the system is followed. It can be argued that this method results in unrealistically large concentrations of defects. While this may be true when averaged over time and space, statistically it can be expected that such high defect concentrations will on occasion be generated over some spatial region. Moreover, the resulting defect structures, unlike rapidly reco mbining point defects, can have a significant role in determining the true radiation performance of the material. This approach does not replace cascade simulations aimed at elucidating the effects of experimental conditions, such as the mass, energy and dose of the incident species on initial defect production. Rather, it is complementary in that it aims to capture more comprehensively the wide range of defect processes and kinetically evolving defect

PAGE 122

122 environments that are rather insensitive to the detail s of specific cascade events. This approach allows to observe defect evolution under different defect environments and reveals the governing mechanism in the radiation tolerance. In application to MgO as a model material, it is demonstrated that this appr oach replicates the defect clusters seen in the full cascade simulation followed by the TAD simulation.107 110 More importantly, it exposes a novel, previously unidentified intricate cluster -formation mechanism that involves not only the original, radiation induced FPs but also formation of new structural FPs as a lattice response that stabilizes the larger clusters. This approach thus offers a simple, computat ionally straightforward and physically transparent way of elucidating the kinetic interplay among defect recombination, annihilation at defect sinks and clustering in the development of radiation damage. 7.4.3 Materials Under Study In this study, MgO and UO2 are chosen to understand the radiation damage in ionic materials. MgO is chosen as a model material to prove the fidelity of the new method. It is motivated by the availability of previous simulations by Uberuaga et al.107 110 of relatively low PKA -energy (0.4 5.0 keV) cascade simulations that revealed that most of the FPs formed in the ballistic phase recombine during the kinetic phase. T hose FPs that did not recombine were seen to evolve kinetically to form highly -stable clusters, their size ranging from di -interstitials to a maximum of seven ions Particularly noteworthy is the fact that these clusters were not formed during the ballistic phase but as a co nsequence of interstitial diffusion during the kinetic phase in which the defects diffused long distances and were not confined to the region affected by the initial ballistic cascade.

PAGE 123

123 From the technological viewpoint, MgO is a widely used as engineering material and is generally very well understood both experimentally and computationally. Furthermore, robust MgO empirical potentials are available that capture point defects phenomenon very well. UO2, on the other hand, is a key component for fission react ors in nuclear -fuel rods. Understanding radiation damage in UO2 is important to enhance the fuel performance. Also, since it has a fluorite crystal structure, it also serves as a model material to understand the radiation damage in the envisioned fluorite related host materials99 for the nuclear waste. This chapter illustrates radiation damage in MgO where, not only the new radiation damage methodology is demonstrated but also the key cluster -formation mechanisms in MgO are r evealed, which have not been known before. In Chapter 8, radiation damage and clustering processes in UO2 are revealed. 7.5 Interatomic Potential and MD Simulation Our MgO single -crystal simulation cell consists of 20x20x20 cubic rocksalt unit cells, con taining 64,000 ions. Periodic boundary conditions are applied in all three directions. The time step of 0.5 fs ensures good energy conservation in tests in the NPE ensemble. The short -range (mainly repulsive) interionic interactions are described by Lewis and Catlows112 rigid ion Buckingham -type potential. The longrange electrostatic interactions are evaluated via direct -summation method38 with spherical truncation at the cut off radius of 8.139 The zero-temperature lattice parameter thus obtained is a0=4.1986 All defect simulations are carried out in the NPT ensemble at 1000 K (the melting temperature of MgO for the potential used here is ~ 3250 K).113 At 1000 K, both Mg2+ and O2interstitials are highly mobile, whereas the corresponding vacancies are completely immobile on the MD time scale.

PAGE 124

124 7.6 Results At the beginning of each simulation, a specified number of FPs is randomly distributed throughout the system; to ensure their stability, i.e., to minimize uninteresting a nd spontaneous recombinations, the interstitials and vacancies are initially separated by a distance of at least 5a0. Because they cannot mutually annihilate, there is no constraint on the vacancy -vacancy or interstitial interstitial distances; thus some a re actually quite close to each other. In order to elucidate the nature and mechanism of defect clustering during the kinetic phase, two distinct types of starting configurations are considered. In the first, FPs are only inserted into either the Mg or the O sub lattice, but not both; in the second, equal numbers of FPs are inserted into both sub -lattices. 7.6.1 FPs Only on One Sublattice Initially, 100 FPs of Mg are created and the system is allowed to thermally equilibrate at 1000 K. After 2450 ps, only 18 Mg FPs are present The schematic evolution of the defects on Mg sublattice is shown in Figure 7 3. The initially created 100 Mg FPs e ventually recombine by interstitial -vacancy recombination. No Mg interstitial clustering is observed; as is shown be low, this is due to the absence of oppositely-charged O FPs. Similarly, 100 FPs were created on the O sublattice. As in the case of the Mg FPs, eventually complete annihilation takes place with no clustering occurring on O sub-lattice due to the absence of the Mg FPs (Figure 7 4). In a conventional collision-cascade simulation it is, of course, not possible to selectively create defects on only one of the sublattices; this demonstrates a unique aspect of this approach in understanding the defect evoluti on. These simulations provide a baseline against which to compare the results of simulations of FPs on both sub-lattices.

PAGE 125

125 Figure 7 3. Evolution of FPs on Mg sublattice. Snapshots taken at A) t=0 and B) t=2450 ps showing the evolution at T = 1000 K of 100 FPs initially introduced on Mg sublattice Of the initially 100 Mg FPs in (A), only 18 Mg FPs are left behind in (B). Dark green and orange spheres represent Mg interstitials and vacancies; eventually all Mg FPs recombine and annihilate. F igure 7 3. Evolution of FPs on O sublattice. Snapshots taken at A) t=0 and B) t=2450 ps showing the evolution at T = 1000 K of 100 FPs initially introduced on O sublattice Of the initially 100 Mg FPs in (A), only 26 Mg FPs are left behind in (B). Red and blue spheres represent O interstitials and vacancies. Eventually all O FPs recombine and annihilate. 7.6.2 FPs on Both Sublattices Next, equal numbers of FPs are simultaneously inserted into both the Mg and O sub -lattices. This scenario corresponds to elec tron radiation commonly used to produce single vacancies and interstitials.114 It also resembles the situation at the end of the ballistic phase of a collision -cascade radiation damage simulation. Figure 7 -5 shows two A B A B

PAGE 126

126 snapshots taken again at (a) t=0 and (b) t=2450 ps for a system containing initially 100 FPs on each of the two sub-lattices. Compared to the case of Mg or O FPs only, this annealing simulation reveals much less interstitial -vacancy recombination but, instead, significant interstitial aggregatio n into clusters. This aggregation starts quickly by the formation of highly mobile di and tri interstitials. Identical di and tri interstitial clusters were observed in the lower temperature, combined collision -cascade and TAD simulations by Uberuaga et al.107, 110 Figure 7 5. Snapshots taken at A) t=0 and B) t=2450 ps showing the evolution at T = 1000 K of 100 FPs initially introduced into each of the two sub -lattices. Of the initially 100 Mg and 100 O FPs in (A), 72 Mg and 76 O FPs are still present in (B). While nearly all of the Mg2+O2 di interstitials combine to form larger clusters, only a few of the tri interstitials (of both kinds: Mg2+O2 -M g2+ and O2 --Mg2+O2-) aggregate to form larger clusters; the others maintain their individual tri interstitial identity. Both the tri -interstitial and the larger -sized interstitial clusters are found to be very stable in that they survive to the end of the simulation at 2450 ps. During this annealing simulation, the majority of FP annihilation events by vacancy -interstitial recombination occur very quickly (within the first 50 ps), during A B

PAGE 127

127 which period each of the defect types is fairly evenly distributed in space. Once interstitials aggregate into clusters, the likelihood of further vacancyinterstitial recombination quickly decreases in favor of cluster growth by interstitial capture. Any remaining single interstitials therefore quickly disappear from the system; in fact, after about 2000 ps, virtually no individual interstitials are left. This is quantified in Figure 7 6, which compares the number of Mg2+ interstitials present in a system with FPs only on the Mg sublattice with a system with FPs on bot h sublattices Figure 7 6. Total number of Mg FPs present in the simulation cell vs. time when the system is initialized with Mg defects only (triangles), and both Mg and O defects (squares). This comparison demonstrates that, when the FPs are p resent on only one sub lattice, all the vacancies and interstitials eventually recombine due to the absence of clustering; hence, after 2450 ps only ~20% of the initial FPs remain By contrast, in the system containing FPs on both sublattices, about 80% of the FPs become immobilized by precipitation into highly stable vacancy -interstitial clusters; in fact, after only 1000 ps the FP count levels off (see Figure 7 6).

PAGE 128

128 The stability of these clusters is most likely strongly enhanced by being stoichiometric, and hence charge neutral. Similar simulations with very low concentration of point defects simultaneously on both the sub-lattices also result in di and tri interstitial clusters that maintain their identity for the entire 2500 ps of the simulation. 7.7 Cluster -Formation Mechanism A close up snapshot of one of the larger interstitial clusters, circled in Figure 7 5(b), is displayed in Figure 7 7 Investigation of the dynamic growth mechanism of this representative cluster reveals two distinct, coordinate d mechanisms that control its formation and growth. Figure 7 7. Close up of the all interstitial cluster circled in Figure. 7 5 (B). The cluster forms an interstitial rocksalt lattice that is displaced with respect to the parent lattice. These ar e: ( i ) the conventional mechanism involving diffusion-controlled aggregation and annihilation of the originally introduced interstitials, and ( ii ) the spontaneous formation of new, structural FPs by the coordinated displacements of lattice atoms within the cluster core. The coordinated operation of the two mechanisms is illustrated in Figure 7 8 which captures four sequential hopping events (see arrows) during the growth

PAGE 129

129 and stabilization of a cluster containing initially four Mg2+ and five O2interstiti als together with one Mg and one O vacancy, Figure 7 8(a). (A cluster with this structure was also observed in the earlier simulations107, 110). Arrow 1 indicate s oxygen diffusion to a neighboring interstitial site while arrow 2 captures the annihilation of an oxygen FP by interstitial -vacancy recombination. Figure 7 8. Two distinct mechanisms control the growth of an interstitial MgO lattice with in the host lattice. The interstitial cluster A) contains four Mg2+ (orange) and five O2interstitials (blue) together with one Mg (green) and one O (red) vacancy. Black spheres represent Mg2+ and brown spheres O2 ions in the MgO parent lattice. B) The s pontaneous formation of a new oxygen FP within the cluster contributes to the growth of the cluster while simultaneously stabilizing it energetically. C) The s pontaneous formation of a second oxygen FP completes one cubic unit cell of the interstitial MgO lattice. A B C

PAGE 130

130 Both are conventional diffusion events of interstitials already present at the beginning of the kinetic phase. By contrast, Figure 7 8(b) reveals a novel mechanism involving the spontaneous creation of a new (not initially present) FP as a regul ar lattice O ion hops to an interstitial site (arrow 3), leaving behind a newly created vacancy within the cluster core; this event is then repeated in Figure 7 8(c) ; arrow 4 Two additional ions already present at the back of the interstitial cube in thi s figure act as nucleation sites for the completion of the next, adjacent cube. The net outcome of these events is a complete cube of an interstitial MgO crystal that has been formed and that continues to grow, on the displaced, interstitial lattice. 7.7 .1 Crystallography of the Cluster To elucidate the driving force for the spontaneous creation of new FPs in the cluster core, previous work on the nature of long range forces in rocksalt -structured crystals is recalled that has shown that a finite -sized cr ystallite is terminated such that the basic cubic building blocks remain intact.115 The reason for this is that any break up of these charge neutral, dipole -moment -free, octopolar (MgO)4 building blocks creates long ranged dipolar forces that significantly i ncrease the Coulomb energy of a finite -sized crystal.115 This analysis has been confirmed in numerous experiments116 and additional simulations117 that revealed that rock -salt structured surfaces, indeed, consist of complete octopolar building blocks (Figure 7 9), thereby eliminating any long range dipole moment perpendicular to the surface.115 The crystal structure of the interstitial cluster is same as that of the parent crystal structure, i.e., both have rocksalt crystal structure. In a building block (here (MgO)4) of a rocksalt crystal structure, every Mg ion is followed by an O ion and vice -versa (Figure 7 -

PAGE 131

131 10 (a)). Mg and O occupy 4a and 4b Wyckoff sites respect ively. Inside this building block, there is one empty octahedral site, Wyckoff 8c. Figure 7 9. A (MgO)4 octopolar. Decrease in the Coulombic energy is the driving force for the formation of a dipole -moment free octopolar cluster. It is this site that is occupied by the interstitial (say Mg, Figure 7 10(b)). If one of the crystallographic directions is translated to juxtapose next rocksalt building block, the next octahedral site is now occupied by an opposite ion, O (Figure 7 10(c)). Similarly, by building rocksalt crystal structure in three dimensions, it is found that the interstitial cluster is formed by alternating Mg and O ions that occupy alternating octahedral sites (Figure 7 9). Hence, interstitial cluster also maintains rocksalt crystal structure which is same as the parent crystal structure. This crystallography of interstitial cluster is unique for rocksalt -based materials. In UO2 (Chapter 8), it is found that interstitial crystal structure does not maintain the parent crystal structure d ue to the crystallographic limitation. Whereas in MgO, multiple unit cells are involved in clustering mechanism, and the cluster size is not limited, in UO2, the cluster is limited to only one unit cell.

PAGE 132

132 Figure 7 10. Crystallography of M gO rocksalt lattice. A) A (MgO)4 building block in which Mg occupy 4a and O occupy 4b Wyckoff sites. B) An octahedral interstitial position with Wyckoff site, 8c. C) In two MgO juxtaposed building blocks, alternate octahedral sit es are occupied by opposit e ions, Mg and O. 7.8 Discussion Applied to the observations in Figure 7 8 this means that the mechanism of formation of the interstitial MgO lattice should evolve by progressive formation of complete cubic (MgO)4 building blocks. In view of this, the c reation of the structural FPs appears to be simply a way to form complete octopolar building blocks, thereby lowering the total energy of the interstitial MgO crystallite. It is emphasized, however, that the FP creation mechanism is a cluster -stabilizing effect that does not happen A B C

PAGE 133

133 independently in perfect -crystal regions; it is therefore only indirectly related to the radiation -induced FPs from the ballistic phase. Similar behavior in UO2 is observed where new FPs are created in the oxygen sublattice; h owever, the details of clustering mechanism are different.118 Temperature accelerated dynamics (TAD) simulations107, 110 have extended the reach of PKA simulation s out to macroscopic time scales. However, PKA simulations of one or a few cascades correspond to low defect concentrations, thereby limiting the complexity of the defect structures that can evolve. By injecting a very high density of Frenkel -pair defects, the present method much more readily allows large clusters to develop. Moreover, this high density essentially allows nanosecond simulations to not only produce defect structures only previously seen in much longer simulations, but also produces heretofor e unreported defect structures. Because the initial defect concentration is extremely high, it might be argued that the dynamically generated complex defect structures seen here may in reality be rather uncommon. However, given that single defects and smal l clusters are annihilated as the system evolves, it is precisely these large, structurally complex clusters that will survive and constitute the long term radiation damage in the system. A similar methodology of sprinkling of defects has also been followe d in Monte Carlo (MC) simulation of radiation damage.119, 120 However, such MC simulations require prior knowledge of the specific atomic-level mechanisms and energetics. Such prior information is not required in these MD simulations; indeed, these simulations are capable of providing the mechanistic and energetic input required for the MC simulations.

PAGE 134

134 To conclude, by complementing the cascade simulations, this robust and straightforward approach enables elucidation of the kinetics of radiation damage developm ent unencumbered by the initial damage due to one specific cascade. By following only the kinetic evolution at elevated temperature, it enables observation of defect aggregation up to significant cluster sizes, and with mechanistic detail not previously po ssible. In addition, the ability to selectively create defects on only one of the sub -lattices provides an additional perspective on the nature of radiation damage development. Since it does not require particularly large system sizes, the approach thus pr ovides a relatively simple and direct route for elucidating the competing roles of point defect annihilation and aggregation in the development of radiation damage, including host lattice responses in the process.

PAGE 135

135 CHAPTER 8 RADIATION DAMAGE IN UO2 8.1 In troduction In Chapter 7, the capability of the new approach to perform radiation damage MD simulations in elucidating the evolution of defects in MgO is demonstrated. In this Chapter, same methodology is employed to elucidate point defect evolution in fluo rite based UO2. Fluorite based materials are widely used in the nuclear applications both as nuclear fuel and in the nuclear waste repository. In these applications, the irradiating highenergy particles cause damage by creating point defects, clustering o f which can form defects structures such as cluster, dislocations, voids and stacking faults, ultimately leading to material failure. The damage tolerance of a material not only depends on the number of defects created during irradiation but also on their long time kinetic evolution; understanding the complete evolution of defects is therefore very important. Computer simulations play a key role in complementing experiments by revealing important radiation -induced fast atomistic phenomenon106, 108, 99 121 that are otherwise impossible to capture. UO2 is a commonly used material in both nuclear applications in which radiation effects have been studied, for a long time, experimentally,122, 123, 124, 125 and more recently computationally.126, 127 The simulations have however been limited only to the initial collisional phase of damage; extending the time scale to elucidate computationally accessible complete defect -evolution mechanisms is therefore the next step. Under oxidizing conditions, UO2 readily accommodates O interstitials to form UO2+x. Ultimately, th ese interstitial s lead to a transformation to another distinct phase, U4O9,128 with a structure closely resembling that of fluori te A unit cell of U4O9 has

PAGE 136

136 space group of I 4 3d and consists of a cluster of oxygen interstitials and vacancies known as cuboctahedral (COT) cluster (Figure 8 1).129, 130 In this structure, U atoms occupy positions very close those in the f luorite positions of UO2 while O atoms form a COT cluster inside the unit cell. In comparison to eight oxygen atoms inside a unit cell of fluorite UO2, U4O9 has twelve oxygen atoms situated on the 12(b) position of I 4 3d Fi gure 8 1. U4O9 crystal structure. In this structure, U atoms occupy same positions as in fluorite UO2. O atoms form a cuboctahedral cluster by occupying 12(b) positions in space group I 4 3d (Reproduced from reference 130) Under irradiation of UO2, there can be a local increase in the oxygen interstitial concentration near the damage site. Under the se conditions it is possible that cluster formation might cause a local phase transformation, changing the cr ystal structure and thereby affecting the radiation tolerance. In this work, MD simulations are used to elucidate the evolution of irradiation induced point defects in UO2, including the formation mechanism of COT cluster s in the fluorite matrix

PAGE 137

137 8.2 Sim ulation Methodology While applying the method mentioned in Chapter 7 to MgO,131 it was able to reproduce the clusters that were observed from the combined radiation damage MD TAD simulations.110 Moreover, by focusing entirely on the kinetic phase of radiation damage, this method allowed to better understand the defect recombination and clustering mechanism in MgO. In particular because the defects are created by hand, it allows to choose the sub-lattice to create the defects a nd therefore, to elucidate for example, the defect evolution on one sublattice in the absence of defects on the other sublattice 8.2.1 Interatomic Potential and Molecular Dynamics Method The interactions between the ions are described by a pair -potenti al. Here, the short range atomic interactions are captured using the potential form provided by BushingIda et al.132 with parameters taken from Basak et al.133 In addition to the widely implemented Buckingham -type form for the ionic materials, the Bushing -Ida potential form also includes a Morse term that introduces a covalent component. Due to the assumption of partial covalency, the ionic charges are given non -formal values, which take into account partial charge transfer between the ions. The melting poi nt for this interatomic description of UO2 has been determined to be ~ 3450 50 K,134 which is in fair agreement with the experimental value of 3100 K135, 136. Likewise, the oxygen sub lattice ordering temperature for this potential has been determined to be 2200 K134 which is also in fair agreement with the experimentally obtained value of 2600 K.137 Throughout this work, all simulations have been carried out at a temperature of 1000 K, which is well below the melting and sublattic e melting temperatures. At this temperature, the oxygen interstitials and vacancies have high diffusivities, whereas the uranium counterparts have very high migration energies and

PAGE 138

138 therefore are less mobile; this temperature is also representative of the temperature range of 800 1600 K present in a typical fuel pellet. In order to elucidate the interplay of U and O point defects at elevated temperature, all simulations are carried out on UO2 single crystal containing 96,000 ions. An MD time step of 0.5 fs is used which is small enough to conserve energy for several thousand MD steps in NVE test runs. To capture as much of the defect evolution as possible, all the simulations have been carried out for more than 1 ns. The point defects, i.e., vacancies and in terstitials, are identified by comparing the structure with defects at any time with the initial structure containing no defects. If any lattice site does not contain an atom within a spatial cut off radius of 1 the site is identified as a vacancy. 8.2 .2 Diffusion -controlled K inetic Evolution of Defects From multiple simulations, it is observed that the defect environment on one sublattice significantly influences the point defect evolution on the other. To illustrate this, point defects are randomly c reated either on one sub lattice or both, and let the system equilibrate at 1000 K. A ll the vacancies of each type are separated by at least 5ao from their counterpart interstitials. In the following section, multiple simulation results distinguished by the defects on one or both sub lattices are discussed. Our first simulation was performed on a system containing no defects. The system was equilibrated for a sufficiently long time at 1000 K. As expected, t he analysis of the structure shows no defects at an y stage. This structure is used as the starting point for the simulations that are discussed in detail here.

PAGE 139

139 8.2.3 Frenkel pair Defects on a Single Sublattice In order to provide a benchmark against which to compare the results, first focus is on the evol ution of the system when the FPs are created on only one sub-lattice: either the O sub lattice or the U sub -lattice. Taking the oxygen case first, FPs are created by randomly picking an O atom from its cubic sub lattice site and placing it at an octahedra l interstitial site. At the beginning of the simulation, 200 FPs are created on the O sublattice on the previously equilibrated system. A snapshot of the initial structure with O FPs is shown in Figure 8 2(a). The system is then allowed to thermally equil ibrate at 1000 K. Figure 8 2. Evolution of the point defects present only on the oxygen sub lattice: red for vacancies, blue for interstitials Snapshots taken at A) t = 0 ps, B) t = 40 ps and C) t = 300 ps. Initially 200 FPs were randomly distrib uted on the oxygen sublattice. As time progresses, after 40ps and 300ps, only 78 and 21 O FPs are left behind. The O FPs annihilate by a vacancy intersititial recombination mechanism with no clustering of the defects. The snapshots of the subsequent tim e evolution of the O FPs to t = 40 ps and 300 ps are shown in Figures. 8 2 (b) and (c) respectively. Due to the low migration energies of both the O interstitials and the vacancies, the O FPs have high diffusivities. As evident from the Figures 8 2(b) and (c), the number of O FPs decreases with increasing time. The number of O FPs are shown as red triangles in A B C

PAGE 140

140 Figure 8 3. By 1400 ps, only four O FPs are left in the system. Eventually, all of the O FPs recombine by vacancy-interstitial recombination mechani sm; no interstitial interstitial or vacancy vacancy clustering is observed. This demonstrates that, in the absence of defects on U (or cation) sub lattice, the O sub lattice is completely healed within a few picoseconds, resulting in no longlasting damage to the material. Figure 8 3. Total number of O FPs present when the system is initialized with three different defect conditions: (1) defects present only on O sublattice (red triangles); (2) defects present only on U sublattice (green diamo nds) and (3) defects present on both O and U sub lattices (blue circles); In contrast to (1), in other two cases, new FPs are created in O sub lattice when defects are also present on U sub lattice. In contrast to the above expected result, the simulation of FPs on the U sub lattice alone yield a surprise. 200 FPs are created only on U sublattice. The initial snapshot of the system with defects is shown in Figure 8 4 (a); snapshots at later times are shown in Figures 8 4 (b) and 8 4 (c). Since, the U inter stitials and vacancies have very high migration energies, recombination events are less frequent. Out of the initial 200 U FPs,

PAGE 141

141 160 U FPs are still present after 1700 ps. This is as expected based on the high migration energies of U vacancies and interstit ials. What was not expected, however, was that the high concentration of the defects on the U sublattice nucleates O FPs. This is shown as the diamonds in Figure 8 3. Starting from no O FP at t = 0 ps, there is a gradual increase in the number of O FPs; t his peaks at ~600ps with a total of ~30FPs, and then slowly decreases; however even after 1700ps about 20 O FPS remain ( Figure 8 4 ( b )). Figure 8 4. Evolution of the point defects present only on the uranium sub -lattice. A) t = 0 ps, initially 20 0 FPs were randomly distributed on the uranium sub -lattice. The U vacancies and interstitials are less mobile at this temperature. Therefore, only few recombination events take place. B) At t = 1000 ps, due to the large number of U defects, some new O FPs are created which diffuse with time. C) At t = 1700 ps, the interstitials diffuse and form a cluster shown as circled. (See Figure 8 -7 for details on cluster). Although most of vacancies and interstitials that make up these O FPs remain close to each other (within 2ao), some of them diffuse far enough from their counterpart to suggest their potential for longterm stability. Indeed, some of the defects are no longer isolated but, as seen in Figure 8 4 (c), actually form clusters. A detailed analysis of thes e structures is given below. What is most noteworthy here, however, is the observation that the presence of the defects on the U sub lattice leads to the spontaneous formation of O FPs, whose diffusion controls the longer time scale evolution of the syste m. A B C

PAGE 142

142 8.2.4 FPs on Both U and O Sub -lattices The above simulations are rather artificial as a material under irradiation always creates defects on both sublattices simultaneously. Attention is now turned to a more realistic radiation -damage situation in which the defects are simultaneously present on both sublattices. In this simulation, 200 FPs defects are simultaneously created on both U and O sub lattices. As in the previous simulations, all of the interstitials are created on the empty octahedral sites of the fluorite unit cell. A snapshot of the random defect distribution in the system at t = 0 ps is shown in Figure 8 5 (a). Figure 8 5. Evolution of the point defects simultaneously present on both oxygen and uranium sublattices: orange and bl ack for uranium vacancies and interstitials. Snapshots taken at A) t = 0 ps, B) t = 40 ps and C) t = 1000 ps. Initially 200 FPs were randomly distributed on each sub lattice. As time progresses, most of the O interstitials aggregate to form clusters. The U vacancies and interstitials are almost immobile and do not participate in clustering directly. During equilibration at 1000K temperature, the U interstitials and vacancies are less mobile, as expected based on the simulations on the U sub-lattice alone, with only a few diffusing at all through the lifetime of the simulation (within first 1000 ps, only 15 U FPs have recombined) Most of them remain at their initial positions. Snapshots of the evolving system at t = 40 ps and 1000 ps are shown in Figures 8 5 (b) and 8 -5 (c).In contrast to the results for FPs on the O sub -lattice alone, a substantial number of defects A B C

PAGE 143

143 remain at the end of the simulation. The time evolution of the effective O FPs concentration shows a large initial increase within the first 100 150 ps, followed by a slow decline. Instead of almost complete elimination of the O FPs at the end of the simulation, their number is actually larger than at the beginning of the simulation. This indicates that there are two different phenomenon taking place simultaneously. First, in few cases, random diffusion of the vacancies and interstitials results in their elimination by recombination. Second, more interesting, and responsible for the sharp increase, is the formation of new O FPs. These are created by an O ion spontaneously displacing from the regular lattice site to a nearby interstitial site. The un recombined original O FPs and the newly created O FPs together evolve into clusters: while the O interstitials form COT clusters, the O vacancies for m Schottky defects with the U vacancies. The structures of these clusters are discussed in detail below. Considered together, these three simulations illustrate two important points: First, U FPs influence the radiation tolerance of UO2. The presence of U FPs not only prevents O FPs recombination, but also in certain cases contributes to creation of new ones. This in turn shows that radiation tolerance of UO2 decreases if the U FPs remain in the system. Second, from the methodology viewpoint, this approach of selectively creating defects on different sub -lattices, allows better understand the defect evolution mechanisms. 8.2.5 Vacancy Clustering: Schottky Defects A Schottky defect consists of two O vacancies and one U vacancy. The evolution of vacancies in to Schottky defects takes place in both the system with defects only on U sub -lattice and in the system with defects on both O and U sublattices. The O vacancies for the vacancy -cluster formation in the system with defects on both sub-lattices comes mainl y from the O FPs created at the beginning of the simulation. For the system with

PAGE 144

144 only U defects, the dynamically created O FPs contribute O vacancies to the Schottky defect. In both cases the process involves the diffusion of the mobile O vacancies to esse ntially stationary U vacancies. A snapshot of a Schottky defect is shown in Figure 8 5. Separate static calculations have shown that O -U O vacancy cluster has a much lower energy than O -U di -cluster with a separate O vacancy. Figure 8 6. Schottky-d efect formation. Two O vacancies tend to align in <111> direction around a U vacancy. The system with completely separated vacancies has a higher energy yet. As a result, although the Schottky defect is created by diffusion, it is extremely stable. All of Schottky defects involve vacancies that are near neighbors of each other. While most of the Schottky defects form as VO-VU-VO along <111> ( Figure 8 5 ), some form along <110>. Regardless of the specific structure, the creation of these Schottky defects sequ esters O vacancies, thereby making them unavailable for vacancy interstitial recombination events. As it will be seen, this leads to the possibility of the O interstitials, whose potential recombination partners are in Schottky defects, forming complex def ect structures.

PAGE 145

145 8.2.6 Interstitial Clustering: Cuboctahedral Clusters With the O vacancies segregated into the vacancy clusters, the O interstitials also evolve into clusters. The segregation of the interstitials into clusters is a rapid process and the va cancy and interstitial clustering occur simultaneously, provided the defects are available. In the system with defects on both sub lattices, the interstitial clusters begin to appear within first 40 ps. On the other hand, in the system with defects only on U sub lattice, it takes more than 1 ns before the newly formed O interstitials segregate and form clusters. The O interstitials form COT clusters in both systems. These COT clusters contain twelve O interstitials and eight O vacancies inside a unit cell of UO2. Such interstitial clusters have been previously observed experimentally in the oxygenrich UO2 systems129, 130, 138 and have been found stable from electron ic -structure calculations.139 These simulations allow the dynamics of the formation process itself to be determined. Figure 8 7. T hree different kinds of cuboctahedral (COT) clusters as seen in Figure 8 5 differentiated by filled or empty octahedra l site. The cubic oxygen sub-lattice has 4 O interstitials in addition to 8 original atoms. These extra interstitials displace the lattice oxygen atoms thereby creating, in total, 12 O interstitials and 8 O vacancies. The U lattice remains undisturbed. A) An empty COT cluster, B) O -filled COT cluster and C) U -filled COT cluster. A B C

PAGE 146

146 There are two different COTs, distinguished by the vacant (COT v ) or occupied (COT o ) octahedral site at their center. In the COT o clusters, t he octahedral site is occupied by an oxygen ion. In the simulations, in addition to these two previously described interstitial clusters, it is found that a U ion can also occupy the octahedral site inside the COT cluster (COT -u ). Figures 8 7 (a), (b) and (c) show the COT -v COT o and COT u c lusters, as spontaneously formed in the simulations. In all of these clusters, while the U FCC lattice (green) remains undisturbed. By contrast the concept of an O sub-lattice begins to break down. In particular, rather than 8 O ions occupying the corners of cubic sublattice, there are 12 O interstitial ions and no ions at all occupying the sub-lattice sites. It would seem extremely unlikely that all of the 8 original O ions were displaced from their lattice sites. What actually happens is the addition o f extra 4 O interstitials to a unit cell that previously had an un -defected O sub lattice, leads to the 8 ions on the O sub-lattice being displaced from their positions. The result is that all 12 O ions finally reside in the interstitial positions. A crystallographic analysis of the COT clusters is given elsewhere.129, 130 The formation of such clusters would be extremely unlikely if it required all of the O ions to converge on a single unit cell at the same time. Indeed, the simulations show that the COT forms in a progressive manner, as illustrated in Figure 8 8 for the COT -u cluster; the COT -v and COT -o clusters form in a similar manner. Figure 8 8 (a) is a sn apshot at t = 3ps of a unit cell containing a U interstitial at its center. At this time, all the O sub lattice sites are occupied. At t = 4 ps, one diffusing O interstitial enters the unit cell; as a result two O ions move from their lattice sites to inte rstitial sites, creating the structure shown in (Figure 8 8 (b)). The unit cell now contains two O vacancies and three

PAGE 147

147 O interstitials. In Figure 8 8 (c) at t = 20 ps, another O interstitial enters the unit cell, which creates two further O FPs. As a resul t, there are four new O FPs created in the unit cell by the addition of two O interstitials. Similarly, further addition of two more O interstitials creates four more O FPs. Figure 8 8. P rogressive formation cuboctahedral (COT) cluster, as seen in an MD simulation. A) A fluorite unit cell containing U interstitial atom at the octahedral site at t = 3 ps. B) At t = 4 ps, a diffusing O interstitial enters the unit cell. Two oxygen atoms are knocked from their lattice site resulting in forma tion of two O vacancies and three O interstitials. C), D), and E) are snapshots at t = 20, 120 and 250 ps. For every additional O interstitial entering the unit cell, two oxygen atoms are knocked off from their lattice sites. F) The relative decrease in th e energy of the system as the cluster is formed. Thus, every addition of one O interstitial in the unit cell displaces two O atoms from their lattice sites thereby creating two new O FPs (shown in Figures 8 -8 (d) and 8 8 A B C D E F

PAGE 148

148 (e)). Figure 8 8 (f) shows the ener getics of the cluster formation. It illustrates the system favors cluster formation instead of separated O interstitials. The energy of the system decreases when the distributed O interstitials segregate and progressively form the COT cluster. The other tw o clusters, COT v and COT o, follow the same clustering mechanism. This formation mechanism of a COT cluster therefore explains the sharp increase in the number of the O FPs as was observed in Figure 8 3. For every complete cluster that appears in the syst em, there is an increase of 8 O FPs, with all of these new FPs appearing from the O displacements within a cluster. It has to be kept in mind that the new FPs are formed only because of an extra oxygen ion entering the unit cell; their independent formatio n has not been observed. Figure 8 9. Clustering of oxygen interstitials in the absence of vacancies. Initially 4 U and 8 O interstitials were added to the system. The O interstitials diffuse and form cuboctahedral (COT) clusters. One cluster is complete while the second one is in the process of formation. The U interstitials (not shown here) remain immobile. From the results, it is believed that O interstitial -cluster formation is stabilized by formation of Schottky defects. In the absence of U v acancies, the O vacancies can

PAGE 149

149 recombine with the O interstitials. When U vacancies (and U interstitials) are present O vacancies can be pinned in Schottky -defect clusters This forces the O interstitials to either remain isolated or form some sort of cluster. To illustrate this directly, simulations on systems containing only interstitials or only vacancies are perfomed. One simulation contained additional stoichiometric U and O interstitials (four extra U and eight extra O interstitials). A snapshot of the system at t = 140 ps is shown in Figure 8 9 ( U interstitials not shown) The initially well -separated O interstitials diffuse and evolve into COT clusters. In Figure 8 9, one of the interstitial clusters is complete and the other one is in the process of construction. Figure 8 10. Snapshots showing clustering of U and O vacancies. A) At t = 0 ps, 4 U and 8 O vacancies are created in the system. B) O vacancies diffuse and cluster around U vacancies forming a <111> oriented cluster. Two <111> aligned clusters have formed by t = 280 ps. A similar scheme is applied to understand the evolution of the independent vacancies. Four U and eight O vacancies are created in the system. The snapshots of the initial structure at t = 0 ps and the equilibrat ed structure at t = 280 ps are shown in Figure 8 10. Initially, all the vacancies are created randomly and sufficiently far away from each other A B

PAGE 150

150 to prevent vacancy clustering. With time, the O vacancies segregate around the U vacancies to form Schottky def ects. In Figure 8 10, two such clusters have already formed and one is in the process of formation. 8.2.7 Evolution of Disproportionately Defected System s Up to this point, the point defect clustering process and the independent evolution of vacancies an d interstitials is illustrated. Here, the interdependence of the evolutions of O FPs on that of U FPs is presented. Figure 8 11. O FP evolution in a system initially containing 200 O FPs and fixed U FPs. After 500 ps, the number of O FPs levels off to ~ 60. Another 200 O FPs are added at 500 ps and 1000ps shown by two peaks. In both cases, the number of O FPs again level off to an equilibrium value. Evolved system at 1000ps was also separately equilibrated for 500ps. This evolution is shown with black triangles. Each case leads to almost the same number of O FPs at t=1500ps. The aim is to understand the O FPs evolution in a disproportionately defected system, i.e., the evolution of a system with a large number of O FPs in the presence of a low concentration of U FPs. To do this, 200 O FPs and only 50 U FPs are created. The evolution of the O FPs with time is shown in Figure 8 11. At 200 ps there is a significant

PAGE 151

151 decrease in the number of O FPs, which then levels off by about 500ps suggesting a c ertain concentration level of O FPs for U FPs. The system is perturbed with the hope of stimulating further interactions by creating a further 200 O FPs, as evidenced by the first peak in the plot. The system again equilibrates for another 500 ps. Most of the O FPs again recombine and annihilate leaving behind almost the same number of O FPs as were observed after first 500 ps. This shows that the additional 200 O FPs did not substantially disturb the evolution of the defects in the system and the O FPs con centration determined by the concentration of U FPs. At this point (1000 ps), two different simulations on the equilibrated system are performed. In the first, yet another 200 O FPs to the system are added again. In the second, the 1000 ps system is simp ly allowed to equilibrate independently without any addition of O FPs (evolution shown by black triangles). After an additional 500ps, the two systems have almost identical numbers of O FPs. Thus, it illustrates that even a frequent addition of the O FPs does not alter the final evolution of the system to a fixed concentration. Hence, it is concluded that the equilibrated O FPs concentration depends entirely on the U FPs concentration in the system. To confirm this conclusion further, six different simulat ions are performed with systems initially containing 100 O FPs but varying number of U FPs, i.e., 5, 10, 20, 30, 50 and 70. Figure 812 gives the time dependence of the number of O FPs for each system. In Figure 8 12, as the initial U FPs concentration is increased, the number of surviving O FPs also increases. For illustration, the systems containing 5 and 70 U FPs are compared. At t = 0 ps, both the systems contain 100 O FPs. However, with time, the number of O FPs differs significantly between two system s. For the system containing 5 U FPs, the O FPs number drops to 38

PAGE 152

152 by t = 100 ps and then to 25 at t = 400 ps. Due to a very low concentration of U FPs in this system, most of the O vacancies and interstitials annihilate by recombination and only few are l eft behind. The number of O FPs levels of by t = 400 ps and the data points for 400ps and 500ps overlap illustrating that no further O FP recombination takes place. Figure 8 12. N umber of O FPs (time progressing) as a function of number of initia lly created U FPs. Initially, 100 O FPs were created on each of six different systems containing 5, 10, 20, 30, 50 and 70 U FPs respectively. After every 100 ps time period, the number of O FPs are shown. There are higher number of O FPs alive after 500 ps in the system with initially higher U FPs. There are ~95 O FPs present in the system with initially 70 U FP (data on extreme right) in contrast to only ~30 O FPs in the system with 5 U FPs (data on extreme left). The system with 70 U FPs displays a differ ent evolution. In this case, first inspection shows that the data points for all six different time periods are concentrated very close to one another and near to the high end of the plot. This illustrates that not many recombination events take place. By t = 100ps, only 15 O FPs have recombined and 85 survive. Even with further progressing time, the number of FPs does not decrease; rather there is an increase from 85 at t = 100ps to 96 at t = 400ps. Unlike the recombination

PAGE 153

153 mechanism that was dominant in t he 5 U FPs (low concentration) system, here the system forms new FPs. The O FPs in the high U FPs system evolve into clusters. This shows that an increasing concentration of U FPs tends to stabilize an increased concentration of O FPs in the system. Furthe r, if the O FPs concentration is not sufficient for the given U FPs concentration, the system responds by creating new O FPs as also seen previously in Figure 8 5. The data for rest of the systems in Figure 812 also shows increased O FPs as the U FPs are increased. 8.3 Discussion The simulation approach used here complement s the conventional collision cascade simulations to provide a better understand the kinetic phase of defect evolution. An application of this method to radiation damage is previously pr esented in Chapter 7. Uberuaga et al.110 had used a combin ation of collision -cascade MD and TAD simulations to characterize the evolution of the point defects in MgO. The energy however, limited by the relatively small numbe r of defects created in a single cascade they could only observe very small clusters from multiple simulations. They found that the point defects left after the ballistic phase in the collision cascade evolved into interstitial clusters. On the other hand, this method allows to generate large concentration of the defects in a relatively smaller system. This method could not only reproduce their clusters but also observe much larger clusters. In addition, due to the formation of larger clusters, this method ology revealed the formation and growth mechanism of rocksalt clusters. It is interesting to compare this UO2 system with MgO. As in UO2, in MgO it is observed that if the defects are present only on one sublattice the FPs quickly recombine whereas, th ey cluster if the defects are present on both sub lattices. Further, not only is the defect evolution similar but in both case the clusters are also in part

PAGE 154

154 stabilized by formation of new FPs. In UO2, it is seen that new FPs are formed only on O sub -lattic e which stabilize the COT clusters. In MgO however, both Mg and O ions can displace from the lattice site to form new FPs. One difference arises from the difference in crystallography of the two systems. The interstitial lattice in MgO has the same structu re of parent rocksalt crystal structure T herefore, formation and growth of clusters is supported by new FPs that complete any missing interstitial in the (MgO)4 crystal structure. Since the crystal lattice and interstitial lattices have the same structur es, the extent and shape of the defect cluster is not limited. By contrast, the lattice of interstitial sites in UO2 does not have the symmetry of the parent lattice. As a result, the defect clusters have a defined extent and shape: a COT within a single u nit cell. The MgO interstitial clusters have been found to be highly stable110 even on experimental time scale; the stability of UO2 clusters beyond MD time scale is yet to be established. It is well known that the exper imentally observed COT clusters are found in oxygen rich UO2. It was shown by Murray et al.138 that the extra oxygen ions form COT clusters in systems with O -U ratio greater than 2.13. The observation of COT clusters in stoi chiometric UO2 with no extra O interstitials is rather unexpected. From the results, it is believed that it m ust be due to the high concentration of the interstitials in the system that act as extra O atoms. The hypothesis is supported by simulations on an only interstitial system (Figure 8 9). However, the absence of clustering in the only O FPs environment illustrates that charge imbalance on the U sublattice is also very important. From analysis of all the simulations above, a clear picture emerges th at the defects on the U sub lattice are primarily responsible for the radiation intolerance of the UO2.

PAGE 155

155 Anderson et al .140 have recently used DFT calculations to predict that split -quad interstitials are the most stable O interstitial configurations in the hyperstiochiometric systems, followed by COT clusters. The least stable are the individual interstitials at the octahedral site. The split quad interstitials are not observed in these simulations. Th is disagreement c ould be due to a lack of materials fidel ity in the empirical potential and could be due to the effects of temperature : the DFT calculations were performed at 0 K whereas the MD simulations are performed at very high temperatures. An atomistic simulation is only as good as the potential that des cribes the interatomic interactions. The Basak potential use d for these simulations has non -formal charges lower than the nominal valences of the ions. This introduces a degree of covalency to the system By contrast, the Busker potential uses strictly for mal charge model (i.e., U4+ and O2-) Corresponding simulations of FP evolution do not show the COT cluster formation. Further, with the Busker potential a pre assembled COT cluster quickly disintegrates. Why do these potentials give such different results and which is more physically reasonable? The difference clearly arises from the difference in charges. In particular, the electrostatic repulsion among the like -charged interstitials is much higher for formal charges (Busker) than for partial charges (Bas ak). That both DFT calculations and experiment have shown the presence of COTs s trongly suggests that, a non -formal charge model might be a better choice to model the defect evolution more appropriately. This conclusion is also consistent with the conclusi on of Govers et al .141 that partial charge models for UO2 generally provide higher materials fidelity than formal charge models.

PAGE 156

156 Finally, this study has been performed on single crystal elucidating the underlying atomic interaction mechanisms. A correspondi ng study in the polycrystalline UO2 is required to elucidate the point -defect evolution in the presence of the grain boundaries (GBs). In the polycrystalline system, not only will the point defects annihilate by vacancy -interstitial recombination mechanism but may also be eliminate d at the GBs. In that case, cluster formation might be impeded. An understanding of GB sink/source strength therefore needs to be established to understand the effect of GBs on defect evolution.

PAGE 157

157 CHAPTER 9 EVOLUTION OF POINT DEFECTS IN BCC MOLYBDENU M1 9.1 Introduction The nucleation, annihilation, diffusion, segregation and clustering of vacancy and self interstitial point defects in polycrystalline solids is of fundamental importance to many materials -science processes, such as c ompact sintering for microstructure development,142 irradiation damage in fuel and cladding components in nuclear systems,143, 144 and electromigration in electronic interconnects,145 to name just a few. As is well known, the equilibrium concentration and mobility of a particular type of point defect at any given temperature is determined by its formation and migration energies. However, under non-equilibrium conditions, the kinetic evolution of point -defect populations is driven by the competition between point -de fect production and annihilation at internal sinks or sources (such as grain boundaries (GBs), voids or dislocation cores) and/or external surfaces. Metals are widely used as structural materials inside nuclear reactors. Under these extreme conditions, they undergo damage due to irradiation. While damage due to irradiation in the bulk has been understood very well,121 146, 147 with some relation to grain boundaries (GBs),148 much work needs to be done to understand the behavior of grain boundaries (GBs) from their source and sink strength viewpoint. The GBs are known to act as sources and sinks for defects. However, their response under extreme conditions is not clear. For example, under normal radiation events, GBs while acting a s sources can emit interstitials to equilibrate high concentration of vacancies created in the grain. However, under extreme conditions, i.e., under high damage due to multiple radiation events at high temperature, GBs will also need to respond at the same rate by emitting 1In this chapter, all work was done in collaboration with Dr. Paul Millet at Idaho Nati onal Laboratory apart from Figures 9 3, 94 and part of 95, which are sole works of Dr. Millett. The nanocrystalline structures were also prepared by Dr. Millett.

PAGE 158

158 very high concentration of interstitials. However, intuitively, it occurs, that there must be a limit to which GBs can respond. This argument is also applicable to GBs while acting as sinks. Hence, a foremost question is: Do GBs saturate? In this work, a foundation is laid to elucidate the phenomenon by modeling a very simple crystal system, bcc Molybdenum. Damage due to radiation is modeled by using same methodology as described in Chapter 7. Here the results of molecular -dynamics (MD) s imulations are shown that focus on the dynamic nucleation and annihilation of both vacancies and interstitials at GB interfaces in nanocrystalline structures. The GB source and sink strength for vacancies as well as the sink strength for interstitials is e xtracted. The structures are initialized to be either over or under -saturated with respect to one particular point -defect specie (either vacancy or interstitial). Upon annealing at high temperatures, these point defects then migrate throughout the latti ce and annihilate at GB sinks (or nucleate from GB sources), thus allowing the characterization of the time -evolution of the system to equilibrium. Over -saturated point -defect systems can equilibrate by two mechanisms: (a) mutual vacancy -interstitial recom bination, and (b) elimination at GBs. The source/sink strength of GBs for vacancies and interstitials allows establishing the governing point -defect equilibration mechanism out of the two. Furthermore, the use of fully deterministic MD simulation in this w ork, as opposed to commonlyused kinetic Monte Carlo models which use a rigid lattice,149 is a f urther step toward s a more realistic characteriz ation of source/sink mechanisms. 9.2 Simulation Methodology In the present work, high temperature atomistic simul ations are performed on nanocrystalline structures of bodycentered cubic (bcc) Mo in order to investigate the

PAGE 159

159 grain -boundary source/sink strengths and diffusion kinetics of point defects at nanoscale dimensions. The three dimensional (3D) periodic simula tion cell contains six hexagonally -shaped, columnar grains of identical diameter (d = 20 nm) with a [100] columnar axis. The individual orientations of the grains in the global coordinate system are assigned to produce only high -energy, asymmetric tilt GB s, each with a misorientation angle of 30o. The structure of these GBs is continuously disordered along the interfacial plane. Further details on this nanocrystalline structure can be found elsewhere.150 The zero temperature MD input structures are heated up to the desired temperature in a step -wise fashion and all simulations are performed in the isothermal isobaric ensemble. The Finnis -Sinclair (FS) potential151 for bcc Mo is used throughout the study. The parameters that are of particular interest for dif fusion, the zero temperature relaxed vacancy -formation energy ( Ev f ) and the vacancy migration energy ( Ev m ), have been previously calculated as 2.54 eV and 1.32 eV,152 respectively. These are slightly below the experim ental values of 3.0 0.2 eV and 1.5 0.2 eV.153 The thermodynamic melting temperature for this potential is found to be Tm = 3097 10 K150 using the solid liquid interface velocity method.154 The highest annealing temperature chosen in this study is 2900 K, at which it was verified that GB pre -melting does not occur. 9.3 Vacancy and Interstitial Diffusivities Before considering nanaocrystalline microstrutures, the diffusivities of both vacancies and interstitials in the perfect crystal are determined. Single -crystal structures with 40 40 40 unit cell dimensions (128,000 atoms) and 3D periodic boundary conditions are isothermally annealed at temperatures ranging between 2500 K and 2900

PAGE 160

160 K. To determine the point de fect diffusivities, 20 vacancies or interstitials into the simulation cell are inserted and their self -diffusion is followed. A vacancy is created by randomly removing an atom from a regular lattice site, and, conversely, a self -interstitial is created by adding an atom to an interstitial site. It is furthermore ensure that, initially, every point defect is separated by a distance of at least 3ao from any other point defect. The self -diffusivity of the atoms DA, is determined from mean -squared displacem ent (MSD) calculations under zero -stress conditions using the Equation 2 19. From this, the vacancy diffusivity ( Dv) and interstitial diffusivity ( Di) can be obtained using the Equations 2 20 and 2 21. In Equations 220 and 2 21, fv and fi are correlation factors with values of 0.732155 and 0.300,156 respectively. For the 20 point defects present in each simulation, the concentrations of cv = ci = 20/128,000 = 1.56 104. Figure 9 1. Mean -squared displacement for single -crystal structures (40 40 40 ao) containing either 20 vacancies (blue squares) or 20 self interstitials (green circles) at T = 2900 K. The vacancy and interstitial diffusivities are calculated to be 0.021 nm2/ps and 0.081 nm2/ps, respectively. The MSD data for both structures is plott ed in Figure 9 1 for T = 2900 K. The least -squares fits to the data give vacancy and interstitial diffusivities of 0.021 nm2/ps and

PAGE 161

161 0.081 nm2/ps, respectively (2.1 108 and 8.1 108 m2/s) The diffusivities at each temperature are then plotted on an Arrh enius diagram for interstitials and vacancies (Figure 9 2), revealing activation energies of Eim = 0.49 eV and Evm = 1.06 eV. The migration energy obtained for vacancy diffusion is somewhat lower than the value of 1 .32 eV determined from the zero temperature calculations;152 this can perhaps be attributed to the thermal expansion of the lattice Figure 9 2. Arrhenius plot for Di and Dv in the temperature range bet ween 2500 K T 2900 K. The activation energies for interstitial and vacancy migration are found to be 0.49 eV and 1.06 eV, respectively. 9.4 Grain Boundaries as Sources and Sinks 9.4.1 Grain Boundaries as Vacancy Sources In order to investigate the GB source strength for vacancy emission, the nanocrystalline structures are initialized to be clearly under -saturated, with a vacancy concentration of cv o = 0, with respect to the thermal -equilibrium value. This is achieved by equilibratin g the GB structures at T = 2400 K (where GB diffusion is activated, however vacancy nucleation is not), and instantaneously stepping up the temperature to T

PAGE 162

162 = 2900 K. To investigate the vacancy nucleation and diffusion, t he sample is then annealed at 2900 K for approximately 1 ns, while the vacancy distributions are monitored at various times throughout the simulations. Figure 9 3. Progressive snapshots of the vacancy population within a single d = 20 nm grain at T = 2900 K for increasing time: A ) t = 0 ps, B) t = 72 ps, C) t = 205 ps, and D) t = 1005 ps. The microstructure is viewed along the [100] columnar axis. The concentration of vacancies, identified by blue dots, increases with time due to the on going nucleation from GB (brown atoms) sou rces. The time -evolution of the vacancy concentration, as quantified in the lower curve of Figure 9 5, follows an exponential approach to the equilibrium concentration of ~ 103. Figure 9 3 shows progressive snapshots in time of the vacancy population wit hin a single grain interior of the under -saturated structure. GB regions are visualized by the use of the common neighbor analysis,157 which identifies atoms in non -bcc crystalline environments (light brown dots). In addition, vacancies (represented by blue dots) are detected by determining lattice sites within the grain interiors that are not occupied by an atom. The viewing direction is along the columnar [100] axis. Initially, at t = 0 ps A B C D

PAGE 163

163 (Figure 9 3(a)), there are no vacancies. After 72 ps of annealing (Figure 9 3(b)), several vacancies have nucleated from the GBs, but have yet to diffuse extensively toward the grain centers. The vacancy concentration, at this time, has increased to cv = 2.9 104. As time continues (Figure 9 3(c), t = 205 ps), the va cancy concentration increases to cv = 4.9 104, and the vacancies have migrated further into the bulk. Finally, at time t = 1005 ps (Figure 9 3(d)), vacancies appear to be essentially evenly distributed throughout the microstructure with an increase in the concentration to cv = 9.3 104. Figure 9 4. Close up views of the nucleation of two lattice vacancies from a high -energy GB taken at A) t = 60 ps and B) t = 70 ps. The vacancies labeled 1 and 2 have nucleated from the GB and have migra ted towards the interiors of the grains they inhabit. In the latter snapshot, several new vacancies have been nucleated and still reside near the GBs. The emission of a vacancy from a GB source is a thermally activated process with a pre -exponential term that is associated with an attempt frequency.158 Therefore, a GB may produce many vacancies that reabsorb at their nucleation sites before a successful emission occurs in which a vacancy adequately distances itself into the lattice (similar to the creation of a stable Frenkel pair). Figure 9 4 provides a close up view of the successful emission of two vacancies during a 10 ps time interval, demonstrating that the GBs do indeed serve as sources for the increasing vacancy populations within the A B

PAGE 164

164 nanocrystallin e grains. Throughout the simulation, it is observed that many short lived vacancies directly adjacent to GBs (see Figure 9 4(b)) quickly annihilate back into the interface. A visual inspection reveals no qualitative change in GB structure or width throughout the simulation, suggesting that the source strength is constant and independent of past nucleation events. 9.4.2 Source -strength Analysis The equilibration of the vacancy concentration in the grain interiors due to emission from GB interfaces is diffus ion controlled. The most commonly used model for predicting the change in point -defect population due to microstructural source/sink mechanisms is the rate theory.144 According to this model, the time evolution of vacancies and interstitials can be expressed as: dcvdt G Rcvci kv 2Dv( cv cv eq) (9 1) dcidt G Rcvci ki 2Dici (9 2) According to Equations 9 1 and 9 2, the concentration rate of vacancies or interstitials depends on their production rate, G the mutual recombination rate, R and the microstructural source/sink strength, k2. cv and ci represent the concentration of vacancies and interstitials respectively, and Dv and Di represent vacancy and interstitial diffusivities respectively. In mater ial systems, certain equilibrium defect concentration lowers the free energy of the system by balancing the entropy components (Equation 2 13). In equation 9 1, cv eq represents equilibrium concentration of vacancies. The equilibrium con centration of interstitials is generally very low and can be neglected.

PAGE 165

165 In the absence of self -interstitials, and hence vacancy interstitial recombination, the time evolution of the vacancy concentration is given by: dcvdt kv 2Dvcv cv eq (9 3) In this continuum model kv 2 is assumed to be: ( i ) a weighted sum of the sink strengths for each microstructural element (GBs, dislocations, voids) and ( ii) spatially averaged over the entire material. In the work, GBs are the only sources in the simulation cell. Integrating Equation 9 3 and re arranging leads to cv( t ) cv eq cvo cv eq exp kv 2Dvt (9 4) where, c v o is the initial vacancy conce ntration. kv 2 is a geometrical parameter related to the average grain size (or, inversely, the density of GBs) with units of 1/length2.159 One can therefore also obtain the equilibration time, or the characteristic time for the system to evolve from cv o to c v eq as 1 kv 2Dv To quantify the evolution of cv in the system, the total number of vacancies in each grain is counted and divided by the number of lattice atoms in that grain. To reduce statistical noise, these values a re ensembleaveraged over all six constituent grains in the simulation cell (the error bars represent deviations in these data sets). The concentrations for both the under and over -saturated nanocrystals are plotted vs. time in Figure 9 5, revealing an asymptotic approach to cv eq = ( 1.0 0.1) 103. Regardless of the initial concentration, the vacancy distribution in the bulk regions thus approaches its equilibrium value within about a nanosecon d. Furthermore, it is pointed out that the equilibrium value is in qualitative agreement with the experimentally -measured vacancy

PAGE 166

166 concentrations in most bcc and fcc metals of 104 < cv eq < 103 just below the melting temperature.160 Figure 9 5. Evolution of the vacancy concentration in the under and over -saturated nanocrystalline structures at T = 2900 K (see Figures 93 and 9 6). In both cases cv evolves exponentially towards the equilibrium concentration of c v eq = (1 0.1) x 103. The dashed lines represent fits to Equation 9 4, yielding the GB sink and source strengths of kv 2 = 0.184 nm2 and kv 2 = 0.179 nm2, respectively. The approach to equilibrium is found t o be exponential, as evidenced by the fits of Equation 9 4 to the data, represented by the dashed lines. Inserting the value of Dv = 0.021 nm2/ps (see Figure 9 1) yields a GB source strength of kv 2 = 0.179 nm2. In addition, the chara cteristic equilibration time is found to be = 1/( kv 2 Dv) = 266 ps. It is found that this time corresponds to the instant at which cv has increased to 63% of the difference between cv eq and cv o (obtained from Equation 9 4) for t = )

PAGE 167

167 9.4.3 Grain Boundaries as Vacancy Sinks Grain boundaries, in addition to being vacancy sources, can also operate as sinks for both vacancies and self -interstitials. This is most relevant for particle irradiation i n nuclear materials, in which displacement cascades continuously produce excess concentrations of Frenkel pairs. Following a cascade, these defects will either cluster to form voids and interstitial loops or annihilate by mutual vacancyinterstitial recom bination or be eliminat ed at microstructural sinks (including dislocations, GBs, or free surfaces). In order to quantify the GB sink strength for vacancy defects, a high temperature anneal (at T = 2900 K) of the same nanocrystalline microstructure is pe rformed, however, with an initial super -saturation of vacancies created by randomly inserting an initial concentration of cv o = 2.0 103 (i.e., 2 cv eq ) into the grain interiors Figure 9 6 shows progressive snapshot s of the super -saturated microstructure (see also Figure 9 3). The initial structure (at t = 0 ps) in Figure 9 6 (a) contains a concentration of cv o = 2.0 103. As the annealing proceeds to t = 72 ps (Figure 9 6 (b)), t = 205 ps (Fi gure 9 6 (c)), and finally t = 1025 ps (Figure 9 6 (d)), the concentrations decrease to cv = 1.8 103, 1.6 103, and 1.2 103, respectively. In the absence of any interstitials in the system, t his decrease in vacancy concentration is solely due to a nnihilation at the GB sinks. The corresponding vacancy concentrations vs. time are shown in the upper half of Figure 9 5, revealing again an exponential approach to equilibrium (dashed line; see Equation 9 4). Remarkably, the evolution of the over saturat ed system is practically indistinguishable from that of the under -saturated system, albeit in an inverted manner. The equilibrium concentration and sink strength are found

PAGE 168

168 to be cv eq = ( 1 0.1) 103 and kv 2 = 0.1 84 nm2, which are nearly identical to the under saturated values. The equilibration time for this system is therefore = 1/( kv 2 Dv) = 259 ps (compared to 266 ps for the undersaturated system). Figure 9 6. Progressive sna pshots of the vacancy population in a single grain of the over saturated sample taken at A) t = 0 ps, B) t = 72 ps, C) t = 205 ps, and D) t = 1025 ps (see also Figure 9 3). The initial vacancy concentration is twice that at equilibrium. The evolution is quantified in the upper curve of Figure 95. The good correspondence between the values of GB source and sink strengths suggests that the GBs are equally efficient at nucleating and annihilating vacancies, and that the equilibration process is diffusion -co ntrolled. A B C D

PAGE 169

169 9.4.4 Grain Boundaries as Interstitial Sinks Similar to the above vacancy simulations, the GB sink strength of interstitials is elucidated by annealing a super -saturated interstitial system at T = 2400 K. This lower temperature is chosen due to t he fact that at 2900 K, vacancies nucleate from the GBs at an appreciable rate (see Figure 9 3) and recombine with the initially inserted interstitials, thus interfering with the GB sink -strength determination. Figure 9 7. Progressive snapshots of the oversaturated interstitial distribution within the d = 20 nm nanocrystalline structure. The sample was annealed at the lower temperature of T = 2400 K to avoid the GB nucleation of vacancies (that would interfere with the evolution of the interstitial s). The snapshots were taken at A) t = 0 ps where ci = 4.4 104 and B) t = 779 ps where ci = 1.1 105 and only a couple of interstitials remain in the grain interiors. At this lower temperature, vacancy nucleation is suppressed, at least for a short time (t < 800 ps) while interstitial diffusion dominates, thus allowing a determination of ki 2 Figure 9 7(a) shows the nanocrystalline structure prior to the simulation, containing an initial interstitial concentration of ci o = 4.4 104. After 779 ps of annealing (Figure 9 7(b)) only a couple of interstitials remain in the grain interiors corresponding to a concentration of ci ~ 105. The corresponding rate equation for the time evolution of A B

PAGE 170

170 interstitial concentratio n in a material with microstructural sources and sinks is generally given as144 dcidt ki 2Dici (9 5) Due to high formation energy, the equilibrium self interstitial concentration in metals is typically very low, even at high temperatures, and therefore is neglected (compare Equation 9 5 with Equation 9 3). Figure 9 8. Quantification of the time evolution of the self -interstitial co ncentration, ci, at T = 2400 K (see Figure 97). The dashed line represents a least -squares fit to Equation 9 6, yielding a GB sink strength for self interstitials of ki 2 = 0.146 nm2. Consistent with this no measurable equilibrium concentration in the simulations is observed (the final structure contains only a couple of interstitials among the six grains, which may still be too many). The time -dependent concentration is obtained by integrating Equation 9 5 to produce

PAGE 171

171 ci( t ) ci oexp ki 2Dit (9 6) As expected, this relation is similar to the vacancy evolution governed by Equation 9 4. Figure 9 8 shows this exponential approach of ci to zero. Fitting to Equation 9 6 (dashed line) y ields the GB sink strength for interstitials as ki 2 = 0.146 nm2. The equilibration time for this simulation is therefore = 1/( ki 2 Di) = 127 ps, using an interstitial diffusivity of Di = 0.054 nm2/ps (extrapolated f rom Figure 9 2(a)). Thus, the GB source strength for vacancies ( kv 2 = 0.179 nm2), the GB sink strength for vacancies ( kv 2 = 0.184 nm2), and the GB sink strength for interstitials ( ki 2 = 0.146 nm2) are all relatively similar 9.4.5 Discussion on GB Source/Sink Strength The vacancy and interstitial concentrations given in Equations 9 3 and 95, and those presented in the results above, are spatially averaged over the entire grain interiors (i.e ., they do not capture spatial gradients within the grain interiors). As a result, the GB source/sink strengths for vacancies ( kv 2 ) and interstitials ( ki 2 ) are simply geometrical constants that only depend on the gra in size and distribution. Brailsford and Bullough159 determined an analytical solution for these source/sink strengths by assuming a periodic array of grains (similar to the microstructures considered here). For the regions between the GBs, Ficks law was solved subject to the boundary condition that cv = cv eq at the GBs, and the vacancy concentration was spatially averaged over the domain throughout time. Their solution for GB source/sink stre ngth is kv 2 ki 2 15 rg2 (9 7) where, rg = d /2 is the grain radius. For d = 20 nm, this yields kv 2 = ki 2 = 0.15 nm2,. This

PAGE 172

172 value lies between the G B strength for vacancies ( kv 2 = 0.184 nm2 for sinks, and kv 2 = 0.179 nm2 for sources) and that for interstitials ( ki2 = 0.146 nm2). 9.5 Frenkel pair Annihilation Until now, the source/sink stre ngth of GBs for either point defects is elucidated (in the absence of another) using the third term in Equation 91 (or 9 2). The attention is now turned to the annihilation mechanisms in the simultaneous presence of both vacancies and interstitials. In Eq uation 91 (and 9 2), while the first term, the defect generation term, is constant throughout the simulations, the second term representing the annihilation by mutual recombination mechanism and the third them representing the annihilation at GBs are chan ging constantly. In this section, the effect of the last two terms is disentangled and the governing point -defect annihilation mechanism is determined. 9.5.1 FP Annihilation by Mutual Recombination The second term in Equation 9 1 depends on the concentrati on of point defects and the mutual recombination rate, R To extract R single -crystal simulation is adopted in which 50 FPs are randomly creased and their evolution is followed. Since the concentration of vacancies and interstitials, at any point of time is same, it can be assumed that: dcFPdt dcidt dcvdt Rcicv Rcv 2 (9 8) Time progressing snapshots of the structures at, t = 20 ps and t = 307 ps are shown in Figures 9 9(a) and 9 9(b) respectively. At any instant, an approaching vacancy and interstitial c an recombine and annihilate. Since there are no other sinks in a single crystal, the elimination is only by pure recombination. A time -dependent mutual recombination profile of FPs is shown in Figure 9 10. The analytical solution of the equation is a power

PAGE 173

173 law represented by cv = 1/[R(t) +1/ cv o] where cv o is the initial concentration of the vacancies (or FPs). Fitting with least -square method yields the recombination rate ( R ) as 42.58 ps1. Figure 9 9. Frenkel pair distribution in a single crystal at two different MD times, A) t = 20ps and B) t = 307ps. 50 FPs were randomly-distributed throughout the system separated by a minimum recombination radius of 2.0ao. Figure 9 10. Frenkel pair concentration vs time for the structure shown in Fig ure 9 9 representing the FPs recombination rate (R) in a single crystal. From the least -squares fit to Equation 9 8, R is found to be 42.58 1/ps. A B

PAGE 174

174 The recombination rate can also be related to the diffusion coefficient of the species by a simple relationsh ip:161 R 4r ( Dv Di) (9 19 where, r is the recombination radius (it is the radius within which interstitials and vacancies can corresponding values, the recombination rate of vacancies and interstitials obtained from Equation 9 9 is 41.61 ps1. This is in very good agreement with that obtained from the simulations (R = 42.58 ps1, Figure 9 10). The recombination radius, r can be backed out from Equation 9 9 by inserting R (= 42.58 ps1). Thus, a theoretically obtained value of recombination radius, r is 1.65 ao which is also a very reasonable value for crystalline systems. 9.5.2 Frenkel Pair Annihilation in a Nanocrystalline Structure Recapitulating the results, the source/sink strength of the GBs is calculated by elucidating the individual effect of third term in Equation 9 1. In this calculation, the second term is disregarded by allowing no recombination events. In the previous section, only the second term was focused by allowing recombination events in the absence of GBs. In this section, all processes are allowed to occur simultaneously, i.e., by inserting bot h vacancies and interstitials in the GBs structures. High concentration of point defects evolves to an equilibrium concentration by both mechanisms i.e., vacancy interstitial mutual recombination and elimination at the GBs represented by second and third t erms respectively. The individually obtained results from previous calculations will help in disentangling these two mechanisms when they occur simultaneously, thereby determining the dominant mechanism of point -defect annihilation.

PAGE 175

175 This simulation is star ted with a nanocrystalline GB equilibrated structured (Figure 9 6(d)). To allow both mechanisms to occur simultaneously, FPs of concentration, cFP = 1.2 103 are added on top of the initially present vacancies of concentration, cv = 1.2 103 The evol ution of the FPs at T = 2900 K is followed. Starting with a high concentration of FPs at t = 0 ps, shown in Figure 911 (a), the system evolves by FPs elimination as shown in Figure 9 11 (b). As expected, by t = 922 ps (Figure 911(b)), all the initially a dded FPs have eliminated and the system is left behind with the equilibrium concentration of vacancies. As discussed above, the vacancies and interstitials would have eliminated by two mechanisms. Figure 9 11. Evolving nanocrystalline structure s over -saturated with the FPs. Initially, the system is annealed at 2900 K to allow equilibration of the vacancies. Then, FPs are randomly added by maintaining a recombination radius of 2.0 ao. A) At t = 0ps, cFP = 1.2 103, B) t = 922ps, all FPs have re combined. In order to determine the dominant mechanism, evolution of FPs in the single crystal to those in the nanocrystalline structures are compared. If the rate of elimination of FPs in A B

PAGE 176

176 single crystal happened to be same as in GBs structure, it would il lustrate that mutual recombination is the dominant mechanism. Both systems are damaged with the same concentration of FPs. Figure 9 12 shows a comparison between the evolution of defects in two systems. Hollow circles represent data for the point defect s in nanocrystalline structures, whereas, filled circles with dotted line represent data for the point defects in single crystal. Figure 9 12. Comparison of the vacancy and interstitial concentrations throughout time for the nanocrystalline structure (hol low dots) and the single crystal (filled dots and dashed line). Initially, both the systems are equilibrated with the equilibrium concentration of the vacancies in addition to an excess concentration of FPs. As evidenced by the overlapping data sets, this plot illustrates that at this temperature and grain size, the effect of GB annihilation is negligible in comparison to FP recombination in the grain interiors. The evolution of the vacancies (top curves) and the interstitials (bottom curve) is represente d separately. The overlapping data for single crystal and nanocrystalline

PAGE 177

177 structures for both vacancies and interstitials illustrates that the rate of elimination of both the vacancies and interstitials is same. This elucidates that in both systems, mutual point defect annihilation mechanism is prevalent, and the point defects eliminate before reaching the GBs. In other words, GBs do not play significant role in the elimination of the point defects. 9.6 Conclusion In conclusion, using MD simulations, the p oint defect kinetics in model microstructures for nuclear metallic materials is elucidated. These simulations have illustrated that they are well -suited to determine the GB source/sink strengths for vacancies and interstitials in systems damaged due to irr adiation. They represent a first step towards more refined simulations that can offer atomic -scale insight into the accumulation, evolution, and /or mitigation of irradiation induced damage in nanostructured materials. By being able to capture atomistical ly, the formation of the vacancies at the GBs, the simulations have been able to disentangle the governing components in the rate theory models and have elucidated the dominant point -defect elimination mechanism. At this temperature and grain size, the eff ect of GB annihilation is negligible in comparison to FP recombination in the grain interiors. An obvious next step will be to understand the effect of the grain size and temperature on the defect elimination process. A crossover between the mutual recombi nation mechanism and GB elimination mechanism is expected. Other outstanding questions that can now be addressed include whether GB sinks can become exhausted after absorbing a high quantity of defects and the importance of defect fluxes in temperature gra dients.

PAGE 178

178 CHAPTER 10 SUMMARY AND FUTURE W ORK The evolution of point defects in fluorite -based materials, Bi2O3 and UO2, used in SOFCs and nuclear applications respectively has been elucidated. In addition, point defect evolution in rocksalt MgO and the role of grain boundaries in BCC Mo was also illustrated. 10.1 -Bi2O3 10.1.1 Cation Polarizability In Bi2O3, the intrinsically present oxygen vacancies play significant role in oxygen diffusion. Bi2O3 has a high concentration of vacancies and has high ionic polarizability, which sustain its high oxygen diffusivity. When doped with cations that have radius and polarizability different from Bi the diffusivity decreases. Experimentally, the effect of either is not clearly understood. Using MD simulations it was found that ionic polarizability is key for high oxygen diffusivity and that the ionic radius does not have significant effect. Furthermore, it was found that the cation polarizability alone has an effect on oxygen diffusion; the oxygen polarizabil ity does not effect diffusion. Moreover, a gradual increase in the cation polarizability increases oxygen diffusion. For a cation polarizability lower than a certain value, the vacancies order by aligning in combined <110> <111> directions. Once this vacancy -ordered structure is formed, oxygen diffusion cease s Hence, higher diffusivity can be achieved by using higher polarizable dopants. 10.1.2 Oxygen Diffusion Mechanism In a unit cell of -Bi2O3, there are multiple directions in which an oxygen atom can diffuse to occupy a vacant site. By following time averaged path of oxygen atoms, in

PAGE 179

179 agreement with experiments, it wa s found that oxygen diffusion preferentially occurs in <100> direction. However, in contrast to experiments that indicate oxygen passes th rough the empty octahedral site in the center of the unit cell, simulations showed that the oxygen follows predominantly a straight path. This apparent discrepancy could be due to different systems under study. While the experiments were performed on doped systems, simulations were performed on pure systems. A change in oxygen diffusion mechanism in the presence of polarizable ions might be expected. 10.1.3 Bonding Using DFT calculations, it wa s found that Bi atoms surrounded by <111> -vacancy ordering have all six Bi O bonds of same length, whereas, with those surrounded by <110> vacancy ordering have three pairs of different bond lengths. In <111> Bi, the electronic charge o n the Bi is equally distributed among surrounding oxygen atoms and there is no pref erential charge distribution. In contrast, <110> Bi has different charge distribution among six surrounding oxygen atoms revealed by different bond lengths. Furthermore, some of the <110> Bi charge is also directed towards the vacancies: this is the 6s2 l one pair charge. Different Bi O bonds also have different bond strengths. <110> Bi is predominantly ionic with valence charge of 2.84e whereas, <111> Bi has a covalent valence charge of 2.12e It was speculated that this directional and difficult to brea k covalent bond could be the limiting factor in oxygen diffusion. Whether this bond could be tailored by adding different dopants to enhance oxygen diffusivity remains to be seen. 10.1.4 Crystallographic Analysis of Structure Due to the ordering of vacanci es and the structural distortion that it undergoes, a 1x1x1 unit cell is no longer able to define the complete structure of Bi2O3. This

PAGE 180

180 structure has a lower symmetry than that of an ideal fluorite structure; the new space group is Fm 3 and the lattice parameter is 11.19761 which is twice that of the original unit cell. The structure is highly distorted with displacements shown by both anions and cations. The anions are displaced by 0.4 in a direction close to but not exactly <11 0>. Similarly, all <110> Bi cations are displaced by 0.3036 in a direction opposite to the plane of the <110> oriented vacancies. However, due to symmetric orientation of vacancies in <111> direction the <111> Bi cations are n ot displaced. 10.2 MgO, UO2 and Mo 10.2.1 Methodology To understand a complete evolution of point defects that are accessible on experimental time scale, a new methodology was designed that focuses on the long term kinetic evolution of the defects. While the conventional methodology focuses more on the point defect creation during initial phase of radiation damage, i.e., ballistic phase, this new methodology focuses on defect annihilation during the kinetic phase. To illustrate the capability of new methodology, MgO was chosen as a model material. It was shown that this methodology reproduces interstitial clusters observed from combined conventional MD TAD simulations, thus validating this methodology. 10.2.2 MgO Interstitial Clustering It was found that if the FPs are present only on one sub latice, i.e., either Mg or O, the corresponding vacancies and interstitials recombine and annihilate. In contrast, if the FPs are present simultaneously on both sub lattices, then instead of recombining, the interstitials segregate to form cluste rs. This provided a possible reasoning behind loss of radiation tolerance of MgO.

PAGE 181

181 It was further found that the clusters maintain the crystal symmetry of the parent lattice, i.e., the interstitial clusters have rocksalt crystal structure. During the cluste r growth, every alternating octahedral interstitial site is occupied by oppositely charged ion, forming a basic unit of (MgO)4. Furthermore, at any instance, if a basic unit lack an atom for completion, the lattice atoms respond by creating new FPs, thus s tabilizing the cluster. Hence, it was found for the first time that there are two distinct diffusion mechanisms for cluster growth, (a) conventional atomistic diffusion involving aggregation, and, (b) formation of new structural FPs. 10.2.3 UO2 Clustering While in MgO, interstitials are the only species that are involved in the cluster formation, in UO2, both interstitials and vacancies form clusters. There are some characteristic similarities in clustering processes in UO2 and MgO, despite the different crystal structures. Similar to MgO, in UO2, it was found that if FPs are present only on anion (O) sub -lattice, they recombined and annihilated. However, in contrast, if FPs are present only on cation (U) sub -lattice, unlike recombination in MgO, U FPs sur vive and nucleate new O FPs. The formation of new FPs was similar in MgO but the underling mechanisms are completely different. In MgO, new O FPs were formed to complete (MgO)4 clusters, whereas in UO2, new O FPs were formed perhaps to balance U defects. H owever, the complete mechanism is not known. I n both materials, if the defects are present simultaneously on both sub-lattices, very few recombination events take place. In UO2, while the O vacancies form Schottky defects by forming VO U I VO clusters, oxygen interstitials form cuboctahedral clusters (COT). The oxygen vacancies sequestration by uranium vacancies was found to be the mechanism of COT cluster

PAGE 182

182 formation. Hence, this illustrated that in the presence of defects on U sub lattice, ne w defects on O sub lattice are created, thereby elucidating that U defects are the ones that caused damage to UO2. The radiation tolerance of UO2 can be enhanced by controlling defects on U sub lattice. 10.2.4 Defect Annihilation in Nanocrystalline Mo As mentioned above, the point -defect evolution profile can change dramatically in the presence of GBs. A simple monotomic bcc Mo was chosen as a model material to elucidate the evolution of point defects in the presence of GBs. The MD simulations were found t o be able to capture the formation and elimination of defects from GBs, thus illustrating the source/sink nature of GBs. By inserting a high concentration of point defects, the source/sink strength of GBs was calculated. Further, while in single crystals, point defect elimination is only by mutual recombination, in the nanocrystalline structure, the point defects could also eliminate at the GBs. It was found that for the studied temperature and grain size, point defect elimination by mutual recombination is more prevalent and GBs do not have much significant role. 10.3 Conclusions and Future Work 10.3.1 -Bi2O3 As discussed in Chapter 4, doping with lanthanides decreases conductivity of Bi2O3 and results in the ordering of vacancies. In the -Bi2O3 literature, it appears that focus has been placed on improving the conductivity by engineering various str ategies such as double doping or varying dopant concentration.162 Little has been said about the site preference of dopants, which, is known to be one of the most significant factors in achieving high conductivity, especially in pyrochlore materials.163, 164 The crystal

PAGE 183

183 structure of pyrochlore materials closely resembles that of fluorite materials. Hence, it can be expected that dopants could prefer particular sites like in pyrochlores. In pyrochlore materials, for example, Gd2Ti2O7, if both cations occupy their regular sites, oxygen diffusion is limited even upon the availability of intrinsic oxygen vacancies. However, in the presence of cation antisites, oxygen diffusion is enhanced due to the creation of new vacancies, 48f which facilitate oxygen diffusion thr ough a 48 f hopping mechanism.163 While antisite formation and the associated hopping mechanisms are characteristic features of pyrochlore materials, the underlying mechanism of site preference, i.e., preference in <110> or < 111> vacancy ordered oxygen sub-lattice, could be applicable to Bi2O3. The site preference of dopants is not known, neither from experiments nor simulations. While -Bi2O3 has been excluded from such investigations, other fluorite -based materials have been studied both experimentally and theoretically. The current understanding is that the site preference strongly depends on the size of the cation. In bixbyite materials, it was found that those cation dopants with radii smaller than the host cation occupy 24b sites, whereas, those with radii larger than host occupy 8 a sites.165, 166 These two Wyckoff sites are the same as the <110> and <111> vacancy -ordered Bi environments respectively. Recently, in -Bi2O3, double doping has been used as viable option to achieve high conductivity while keeping dopant concentration to minimum.162 It remains to be seen, both in single -doped and double -doped Bi2O3, which sites the dopan ts prefer. Based on the above observation, one would expect dopants to occupy <110> vacancy ordered Bi site, as all dopants, except La, are smaller than Bi. In double doped system, even though the two dopants would compete for one site, one would still

PAGE 184

184 exp ect both dopants to be accommodated in <110>, in part because there are three times of these sites than <111> sites. This indicates that, regardless of dopant concentration, all dopants would occupy <110> sites. Keeping this in mind, it appears that in doped Bi2O3 experiments, the observation of loss in conductivity corresponds to materials with dopants sited at <110>. If this is indeed the case, it can be speculated that, with dopants occupying <110> and Bi occupying <111> (recall covalent character bonding), one of the ways that conductivity might be enhanced could be by breaking such symmetrical structure, i.e., by adding dopants in <111> to create structural distortions. Under such an assumption, dopants of size bigger than Bi, may hold the key in breaki ng such the symmetry by sitting at <111>. Computer simulations will be very useful to clarify this process. However, it appears that physically realizing this mechanism may not be straightforward. The Bi3+ cation has an ionic radius of 131 pm16 which is the biggest among all of the dopants of valence +3. To add larger sized dopants than Bi, one might have to look for options other than just the widely used +3 smaller lanthanide dopants. In addition to the above questions, there are still longstanding issues that are important to resolve in order to gain a complete understanding of Bi2O3. A preferential oxygen diffusion mechanism is known from experiments. However, it would be interesting to identify any difference in the vi cinity of dopants. Lanthanide dopants have lower polarizabilities than Bi; an altered diffusion direction thus may not be ruled out. Due to uncharged vacancies and charge -balanced cation site in doped Bi2O3, association of vacancies due to charge imbalan ce is ruled out. However, the difference in

PAGE 185

185 the size of the dopants to Bi would create stresses in the system. Therefore, the effect of stress on oxygen vacancies due to different dopant size can be understood in isolation. 10.3.2 Radiation Damage In thi s work on UO2 and MgO, evolution of irradiation -induced point defects was studied in single -crystal systems. The obvious next step is to understand their evolution in the presence of GBs and their coupling to microstructure. It can be expected that GBs wil l tend to maintain equilibrium concentration of point defects by acting as sources and sinks. Generally, while the concentration of interstitials is negligible due to their high formation energies, there is always a vacancy concentration in bulk. In UO2, t he U vacancy formation and migration energies are relatively higher than their O counterparts. For the temperatures of interest (1000 1600 K), U vacancy formation from GBs is not to be expected. By contrast, GBs can frequently form O vacancies and interst itials to maintain their equilibrium concentration. In irradiation induced point defects, it was shown in Chapter 8 that O vacancies do not annihilate if U defects survive; instead they form Schottky defects by combining with U vacancies. Under high concen tration of U point defects, O vacancies and interstitials easily exceed their equilibrium concentration. To maintain equilibrium concentration of defects, GBs may emit corresponding defects. Therefore, it is anticipated that three mechanisms could operate simultaneously, (a) formation of O interstitials from GBs to eliminate excess O vacancies, (b) drive for excess interstitials (clustered in COT) from bulk towards GBs. Through these mechanisms, GBs may drive for equilibrium concentration of O defects in the bulk. Depending upon the attempt frequency, GBs may also annihilate U defects, however, higher formation energies might limit that phenomenon. Meanwhile to maintain a high concentration of O defects, another mechanism may also prevail, i.e., (c) U defect s

PAGE 186

186 attempting to maintain fixed concentration of O defects. Hence, two competing mechanisms, point -defect annihilation driven by GBs and point -defect creation driven by U defects could prevail; understanding which one is dominant is essential. In MgO, the presence of GBs is also expected to affect the dynamics of the system. Diffusing Mg and O interstitials form interstitial clusters, whereas, their counterpart vacancies are essentially immobile. In the presence of GBs, in pursuit of maintaining equilibrium concentration of defects, it can be speculated that GBs may emit both Mg and O interstitials into bulk for reducing the concentration of Mg and O vacancies. However, at this point, it is difficult to foresee the fate of the already formed interstitial clus ters. Due to these speculations, it is very desirable to perform such studies, which are quite possible using computer simulations. The role of GBs can be further understood during frequent damage of material. Frequent damage would create large number of defects, to compensate for which, GBs would more frequently act as sources or sinks. It needs to be determined that under such circumstances, is the ability of GBs to act as source/sink is limited? In other words, do GBs become saturated with defects? Not all answers can be gained from experiments; therefore, simulations play a significant role in understanding radiation damage. Conventional radiation damage PKA MD simulation however will be incapable in unraveling all. Besides, the large system sizes requ ired, not only are PKA simulations highly time consuming, which limits observing such cluster formation, but they also require multiple trajectories of PKA to capture global defect behavior. They have not even been able to reveal complete point defect evol ution in single crystals, thus unraveling response of GBs alone from PKA

PAGE 187

187 simulations will be an enormous task. In contrast, by using the methodology described in this work, with the point defect clustering already understood from single crystals, such clus ters can now be easily adapted to polycrystalline materials. Hence, while the PKA simulation methodology appears to have limited value for single crystals, in polycrystalline materials, it may be completely useless. Apart from the above obvious questions, there are some other less understood problems, which can be understood in combination with this study. The stability of COT clusters is not known clearly. Although it was speculated previously that with higher temperatures, the COT clusters might disintegr ate,130 there is absolutely no understanding on their long -time stability. While they have been observed to be stable on MD time scale, their stability on experimentally accessible time scale needs to be established. The ef fect of non-stiochiometry is another standing issue. The evolution of point defects due to formation of U3 ,U5 in single crystals and later in polycrystals, needs to be understood. Further, there are stress and temperature gradients inside the nuclear reactor. A complete understating of evolution of point defects, or clusters, under gradients is also lacking. Similarly, deformation of the material under irradiation needs to be addressed. Deformation studies on un irradiated UO2 have been performed recently and creep mechanism has been understood.134 It will be interesting to observe any significant change in deformation mechanism in the radiated UO2. Finally, atomistic studies should be able to provide information to mesoscale level simulation methods such that evolution of large -sized defects such as voids and gas bubbles can be understood in the irradiated UO2.

PAGE 188

188 LIST OF REFERENCES 1 International Energy Agency (Key World Energy Statistics, IEA, Paris, 2006). 2 Ene rgy Information Administration (U.S. Department of Energy, International Energy Annual, Washington D.C., 2005). 3 V. S. Arunachalam and E. L. Fleischer, MRS Bulletin 33, 264 (2008). 4 M. S. Dresselhaus, G. W. Crabtree and M. V. Buchanan, MRS Bulletin 30, 5 18 (2005). 5 V. Vaitheeswaran, I. Carson, Zoom: The Global Race to Fuel the Car of the Future (Twelve Publishing, 2007). 6 D. M. Etheridge, L. P. Steele, R. L. Langenfelds and R. J. Francey, Journal of Geophysical Research -Atmospheres 101, 4115 (1996). 7 T he Intergovernmental Panel on Climate Change, (IPCC Fourth Assessment Report, AR4, IPCC, 2007). 8 G. A. Meehel, W. M. Washington, W. D. Collins, J. M. Arblaster, A. Hu, L. E. Buja, W. G. Strand and H. Teng, Science 307, 1769 (2005). 9 J. A. Turner, Scienc e 285, 687 (1999). 10 Intergovernmental Panel on Climate Change, Special Report on Carbon Dioxide Capture and Storage (Cambridge University Press, Cambridge, 2005). 11 L. B. Lave, MRS Bulletin 33, 291 (2008). 12 B. Hayman, J. W. Heinen and P. Brndsted, MRS Bulletin 33, 343 (2008). 13 D. Ginley, M. A. Green and R. Collins, MRS Bulletin 33, 355 (2008). 14 MRS Bulletin, Special Issue April 33, (2008). 15 Y. Chimi, A. Iwase, N Ishikawa, A. Kobiyama, T. Inami and S. Okuda, Journal of Nuclear Materials 297, 355 (2001). 16 M. W. Barsoum, Fundamentals of Ceramics (Institute of Physics Publishing, Philadelphia, 2003). 17 D. R. Askeland and P. P. Phule, The Science and Engineering of Materials, (Fourth Edition, Thompson, USA, 2003). 18 T. Hahn, International Tables for C rystallography Volume A (Springer, Fifth Edition, Boston, 2002).

PAGE 189

189 19 W. A. Goddard III, Handbook of Materials Modeling, Edited by S. Yip, 2707 (Springer, Netherlands, 2005). 20 J. Lemaitre, Journal of Engineering Materials and Technology Transactions of the ASME 107, 83 (1985). 21 J. W. Evans, Handbook of Materials Modeling, Edited by S. Yip, 1753 (Springer, Netherlands, 2005). 22 L. Q. Chen, Annual Reviews of Materials Research 32, 113 (2002). 23 Y. Mishin, Handbook of Materials Modeling, Edited by S. Yip, 459 (Springer, Netherlands, 2005). 24 J. D. Gale, Handbook of Materials Modeling, Edited by S. Yip, 479 (Springer, Netherlands, 2005). 25 A. D. Mackerell, Handbook of Materials Modeling, Edited by S. Yip, 509 (Springer, Netherlands, 2005). 26 C. J. Cramer, Essent ials of Computational Chemistry Theories and Models (Wiley, UK, 2005). 27 L. Verlet, Physical Review 159, 98 (1967). 28 R. W. Hockney, Methods in Computational Physics 9 136 (1970). 29 D. Potter, Computational Physics (Wiley, New York, 1972). 30 D. Beeman, Journal of Computational Physics 20, 130 (1976). 31 C. Gear, Numerical Initial Value Problems in Ordinary Differential Equation (Prentice Hall, Englewood Cliffs, NJ, 1971). 32 T. Schlick, Molecular Modeling and Simulation (Springer, Berlin, 2002) 33 M. E. Tu ckerman and G. J. Martyna, Journal of Physical Chemistry B 104, 159 (2000). 34 M. P. Allen and D. J. Tildesley, Computer Simulations of Liquids (Oxford University Press, Oxford, 1987). 35 J. Li, Handbook of Materials Modeling, Edited by S. Yip, 565 (Springer Netherlands, 2005). 36 C. A. Coulson, Nature 195, 744 (1962). 37 S. Khan, M. Islam and D. R. Basetb, Journal of Materials Chemistry 8 2299 (1998).

PAGE 190

190 38 D. Wolf, P. Keblinski, S. R. Phillpot, J. Eggebrecht, Journal of Chemical Physics 110, 8254 (1999). 39 J E. Lennard Jones, Proceedings of Royal Society 106, 441, 463, 709 (1924). 40 J. E. Lennard Jones, Proceedings of Royal Society 109, 476 (1925). 41 R. A. Buckingham, Proceedings of Royal Society A 168, 264 (1938). 42 F. London, Transactions of the Faraday Society 33, 8 (1937). 43 B. G. Dick and A. W. Overhauser, Physical Review 112, 90 (1958). 44 H. C. Andersen, Journal of Chemical Physics 72, 2384 (1977). 45 M. Parrinello and A. Rahman, Journal of Applied Physics 52, 7182 (1981). 46 M. Parrinello and A. Rahman Physical Review Letters 45, 1196 (1980). 47 M. Parrinello and A. Rahman, Journal of Chemical Physics 76, 2662 (1982). 48 S. Nos Molecular Physics 52, 255 (1984). 49 H. J. C. Berendsen, J. P. M. Postma, W. F. Vangunstteren, A. Dinola and J. R. Haak, Journal of Chemical Physics 81, 3684 (1984). 50 S. A. Adelman, J. D. Doll, Journal of Chemical Physics 64, 2375 (1976). 51 S. Nos, Journal of Chemical Physics 81, 511 (1984). 52 H. J. C. Berendsen, J. P. M. Postma and W. F. Vangunsteren, Journal of Chemical Phys ics 81, 3684 (1984). 53 J. Skinner and J. Kilner, Materials Today 3 30 (2003). 54 T. H. Etsell and S. N. Flengas, Chemical Reviews 70, 339 (1970). 55 N. J. Maskalick and C. C. Sun, Journal of Electrochemical Society 118, 1386 (1971). 56 H. Inaba and H. Taga wa, Solid State Ionics 83, 1 (1996). 57 S. Omar, Enhanced Ionic Conductivity of Cerai -Based Compounds for the electrolyte Applications in SOFCs (Dissertation, University of Florida, 2008). 58 J. W. Fergus, Journal of Power Sources 162, 30 (2006). 59 T. Takahashi, H. Iwahara and T. Arao, Journal of Applies Electrochemistry 5 187 (1975).

PAGE 191

191 60 D. M. Smyth, The Defect Chemistry of Metal Oxides (Oxford University Press, New York, 2000). 61 A. Bogicevic and C. Wolverton, Physical Review B 67, 024106 (2003). 62 N. Jiang E. D. Wachsman and S. H. Jung, Solid State Ionics 150, 347 (2002). 63 N. Sammes, G. A. Tompsett, H. Nfe and F. Aldinger, Journal of European Ceramics Society 19, 1901 (1999). 64 B. Aurivillius and L. G. Sillen, Nature 155, 305 (1945). 65 C. N. R. Rao, G. V. S. Rao and S. Ramdas, Journal of Physical Chemistry 73, 672 (1969). 66 N. Jiang and E. D. Wachsman, Journal of the American Ceramic Society 82, 3057 (1999). 67 E. D. Wachsman, S. Boyapati and N. Jiang, Ionics 7, 1 (2001). 68 P. W. M. Jacobs and D. A. Mac D naill, Solid State Ionics 23, 279 (1987). 69 R. D. Shannon, Acta Crystallographica A 32, 751 (1976). 70 K. Shirao, T. Iida, K. Fukushima and Y. Iwadate, Journal of Alloys and Compounds 281, 163 (1998). 71 R. D. Shannon, Journal of Applied Physics 70, 348 (1993). 72 S. Boyapati, E. D. Wachsman and B. C. Chakoumakos, Solid State Ionics, 138, 293 (2001). 73 S. Boyapati, E. D. Wachsman and N. Jiang, Solid State Ionics 140, 149 (2001). 74 B. T. M. Willis, Acta Crystallographica 18, 75 (1965). 75 E. D. Wachsman, S. Boy apati, M. J. Kauffman and N. Jiang, Journal of the American Ceramic Society 83, 1964 (2000). 76 A. Walsh, G. W. Watson, D. J. Payne, R. G. Edgell, J. H. Guo, P. Glans, T. Learmonth and K. E. Smith, Physical Review B 73, 2351041 (2006). 77 G. Kreese and J. H afner, Physical Review B 49, 14251 (1994). 78 G. Kreese and J. Furthmuller, Computational Materials Science 6, 15 (1996). 79 G. Kreese and J. Furthmuller, Physical Review B 54, 11169 (1996).

PAGE 192

192 80 J. P. Perdew and Y. Wang, Physical Review B 45, 13244 (1992). 81 P E. Blochl, Physical Review B 50, 17953 (1994). 82 H. A. Harwig and A. G. Gerards, Thermochemica Acta 28, 121 (1979). 83 N. I. Medvedeva, V. P. Zhukov, V. A. Gubanov, D. L. Novikov and B. M. Klein, Journal of Physics and Chemistry of Solids 57, 1243 (1996). 84 P. D. Battle, C. R. A. Catlow, J. Drennan and A. D. Murray, Journal of Physics C: Solid State Physics 16, L561 (1983). 85 G. Henkelman, A. Arnaldsson, and H. Jonsson, Computational Materials Science 36, 254 (2006). 86 E. Sanville, S. D. Kenny, R. Smith, a nd G. Henkelman, Journal of Computational Chemistry 28, 899 (2007). 87 A. Laarif and F. Theobald, Solid State Ionics 21, 183 (1986). 88 F. S. Galasso, Structure and Properties of Inorganic Solids (Pergamon, Oxford, 1970). 89 D. S. Aidhy, S. B. Sinnott, E. D. Wachsman, S. R. Phillpot and J. C. Nino, Journal of the American Ceramic Society 91, 2349 (2008). 90 P. D. Battle, C. R. A. Catlow, and L. M. Moroney, Journal of Solid State Chemistry 67, 42 (1987). 91 D. J. M Bevan and R. L. Martin, Journal of Solid State C hemistry 181, 2250 (2008). 92 P. D. Battle, C. R. A. Catlow, J. W. Heap, and L. M. Moroney, Journal of Solid State Chemistry 63, 8 (1986). 93 D. S. Aidhy, S. B. Sinnott, E. D. Wachsman, S. R. Phillpot and J. C. Nino, Journal of Solid State Chemistry, doi:10. 1016/j.jssc.2009.02.019, (2009) 94 E. Rutherford, Radioactivity in The Encyclopedia Britannica (The Encyclopedia Britannica Company, New York, Eleventh Edition, 1910). 95 World Association of Nuclear Operators, Performance Indicators (London, UK, 2006). 96 R. Ewing, W. Weber and J. Lian, Applied Physical Reviews 95, 5949 (2004). 97 R. Ewing, The Canadian Mineralogist 43, 2099 (2005). 98 I. Farhan, H. Cho and W. J. Weber, Nature 445, 190 (2007).

PAGE 193

193 99 K. E. Sickafus, L. Minervini, R. W. Grimes, J. A. Valdez, M. Ishimaru, F. Li, K. J. McClellan and T. Hartmann, Science 289, 748 (2000). 100 J. A. Angelo, Nuclear Technology (Greenwood Press, Connecticut, 2004). 101 D. Bodansky, Nuclear Energy Principles, Practices, and Prospects (Springer, New York, 2004). 102 Hj. Matz ke, P.G. Lucuta and T. Wiss, Nuclear Instruments and Methods in Physics Research Section B 166-167, 920 (2000). 103 W.J. Weber, R. C. Ewing, C. R. Catlow, T. D. de la Rubia, L. W. Hobbs, C. Kinoshita, Hj. Matzke, A. T. Motta, M. Nastasi, E. K. H. Salje, E. R Vance and S. J. Zinkle, Journal of Materials Research 13, 1434 (1998). 104 M. T. Robinson and I. M. Torrens, Physical Review B 9 5008 (1974). 105 R. Stoller, Journal of Materials Journal of Minerals Metals and Materials Society 48, 23 (1996). 106 T. D. de l a Rubia, H. M. Zbib, T. A. Kharaishi, B. D. Wirth, M. Victoria and M. J. Caturla, Nature 406, 24 (2000). 107 B.P. Uberuaga, R. Smith, A.R. Cleave, F. Montalenti, G. Henkelman, R.W. Grimes, A.F. Voter and K.E. Sickafus, Physical Review Letters 92, 115505 (200 4). 108 Y. Matsukawa and S.J. Zinkle, Science 318, 959 (2007) 109 D.J. Bacon and T.D. de la Rubia, Journal of Nuclear Materials 216, 275 (1994). 110 B.P. Uberuaga, R. Smith, A.R. Cleave, G. Henkelman, R.W. Grimes, A.F. Voter and K.E. Sickafus, Physical Review B 71, 104102 (2005). 111 L.V. Brutzel, J.M. Delaye, D. Ghaleb, Philosophical Magazine 21, 4083 (2003). 112 G. V. Lewis and C. R. A. Catlow, Journal of Physics C 1 8, 1149 (1985). 113 R. Ferneyhough, D. Fincham, G. D. Price and M. J. Gillan, Modelling and Simulat ion in Materials Science And Engineering 2 1101 (1994) 114 C. J. Meechan and J. A. Brinkman, Physical Review 103, 1193 (1956). 115 D. Wolf, Physical Review Letters 68, 3315 (1992). 116 M. Gajdaraziaka Josifovska, R. Plass, M. A. Schofield, D. R. Giese, R. Sha rma, Journal of Electron Microscopy 51, S13 (2002).

PAGE 194

194 117 F. Finocchi, A. Barbier, J. Jupille, C. Noguera, Physical Review Letters 92, 136101 (2004). 118 D. S. Aidhy, P.C. Millett, D. Wolf and S.R. Phillpot, (unpublished). 119 E. Vincent, C.S. Becquart, C. Domain Nuclear Instruments and Methods in Physics Research Section B 255, 78 (2007). 120 F. Soisson, Journal of Nuclear Materials 349, 235 (2005). 121 W. J. Phythian, R. E. Stoller, A. J. E. Foreman, A. F. Calder and D. J. Bacon, Journal of Nuclear Materials 223, 245 (1995). 122 D. R. Olander, Fundamental Aspects of Nuclear Reactor Fuel Elements (Technical Information Center, 1976). 123 Hj. Matzke and J. L. Whitton, Canadian Journal of Physics 44, 995 (1966). 124 Hj. Matzke, Radiation Effects 64, 3 (1982). 125 W. J. We ber, Radiation Effects and Defects in Solids 83, 145 (1984). 126 L. V. Brutzel and M. Rarivomanantsoa, Journal of Nuclear Materials 358, 209 (2006). 127 L. V. Brutzel and E. Vincent -Aublant, Journal of Nuclear Materials 377, 522 (2008). 128 B. T. M. Willis, Act a Crystallographica A34 88 (1978). 129 D. J. M. Bevan, I. E. Grey and B. T. M. Willis, Journal of Solid State Chemistry 61, 1 (1986). 130 R. I. Cooper and B. T. M. Willis, Acta Crystallographica A60 322 (2004). 131 D. S. Aidhy, P. C. Millett, D. Wolf, S. R. Ph illpot and H. Huang, Scripta Materialia, 60, 691 (2009). 132 Y. Ida, Physics of the Earth Planetary Interiors 13, 97 (1976). 133 C. B. Basak, A. K. Sengupta and H. S. Kamath, Journal of Alloys and Compounds 360, 210 (2003). 134 T. G. Desai, P. C. Millett and D. Wolf, Acta Materialia 56, 4489 (2008). 135 R. A. Hein and P. N. Flagella, GEMP 578. Schenectady: General Electric Company (1968). 136 J. P. Hiernaut, G. J. Hyland and C. Ronchi, International Journal of Thermophysics 14, 259 (1993).

PAGE 195

195 137 J. Ralph, J. Chem. Soc., Faraday Transactions II 87, 1253 (1987). 138 A. D. Murray and B. T. M. Willis, Journal of Solid State Chemistry 84, 52 (1990). 139 H. Y. Geng, Y. Chen and Y. Kaneta, Physical Review B 77, (2008) 180101 (R). 140 D. A. Anderson, J. Lezama, B. P.Uberuaga, C. Deo and S. D. Conradson, Physical Review B 79, 024110 (2009). 141 K. Govers, S. Lemehov, M. Hou and M. Verwerft, Journal of Nuclear Materials 366, 161 (2007). 142 R.M. German, Sintering Theory and Practice (John Wiley, New York, 1996). 143 I.M. Robertson, R.S. Averb ack, D.K. Tappin and L.E. Rehn, Journal of Nuclear Materials 216, 1 (1994). 144 N.V. Doan and G. Martin: Physical Review B 67, 134107 (2003). 145 C.A. Volkert, Encyclopedia of Materials: Science and Tech. (Elsevier, 2001). 146 D. J. Bacon, Journal of Nuclear Ma terials 323, 152 (2003). 147 Y. N. Osetsky, D. J. Bacon and A. Serra, Philosophical Magazine Letters 79 273 (1999). 148 M. Samaras, P. M. Derlet, H. Swygenhoven and M. Victoria, Philosophical Magazine 83, 3599 (2003). 149 L. Malerba, C.S. Becquart and C. Doma in, Journal of Nuclear Materials 360, 159 (2007). 150 P.C. Millett, T. Desai, V. Yamakov and D. Wolf, Acta Materialia 56, 3688 (2008). 151 G.L. Ackland and R. Thetford, Philosophical Magazine A 56, 15 (1987). 152 J.M. Harder and D.J. Bacon, Philosophical Magaz ine A 54, 651 (1986). 153 K. Maier, M. Peo, B. Saile, H.E. Schaefer and A. Seeger, Philosophical Magazine A 40, 701 (1979). 154 S.R. Phillpot, J.F. Lutsko, D. Wolf and S. Yip, Physical Review B 40 2831 (1989). 155 J. R. Manning, Diffusion Kinetics for Atoms in Crystals (Van Nostrand, Princeton, 1968). 156 N. Soneda and T. Diaz de la Rubia, Philosophical Magazine A 78, 995 (1998).

PAGE 196

196 157 J.D. Honeycutt and H.C. Andersen, Journal of Physical Chemistry 91 4950 (1987). 158 A. P. Sutton and R. W. Balluffi, Interfaces in Crys talline Materials (Clarendon Press, Oxford, 1995). 159 A.D. Brailsford and R. Bullough, Philosophical Transactions of The Royal Society of London Series A 302, 87 (1981). 160 R.W. Siegel, Journal of Nuclear Materials 69& 70, 117 (1978). 161 T.R. Waite, Physical Review 107, 463 (1957). 162 S. H. Jung, E. D. Wachsman and N. Jiang, Ionics 8 210 (2002). 163 P. J. Wilde and C. R. A Catlow, Solid State Ionics 112, 173 (1998). 164 L. Minervini, R. W. Grimes and K. E. Sickafus, Journal of the American Ceramic Society 83 1873 (2000). 165 C. R. Stanek, K. J. McClellan, B. P. Uberuaga, K. E. Sickafus, M. R. Levy and R. W. Grimes, Physical Review B 75, 134101 (2007). 166 C. R. Stanek, K. J. McClellan, M. R. Levy and R. W. Grimes, IEEE Transactions on Nuclear Science 55, 1492 (2008).

PAGE 197

197 BIOGRAPHICAL SKETCH Dilpuneet Singh Aidhy was born in 1982 at Rajpura, Pun jab, India. He received his B. E. degree in metallurgical engineering from Punjab Engineering College, Chandigarh, India in 2004. He joined the Department of Materials Science and Engineering at University of Florida in Fall 2004 to pursue his Ph.D. with P rof. Simon R. Phillpot. During his graduate studies, he interned at Idaho National Laboratory for one year where he worked with Dr. Dieter Wolf. Dilpuneet is expecting his PhD in the spring of 2009 and has accepted a post -doctoral position at the Departmen t of Materials Science and Engineering at Northwestern University.