Vacancies and Vacancy-Mediated Self Diffusion in Cr2O3: A First

Danny Broberg , Bharat Medasani , Nils E.R. Zimmermann , Guodong Yu , Andrew Canning , Maciej Haranczyk , Mark Asta , Geoffroy Hautier. Computer Physi...
0 downloads 0 Views 85MB Size
Subscriber access provided by UB + Fachbibliothek Chemie | (FU-Bibliothekssystem)

Article

Vacancies and Vacancy Mediated Self Diffusion in CrO: A First Principles Study 2

3

Bharat K Medasani, Maria L. Sushko, Kevin M. Rosso, Daniel K. Schreiber, and Stephen M. Bruemmer J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.7b00071 • Publication Date (Web): 05 Jan 2017 Downloaded from http://pubs.acs.org on January 9, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry C is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Vacancies and Vacancy Mediated Self Diffusion in Cr2O3: A First Principles Study Bharat Medasani,∗,† Maria L. Sushko,† Kevin M. Rosso,† Daniel K. Schreiber,‡ and Stephen M. Bruemmer‡ †Physical and Computational Sciences Directorate ‡Energy and Environment Directorate, Pacific Northwest National Laboratory, Richland, WA 99354, USA E-mail: [email protected]

1 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract Charged and neutral vacancies and vacancy mediated self diffusion in α-Cr2 O3 were investigated using first principles density functional theory (DFT) and periodic supercell formalism. The vacancy formation energies of charged defects were calculated using the electrostatic finite-size corrections to account for electrostatic interactions between supercells and the corrections for the bandgap underestimation in DFT. Calculations predict that neutral oxygen (O) vacancies are predominant in chromium (Cr)-rich conditions and Cr vacancies with -2 charge state are the dominant defects in O-rich conditions. The charge transition levels of both O and Cr vacancies are deep within the bandgap indicating the stability of these defects. Transport calculations indicate that vacancy mediated diffusion along the basal plane has lower energy barriers for both O and Cr ions. The most favorable vacancy mediated self diffusion processes correspond to the diffusion of Cr ion in Cr3+ charge state and O ion in O2− state, respectively. Our calculations reveal that Cr triple defects comprised of Cr in octahedral interstitial sites with two adjacent Cr vacancies along the c-axis have a lower formation energy compared to that of charged Cr vacancies. The formation of such triple defects facilitate Cr self diffusion along the c-axis.

INTRODUCTION Chromium is a critical alloying element for corrosion resistance of Fe- and Ni-base alloys in a wide variety of aqueous and gaseous environments. 1,2 In these alloys, Cr often forms a protective oxide coating in the form of Cr2 O3 , thereby minimizing chemical corrosion in the bulk alloy. Given that bulk Cr2 O3 is an anti-ferromagnetic semiconductor with a high N´eel temperature, it is also considered to be a promising material for spintronics and other magneto-electric applications. 3–5 Material degradation in corrosion-resistant, high-Cr alloys involves various mechanisms including localized corrosion/oxidation and stress corrosion cracking. 6–12 Self diffusion in

2 ACS Paragon Plus Environment

Page 2 of 51

Page 3 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Cr2 O3 plays an important role in the protective nature of the chromia coating, 13–15 by influencing the growth of the chromia scales as well as by controlling the supply of oxygen to the oxide/metal interface. Diffusion in oxides is usually mediated by point defects such as vacancies or interstitials. 16 Although self diffusion coefficients in Cr2 O3 under various conditions have been measured experimentally, 15,17–21 the diffusion coefficients reported in these studies span a wide range of values. For example at ∼1400◦ K, the reported diffusion coefficients range from 1 × 10−8 to 1 × 10−19 cm2 s−1 . In addition, the mechanism for diffusion is not well understood. Some studies concluded that Cr is more mobile than O, based on the study of the outward growth of Cr2 O3 scales, 15,22 while Sabioni et al. reported a higher diffusion coefficient for O. 19 The broad variability of the results from different experiments is often attributed to the morphology of the crystal, specifically to the presence of grain boundaries, which provide fast routes for diffusion. Another source for the variability is the presence of impurities in Cr2 O3 , in particular, divalent impurities such as Mg. These difficulties preclude conclusive mechanistic insights of defect mediated diffusion in chromia from experiments alone. Complementing the experimental studies, theoretical studies enable a fundamental understanding of the role of various native point defects in self diffusion processes and associated physical properties of metal oxides including their chemical, electrical, optical, and mechanical behaviors. 23 In addition self diffusion of native point defects influences the effectiveness of doping through self compensation 24 rendering the knowledge of self diffusion processes important in materials design. 4,5 To date few theoretical studies have been reported on point defects and diffusion in Cr2 O3 . 25–27 Catlow et al. and Vaari studied native defects and self-diffusion in Cr2 O3 using empirical pair-potential calculations. 25,26 Whereas Catlow et al. 25 studied the defect energetics and diffusion barriers using Mott-Littleton approach 28 and thermodynamic arguments, Vaari 26 modeled the diffusion coefficients at different defect concentrations using molecular dynamics. While providing important insights, molecular dynamics simulations

3 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

do not fully capture local atomic and electronic structure surrounding defects. To overcome these difficulties of classical approaches, first principles density functional theory (DFT) has emerged as a reliable theoretical method for materials modeling. 29,30 DFT with periodic supercell methodology provides rigorous approaches to study point defects and diffusion mechanisms. 31,32 Using DFT, Lebreau et al. 27 investigated neutral vacancies, Frenkel defects and the associated diffusion barrier energies in Cr2 O3 . However, defect charge state plays an important role in the formation energies of the defects, barrier energies, and the resulting activation energies for self diffusion. 31,33,34 Hence, a complete picture for self-diffusion processes in Cr2 O3 , which is a bandgap material, should include vacancies and interstitials in multiple charge states. The present study also encompasses charged vacancies, making the next step towards a comprehensive model for self diffusion in chromia. Using DFT+U level of theory, 35 we investigated Cr and O vacancies in the charge states ranging between [-3, 0] and [0, 2], respectively. The notation used in this work to represent the charge of a vacancy is different from the Kroger-Vink notation, and the correspondence between the two notations can be found in Table 1. Vacancy formation energy, Ef , an important parameter for the evaluation of diffusion coefficient and speciation of defects is computed using the periodic supercell methodology 31 within the thermodynamic formalism developed by Zhang and Northrup. 36 By evaluating the defect formation energies at compositional deviations from stoichiometric Cr2 O3 that are at the phase boundaries of Cr2 O3 , we estimated the possible range of vacancy formation energies. Special care was taken to correct for the intrinsic errors of periodic DFT+U in the calculation of formation energies for charged defects. The need for correction arises from the electrostatic interaction of the vacancy with its periodic images and with the neutralizing background charge. 31 Freysoldt et al. 37 proposed a posteriori correction method that was found to result in well converged defect formation energies even for small supercell sizes provided the defect charge was sufficiently localized. We used an extension to the correction method of Freysoldt et al. proposed by Kumagai and Oba 38 for anisotropic systems. Another

4 ACS Paragon Plus Environment

Page 4 of 51

Page 5 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

correction is needed due to the under-prediction of the bandgap by DFT+U. In this work we used the DFT+U correction scheme proposed by Janotti and Van de Walle 39 to correct the vacancy formation energies. This correction scheme was previously tested against the results obtained using a more accurate HSE06 hybrid functional 40 for defects in ZnO. The pathways for the diffusion processes are explored with climbing image nudged elastic band (cNEB) 41 to compute the barrier energies for all the charge states considered in this work. Using the corrected formation energies and barrier energies, we demonstrated the agreement between experimental and theoretical coefficients for O diffusion.

METHODS Density Functional Calculations All calculations in this study were performed within the DFT framework using the Vienna ab initio simulation package (VASP). 42–44 We used the projector augmented wave (PAW) method 45,46 and the PBEsol 47 generalized-gradient approximation (GGA) for the exchangecorrelation functional. The rotationally invariant DFT+U approach by Dudarev et al. 48 as implemented in VASP was used to treat the 3d electrons of Cr with an effective on-site Coulomb interaction parameter of Uef f = U − J = 3.7 eV. The choice of U parameter is discussed in detail in the Results section. A cut-off value of 400 eV was used for the planewave basis set in all calculations based on the convergence for bulk Cr2 O3 . Spin polarization with anti-ferromagnetic (AFM) ordering of Cr spins was used. Gaussian method with a width of 0.02 eV was used for electronic smearing. The unit cell Cr12 O18 with ferromagnetic ordering was the starting point for our calculations. It was obtained from the Materials Project 49 database and re-optimized with Gamma centered 6 × 6 × 3 grid and AFM ordering. The defect calculations in this study were performed within the supercell formalism. The atomic positions in the defect supercells were relaxed at constant volume until the individual forces on each atom were smaller than 5 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 51

0.01 eV/˚ A. Based on the convergence tests of defect energies with the supercell size, the 2 × 2 × 1 supercells (108 atoms) with 2 × 2 × 2 Monkhorst-Pack k -point grid were used for defect calculations.

Vacancy Formation Energies In the thermodynamic formalism proposed by Zhang and Northrup, 36 the formation energy of a vacancy of type X with charge q is given by

Ef (X q ) = Etot (X q ) − Etot [bulk] + µX + qEF + Ecorr ,

(1)

where Etot (X q ) is the total energy of defect supercell, Etot [bulk] is the total energy of bulk supercell, µX is the chemical potential contribution due to the removal of the atom, EF is the Fermi level and Ecorr is the correction term. In this work, finite-size electrostatic and bandgap corrections are considered in sequential manner. To account for the anisotropy of Cr2 O3 , we used an extension proposed by Kumagai and Oba 38 to the electrostatic correction scheme developed by Freysoldt et. al 37 to correct for the long range Coulomb interactions between finite size periodic supercells for charged defects. In this approach, the electrostatic corrections can be separated into point charge correction and potential alignment terms as

Ecorr = EP C − q∆VP C,q/b .

(2)

The point charge correction EP C , which represents the interaction energy of the defect with its periodic images, can be written as

EP C

1 = 2

Z (−VP C,q ) qδ(r)dr,

(3)



where VP C,q is the Madelung potential. The potential alignment term, ∆VP C,q/b , is the 6 ACS Paragon Plus Environment

Page 7 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

potential difference between the defect-induced potential, Vq/b = Vdef ect,q − Vbulk , and the Madelung potential, VP C,q and corresponds to the interaction energy of the defect with the neutralizing background charge, which is short range in nature. The potential alignment is computed by averaging over the potentials at atomic sites far from the defect. The electrostatic corrections were obtained with Python Charged Defects Toolkit (PyCDT), a software package developed by one of the authors. 50 Transition states are defined as thermodynamic transition level at which the defect changes its charge state. In particular, the transition state (q/q 0 ) between two defect charge states q and q 0 is given by 0

Ef (X q )|(EF =0) − Ef (X q )|(EF =0) . X (q/q ) = q0 − q 0

(4)

When the transition states are in either conduction or valence bands, then those states are considered to be in resonance with the corresponding bands. If the transition states are within the bandgap, but very close to either the valence band maximum (VBM) or conduction band minimum (CBM), the transition levels are classified as shallow states. If the transition levels are close to the middle of the bandgap, then they are considered as deep levels. Generalized gradient approximation (GGA) in DFT severely under-predicts the bandgaps due to the well known self-interaction error. 51 GGA+U partially corrects the bandgap error by minimizing the repulsion between the Cr-3d electrons with the core electrons. As a result, the valence band is shifted downwards and the conduction band is shifted upwards, thereby increasing the bandgap. Due to the underestimation of the Cr2 O3 bandgap by GGA and GGA+U, the gap between the defect electronic states in the bandgap and the VBM also tends to be underestimated. The error in the position of the defect states is greater when the defect states are closer to the CBM. When the defect states are occupied by electrons, the errors in the defect state positions introduce errors in defect formation energies. One

7 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 51

of the approaches for minimizing the error in the position of defect states is the so-called scissor operator method. 52 In this procedure the CBM is shifted upwards until the bandgap reaches the experimental value, and the defect states closer to the CBM are also rigidly shifted. When the defect energy states within the bandgap exhibit a mixture of conduction band and valence band characters, the amount of shift depends on the relative valence band to conduction band character 39 rendering the scissor operator correction insufficient. To further increase the accuracy in the position of transition levels predicted by GGA+U with respect to the VBM, we utilized a scheme proposed by Janotti and Van de Walle. 39 In this scheme, the modified defect transition levels (q/q 0 ) are defined as

(q/q 0 ) = (q/q 0 )GGA+U +

 ∆ Egexpt − EgGGA+U , ∆Eg

(5)

with (q/q 0 )GGA+U − (q/q 0 )GGA ∆ = , ∆Eg EgGGA+U − EgGGA

(6)

and Eg representing the bandgap. When the defect states within the bandgap are occupied, the corrections to the position of transition levels contribute to corrections for defect formation energies. The resulting changes in formation energies are given by Egexpt − EgGGA+U n∆, ∆E = GGA+U Eg − EgGGA f

(7)

where n is the occupancy of defect levels.

Migration Barrier Energies We used climbing image nudged elastic band method 41 as implemented in VASP with 5 images to calculate the barrier energies for the elementary vacancy migration. The images were relaxed until the forces on individual atoms were smaller than 0.05 eV/˚ A. For migration barrier calculations, we used Gamma centered 2x2x1 Monkhorst-Pack k -point grid.

8 ACS Paragon Plus Environment

Page 9 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

To validate the migration pathways, the phonon frequencies of transition states and the vacancies were computed using Phonopy, 53 a software for phonon calculations. Before computing the phonon spectra, the structures were further locally optimized with a forcebased quasi-Newton algorithm until the individual forces on each ion become smaller than 0.001 eV/˚ A. For phonon calculations, 714 displacement supercells of the same size as the defect supercells were generated for each defect and transition state to compute the force constants and the dynamical matrix. The total energy of each of the phonon supercells was computed within a tolerance of 10−8 eV/˚ A.

RESULTS AND DISCUSSION Bulk Cr2 O3 We first computed the bulk properties of Cr2 O3 to investigate the performance of the DFT+U method with PBEsol functional and validate the choice of the Ueff parameter. In this work we used Ueff = 3.7 eV within Dudarev approach consistent with the recommended value in the Materials Project database. 49 In previous DFT+U studies a range of Ueff values between 3.3 and 5 have been used. In particular, Lebreau et al. 27 used Ueff = 5 eV within the same approach, 27 and Rohrbach et al 54 used U = 5 eV and J = 1 eV within GGA+U correction as proposed by Liechtenstein and Dudarev. 35 Using an ab initio method for calculating U and J for Crm On clusters, Mosey et al. 55 obtained U = 3.3 eV and J = 0.1 eV as converged values with increasing cluster size. Shi et al. 56 investigated the Ueff dependence of magnetic properties of Cr2 O3 using LSDA+U and chose U=4.0 eV and J=0.58 eV, or Ueff = U-J = 3.42 as the optimal values. Our results show good agreement for the calculated geometric, electronic and magnetic properties of Cr2 O3 with previous DFT+U simulations and experimental data (Table 2). While these properties appear not to be very sensitive to the choice of Ueff parameter, the major differences between various Ueff values are manifested in the details of the electronic density of states (DOS) for bulk Cr2 O3 . In our simulations O-p 9 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

and Cr-d orbitals are hybridized in both valence and conduction bands (Figure 1). Orbital hybridization at the edge of the valence band indicates that Cr2 O3 is an intermediate charge transfer Mott Hubbard insulator. This prediction is consistent with previous GGA+U studies of Cr2 O3 . 27,54 In addition there is a small gap in the valence band 1.1 eV below VBM. This gap was found in previous Cr2 O3 studies with screened exchange hybrid functional, although at 1.5 eV below VBM. 57 In contrast, this gap is absent or barely present in some of the previously reported GGA+U studies, 27,54 indicating a higher accuracy of our approach in treating the electronic structure of Cr2 O3 within the GGA+U formalism.

Vacancies The effect of vacancy charge on the relaxation of the nearby ions was studied using full cell relaxation under the constant volume conditions. The relaxed positions and the relaxation vectors of the ions surrounding the vacancies are given in Figures 2 and 3. Due to the anisotropy of Cr2 O3 , the relaxations are also anisotropic as shown in the figures. The magnitude of the relaxation is presented in Tables 3 and 4. For Cr vacancies, the O ions nearest to the Cr vacancy undergo a 5.0% to 5.9% outward relaxation (Table 3). No clear correlation between the magnitude of relaxation and vacancy charge state was found in the first coordination shell. The outward relaxation of the next nearest O ions, however, increases with the vacancy charges from 5.7 to 7.0%. The nearest Cr ions undergo an inward relaxation of -5.2 to -6.5%. The inward relaxation decreases for the second nearest Cr ions with relaxation ranging from -3.4 to -4.5%. In the case of O vacancies, the amount of relaxation of the nearby ions is an increasing function of the defect charge (Table 4). The nearest Cr ion relaxes inward for a neutral vacancy, but relaxes outward for charged vacancies, with the amount of outward relaxation being significantly larger for a +2 vacancy charge. A similar trend is also observed for the next nearest Cr ions. The O ions in the second coordination shell undergo a smaller inward relaxation for all charge states ranging from -0.6 to -2.7%. The electronic DOS plots of the optimized vacancy structures are presented in Figures 4 10 ACS Paragon Plus Environment

Page 10 of 51

Page 11 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

and 5. Simulations reveal that the defect energy states corresponding to both Cr and O vacancies are located within the bandgap. Both Cr and O vacancies are of hybrid O-p and Cr-d character, which is consistent with the hybrid nature of the valence and conduction 0 bands. For Cr neutral vacancies, VCr , and for O vacancies with +2 charge, VO2 , the defect

energy levels are unoccupied. The defect levels of Cr vacancy are located in the lower part of the bandgap. The O vacancy introduces one degenerate level in the lower and another in the upper regions of the bandgap. As the Fermi level is increased from VBM to CBM, the defect levels become occupied and the charge of the defect consequently increases. In addition, the spin-splitting of the defect levels is observed for O vacancies. For VO0 also referred to as F-center, the VO is occupied with two electrons and spin splitting occurs even in the valence −3 and conduction bands. For VCr , the defect levels in the bandgap are completely filled. We

note that in previously reported DFT+U simulations by Lebreau et al., 27 O vacancy levels are within the valence band, and only one Cr defect level is within the bandgap. 27 These differences in the position of the vacancy levels are likely to be associated with a larger U value used in that work, which is manifested in the shift in the relative position of O and Cr bands. Vacancy formation energies computed by ab initio thermodynamics characterizes the stability of defects taking into account conditions in the external environment. The influence of the external environment is accounted for by the chemical potentials of the metal and oxygen in metal oxides and can be calculated using equation 1. The range of chemical potentials is determined by the stable phases formed by the constituent elements of the oxide, oxygen partial pressure (pO2 ) and temperature. In this work, we did not consider the effects of pO2 and temperature, but focused on the 0 K DFT computed phase stability. The limiting stable phases for chromia Cr2 O3 are metallic Cr and CrO2 and the corresponding formation enthalpies provide the limits for the chemical potential values of Cr and O. The upper limit of Cr chemical potential is imposed by the formation energy of bulk Cr. Furthermore, the upper limit of Cr chemical chemical potential imposes a lower limit on the O chemical potential

11 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 51

through the formation enthalpy of Cr2 O3 0 µmax Cr = µCr

 µmin = 1/3 ∆HCr2 O3 − 2µ0Cr . O

(8)

In the above equations, µ0Cr is the chemical potential of bulk Cr, and ∆HCr2 O3 is the formation enthalpy of Cr2 O3 . The upper limit on µO is determined by the equilibrium between CrO2 and Cr2 O3 . At the equilibrium between the two phases, the chemical potentials satisfy the conditions

max µmin = ∆HCrO2 Cr + 2µO max 2µmin = ∆HCr2 O3 . Cr + 3µO

(9)

To determine the limits on the chemical potentials, we computed the Cr-O phase diagram using optimized bulk Cr and O2 molecules using PBEsol/DFT and optimized bulk Cr2 O3 along with CrO2 using PBEsol/DFT+U. The optimized bcc AFM Cr has a lattice constant of 2.81 ˚ A, which is in good agreement with the experimental value of 2.91 ˚ A. The optimized O2 has a bond length of 1.22 ˚ A, which compares well with the experimental O2 bond length of 1.21 ˚ A. In contrast, the accuracy of the enthalpies of formation for various Cr oxide phases computed using DFT+U is insufficient due to two factors. The first source of error is the over binding of O2 molecules by semi-local exchange-correlation functionals. A possible solution is to use experimental formation enthalpy for O2 (µ0O2 ) instead of those computed in DFT when generating the phase diagrams of oxides. 58 The corresponding correction to PBESol/DFT formation enthalpies is equal to +0.45 eV for µO 59 in equations 8 and 9. Another source of error is the requirement for using different flavors of DFT to capture the electronic structure of Cr metal and Cr oxides. 60 In particular when computing ∆HCr2 O3 in equation 8, GGA total energies of Cr and O2 are subtracted from GGA+U total energy of Cr2 O3 . To rectify this inconsistency a correction of -4.95 eV was added to the total energy of Cr2 O3 to make 12 ACS Paragon Plus Environment

Page 13 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

the DFT+U computed ∆HCr2 O3 consistent with the experimental enthalpy of formation. The above corrections are equally applicable to stoichiometric oxides and those with charged defects. However when the formation energies of charged vacancies in bandgap materials are computed with DFT using finite size supercells and periodic boundary conditions, additional corrections for the electrostatic interactions between supercells and bandgap errors have to be applied. 61 We applied these corrections as outlined in the Methods section in a sequential manner; first correcting the Coulomb interaction error and then using the resultant data to correct for bandgap related error (see Figure S1 in the Supplementary Materials for the data on the convergence of the alignment potential in the region far from the defect). The electrostatic errors associated with each defect as predicted by the method developed by Kumagai and Oba 38 are given in Table 5. The data show that the magnitude of finite size electrostatic error is significant for charged vacancies and the magnitude of correction increases with vacancy charge. Interestingly, the electrostatic error associated with GGA+U is higher then that in the GGA. The effect is likely due to the weaker screening of the vacancy charges by the core electrons in GGA+U, which reduces self-interaction error in GGA+U, but leads to stronger electrostatic interactions between the defect charges and hence to higher electrostatic error. The transition levels computed with GGA and GGA+U, both of which include electrostatic corrections and the bandgap error corrected levels, are given in Table 6. The data shows that GGA predicted 2/1 and 1/0 VO transition levels are in resonance with the valence and conduction bands, respectively, while all VCr transition levels are in the lower region of the bandgap. In contrast, GGA+U predicts that the transition levels for O vacancies are within the bandgap and the transition levels for Cr vacancies are shifted up with respect to the VBM compared to the GGA results. Using the GGA and GGA+U transition levels, we computed the corrections to GGA+U transition levels according to the approach proposed by Janotti and Van de Walle. 39 While the amount of correction is nearly equal for all

13 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

VCr transition states, the correction varies significantly among the VO transition states. For O vacancies, the correction pushes both the 2/1 and 1/0 transition levels deeper into the bandgap by increasing the energy level of 2/1 transition and by decreasing the energy level of 1/0 transition. For Cr vacancies on the other hand, the correction increases the energy levels of all possible transitions. Using the corrections associated with 2/0 and 2/1 VO and 0/-1, 0/-2, and 0/-3 VCr transition levels, the vacancy formation energies are corrected to account for the underprediction of bandgap by GGA+U according to Eq. 7. The resulting corrections for each of the charge states of Cr and O vacancies are shown in Table 7. The magnitude of this correction increases with the occupancy of the defect energy levels by electrons. Here, VO0 has an occupancy of 2 and VO2+ has 0. For Cr vacancies, the occupancy is proportional to the magnitude of defect charge. The correction significantly reduces the formation energy of −2 −3 VO0 and increases the formation energies of VCr , and VCr . In the case of Cr vacancies, the

correction is zero for the neutral vacancies, whereas the electrostatic and bandgap corrections −2 −3 are large for VCr and VCr . For O vacancies, the electrostatic correction is equal to zero for

neutral vacancy, but the bandgap-related correction is substantial and equal to -0.88 eV. Based on thermodynamics arguments, it can be expected that O vacancies or Cr interstitials are favorable in O-deficient conditions. And in O-rich conditions, either Cr vacancies or O interstitials can be expected to be the dominant defects. The vacancy formation energies in Cr and O rich conditions are plotted with respect to the position of the Fermi level within Cr2 O3 bandgap (Figures 6 and 7). The results reveal that for Cr-rich conditions and at higher Fermi levels, Cr vacancies become energetically more favorable compared to O vacancies, while for O-rich conditions Cr vacancies are always energetically more favorable. The figures also reveal that VO2+ , which has unoccupied defect levels, has the lowest 0 formation energy among the charge states of O vacancies near VBM. Similarly VCr , whose

energy levels are also unoccupied, has the lowest formation energy among all of the charge states of the Cr vacancies near VBM. The plots in Figures 6 and 7 show that the corrections

14 ACS Paragon Plus Environment

Page 14 of 51

Page 15 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

reduce the energy range over which the intermediate charge state (+1) are energetically favorable for VO while increasing that of the intermediate -2 charge state for VCr . Near CBM, VO0 , with two occupied defect levels has the lowest formation energy among all charge −3 states of O vacancies and VCr with completely filled defect energy levels has the lowest

formation energy among all possible charge states of the Cr vacancies. Our calculations also reveal that -1 charge state is not favored for VCr . The error-corrected vacancy formation energies are given in Table 8 for three cases corresponding to the position of the Fermi level near the VBM, near the CBM and in the middle of the bandgap. We also included the formation energies of Cr and O vacancies reported 0 is comparable with previously reported in the literature for comparison. The Ef for VCr

DFT results by Lebreau et al. 27 However, Ef for VO0 is higher by 0.89 eV arising from the 0 , the amount of correction is corrections to defect formation energy discussed above. For VCr

zero and a direct comparison between our results and previous simulations can be made. The 0 calculated VCr formation energy is 4.55 eV, which is in-line with the previously reported 4.84

eV (Table 8). The 0.29 eV difference in formation energy is likely to be due to the lower Uef f value used in our simulations, which affects the position of defect levels with respect to the valence and conduction bands, defect localization and the amount of structural relaxation around the vacancy.

Vacancy based defect complexes Frenkel Defects A Frenkel defect is formed when an ion/atom is displaced from its sublattice and occupies a neighboring interstitial site. In Cr2 O3 , only two thirds of the octahedral sites inside the oxygen hexagonal sublattice are occupied by the Cr ions and the remaining one third can accommodate interstitials. For Cr Frenkel defects, two configurations of Frenkel defects as proposed by Lebreau et al. 27 were evaluated. Of the two configurations, the one where the vacancy and the interstitial are in the neighboring Cr bi-layers has the lowest formation 15 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

energy of 2.6 eV. The other Frenkel defect, where the vacancy and the interstitial are in the same Cr bi-layer, has significantly higher formation energy of 5.18 eV. These formation energies are in good agreement with those reported by Lebreau et al. 27 The electronic DOS plot of the lower energy Frenkel defect configuration as shown in Fig. 8 reveals that the Cr Frenkel defect creates multiple occupied energy states in the bandgap near the valence band. These defect levels consist of hybridized Cr-3d and O-2p states. The Cr Frenkel defect also creates an unoccupied defect level near the conduction band and is composed mainly of Cr-3d states.

Schottky Defects A Schottky defect is a defect complex with two Cr and three O vacancies. Again two configurations are considered: A) two Cr vacancies in neighboring Cr bi-layers with three O vacancies in the intermediate O basal plane, and B) two Cr vacancies in the same Cr bi-layer and three O vacancies in the adjacent O basal planes. Of these, configuration A has a formation energy of 3.98 eV and configuration B 4.77 eV. Electronic DOS pertaining to configuration A given in Figure 9 shows that Schottky defect generates three spin-degenerate levels in the bandgap. One of them is near the VBM and is occupied and the remaining two are deep levels, in the middle of the bandgap. The occupied defect level consists of hybridized Cr-3d and O-2p states and the unoccupied defect levels are primarily composed of Cr-3d states.

Diffusion We have studied the vacancy-mediated diffusion of Cr and O in different charge states. The diffusion pathways considered in this study include vacancy migration to one of the nearest atoms of the same species. The nearest atoms to the Cr and O vacancies are labeled in Figure 10. These pathways are either in the basal plane or along the c-axis and are equivalent to the ones considered by Lebreau et al. 27 16 ACS Paragon Plus Environment

Page 16 of 51

Page 17 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Cr diffusion Figure 10a illustrates the Cr diffusion pathways studied. In this figure the magenta sphere within the octahedron represents Cr vacancy (denoted as V for convenience) and the nearest ions to the Cr vacancy considered for Cr diffusion are labeled. Of these, Cr site “A” is in the same Cr-bilayer as that of the vacancy and is at a distance of 2.91 ˚ A. There are three equivalent sites of this type. “B” denotes the closest Cr sublattice site along the caxis 2.67 ˚ A away from the defect, and “C” is the next nearest Cr ion along the c-axis at a distance of 4.16 ˚ A . It is separated from the vacancy by two O layers and one Cr bilayer. The V-C path has a vacant octahedral site (denoted as “T.D.” in Figure 10a) in the intermediate Cr bi-layer. Simulations predict that for neutral and charged vacancies the basal plane oriented diffusion pathway has lower barrier energies for all defect charge states indicating that low temperature vacancy-mediated Cr diffusion is more probable along the basal plane (Table 9 and Figure 11). Among the pathways along the c-axis, the barrier energies are the lowest for −2 V-C pathway for all vacancy charge states, except for VCr . Among the three pathways con−3 sidered, the barrier energies are lower for VCr mediated diffusion process. For basal plane

oriented diffusion, the barrier energies are inversely proportional to the magnitude of the defect charge. In contrast, no correlation is found between the magnitude of defect charge −1 and the barrier energies for the diffusion along the c-axis via both pathways and VCr has

the largest barrier energy for the diffusion to the nearest Cr bi-layer along the c-axis. The diffusion process along the c-axis via the empty octahedral site shows some interesting dependence on the charge state of the vacancy. The diffusion plot in Figure 11c shows that the intermediate image, where the diffusing Cr ion is in the octahedral site, is stable with respect to the neighboring images. We call the configuration of the Cr occupied (otherwise vacant) octahedral site plus the two vacant Cr sublattice sites in the neighboring Cr bi-layers a triple-defect, the terminology borrowed from intermetallics literature. 62 To calculate the diffusion along this pathway, we first optimized the supercell of the triple-defect 17 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

configuration using the same convergence criteria as those used for the vacancy calculations. Two separate diffusion calculations were then performed with the vacancies and optimized triple-defect configuration as end points and combined into the diffusion plot. The transition states computed with CI-NEB were validated by vibrational analysis. −1 Phonon DOS plots given in Figure 12 show that all Cr vacancies except for VCr have no −1 imaginary frequencies. The phonon DOS for VCr given in Figure S2 of Supplementary −1 Information shows that VCr exhibits a soft mode confirming the conclusion derived from −1 Figures 6 and 7 that VCr is not a stable defect. The transition states of all Cr migration

pathways have one imaginary frequency (shown as negative frequency in the plot), confirming that they indeed are the transition states. The energy of the neutral triple-defect along the c-axis is higher than that of Cr vacancy. However when charged, the triple-defect configuration has the lower energy indicating it is a more stable configuration compared to the similarly charged Cr vacancy configuration. The relative stability of the Cr triple-defect configuration with respect to the Cr vacancy configuration increases with the defect charge. To understand the reason for higher stability of the triple defect configuration relative to the similarly charged vacancy, we calculated the DOS of the triple-defect configuration (Figure 13). The DOS plot reveals that some of defect levels of the triple-defect configuration have opposite spins when compared to those of the vacancies. As the defect charge increases, the spin-flipped defect levels within the bandgap are occupied and in the process the defect energy levels become near degenerate leading to the overall lower energy. To understand the origin of the spin-flipped defect levels of Cr triple defect, we studied the band decomposed charge density corresponding to the defect levels (Figure 14). The plots reveal that the charge corresponding to each of the near-degenerate defect levels localizes on one of the two vacancies and on the neighboring O ions. Furthermore, the spin of the charge localized at one vacancy is opposite to the spin of the charge at the other vacancy. This spin reversal of the charges at the two Cr vacancies is consistent with the AFM ordering of the Cr at the corresponding Cr sites in bulk Cr2 O3 .

18 ACS Paragon Plus Environment

Page 18 of 51

Page 19 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

A similar defect configuration was discovered by Hine et al. 63 to be more stable than an Al vacancy in Al2 O3 for all charge states. This defect was called an Al split-vacancy and considered it to be a relaxed structure of Al vacancy because of the low barrier energy between the Al triple-defect and the Al vacancy. The latter is not the case in Cr2 O3 . The barrier energy between Cr triple defect and Cr-vacancy is significant, approximately equal to 2.5 eV indicating that Cr triple defect is a distinct type. A similar Cr triple defect can be formed along the basal plane by creating in-plane Cr vacancies at the neighboring Cr sites to the interstitial. The possible configurations of a Cr triple defect are shown in Figure 15. We investigated the feasibility of such Cr triple defects by evaluating the energies of triple defects formed by combinations of vertices on the bipyramid that subtend an angle > 120◦ , and found that the formation energies of such triple defect configurations are much greater than those of the triple defect along c-axis (denoted by (A, A’) in Figure 15). The energy difference between a charged Cr triple defect and a similarly charged Cr vacancy is of the order of 0.1 - 0.5 eV suggesting that a significant number of triple defects is present in chromia relative to Cr vacancies despite fewer lattice sites available for Cr triple defects (half of the Cr sublattice sites available for vacancies). The activation energy for Cr diffusion along the V-C pathway, which is the sum of defect formation and migration barrier energies, is equal for both vacancies and triple defects. This indicates that the jump frequency for vacancies and for triple defects along V-C pathway are also equal. Given that the geometric and correlation factors in the diffusion prefactor are different, the corresponding diffusion coefficients for vacancy and triple defect mediated diffusion along the V-C pathway are expected to be different. At higher temperatures (> 1000◦ K), the differences in the migration barrier energies for the basal plane and the c-axis pathways are relatively small (∼ kT). As a result, Cr ion diffusion along the c-axis mediated by triple defects and vacancies appears to be competitive with Cr ion diffusion along the basal plane mediated by Cr vacancies.

19 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

O diffusion O diffusion pathways considered here are shown in Figure 10b. In this Figure the orange sphere within the tetrahedron represents the O vacancy and the nearest O ions to the O vacancy considered for O diffusion are labeled. Of these, atoms “A” and “B” are in the same basal plane as the vacancy at distances 2.63 ˚ A and 3.02 ˚ A, respectively, while atoms “C” and “D” are the closest O sublattice sites in the neighboring basal plane separated by one Cr bi-layer and are at distances of 2.75 ˚ A and 2.87 ˚ A, respectively. The calculated migration barrier energies are given in Table 10 and the barrier energy plots are shown in Figure 16. The results suggest that the V-A pathway has the lowest barrier energies among the basal plane pathways for all vacancy charges. Pathway V-C is characterized by the lowest barrier energies among the c-axis pathways. Among all the pathways considered, V-A has the lowest migration barriers, indicating that O vacancy is more mobile along the basal plane. The barrier energy decreases with increasing defect charge for all pathways. These results are consistent with those reported for α-Al2 O3 . 34 The computed barrier energies show that both Cr and O vacancies are more mobile along the basal plane. Of the Cr and O vacancies, O vacancies have overall lower barrier energies suggesting that O vacancies are more mobile. The phonon DOS plots in Fig. 17 show that the transition states of all the O migration barriers given in Figure 16 have one imaginary frequency, whereas the corresponding vacancies have no imaginary frequencies.

Experimental Comparison We can estimate the possible range of diffusion coefficients for vacancy-mediated diffusion using the vacancy formation energies at phase boundaries, the diffusion barrier energies presented in this work, and the diffusion prefactors computed by Lebreau et al. 27 For example, the upper limit on the vacancy-mediated O diffusion coefficient can be obtained at Cr-rich conditions and assuming that the Fermi level is fixed near VBM. In these conditions, VO2 are 20 ACS Paragon Plus Environment

Page 20 of 51

Page 21 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

the dominant defects with a formation energy of 0.87 eV. Considering that the lowest barrier energy is equal to 1.22 eV (V-A pathway), the upper limit for O self diffusion coefficient (DO ) is obtained as 7.8 × 10−13 cm2 s−1 at 1173◦ K. The lower limit for DO is obtained at O-rich conditions with the Fermi level in the vicinity of CBM. In this case, the dominant O defects are VO0 . The resulting DO is equal to 3.3 × 10−41 cm2 s−1 at 1173◦ K. The experimental DO values 19 range between 1 × 10−15 and 3 × 10−19 and are well within the theoretical limits. The experimental diffusion coefficients are statistical quantities averaged over many defects with various charge states diffusing via multiple pathways. The Fermi level and the chemical potentials are influenced by the impurities, pressure, temperature, and the details of sample preparation. Hence, some of the above mentioned physical parameters need to be accounted for a direct comparison with experimental data. In principle, one can obtain the Cr and O chemical potentials at different temperatures and O2 partial pressures using NISTJANAF thermodynamic tables 64 and a thermodynamic formalism such as the constrained grand-canonical formalism employed in previous works. 33,65,66 Then the Fermi level can be obtained as a function of pO2 and T after applying charge neutrality constraint. For direct comparison with experimental diffusion coefficients, we evaluated the vacancy mediated O diffusion coefficients at 1573◦ K and 3.8 × 10−13 atm ignoring the contributions of electronic entropy and vibrational enthalpy to the defect concentrations. At these conditions, experimental DO was reported 19 to be 6.7 × 10−17 cm2 s−1 , while our predicted DO is equal to 8 × 10−17 cm2 s−1 . (A detailed analysis on temperature dependent diffusion in chromia will be reported in the future.) The predicted self-diffusion coefficient for O refers to the upper limit for the diffusion coefficients pertaining to different defect charges and the dominant defect for O diffusion was found to be VO0 . The predicted DV +1 is nearly equal to DVO0 at O

−17

6 × 10

2 −1

cm s . A close match between the experimental DO and computed DVO indicates

that O vacancies are responsible for O diffusion in Cr2 O3 .

21 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

SUMMARY To gain a fundamental understanding of the diffusion processes in Cr2 O3 , we computed the electronic, thermodynamic, and diffusion properties of vacancies in Cr2 O3 as a first step. The formation energies and the migration barrier energies were evaluated for Cr and O vacancies in various charge states. Using the thermodynamic formalism proposed by Zhang and Northrup 36 we considered Cr-rich and O-rich conditions. The vacancy formation energies were corrected for spurious electrostatic interaction errors arising in the periodic supercell method and for the errors due to the under prediction of the bandgap by GGA+U. The errorcorrected defect transition levels indicate that charged and neutral Cr and O vacancies are deep defects and induce considerable local structural relaxation. Our calculations also predict the existence of a stable triple-defect configuration, where Cr ion occupies an interstitial octahedral site in Cr bi-layer in combination with vacant Cr sites along the c-axis in the neighboring Cr bi-layers. The computed migration barrier energies reveal that O vacancies are more mobile than Cr vacancies. The preferred path for vacancy-mediated diffusion is along the basal plane for all vacancies studied with VO+2 having the highest mobility. The present study is a step towards a comprehensive picture of diffusion in chromia. Additional data on defects in Cr2 O3 such as self interstitials is needed to satisfactorily estimate various experimentally measurable quantities such as the self-diffusion coefficients. We would like to note that formation energies and hence the diffusion coefficients are environmentdependent properties. The results presented in this study along with data on other defects when combined with environmental variables (such as pO2 and T through thermodynamics) can be used to evaluate self-diffusion coefficients for a wide range of conditions thereby enabling a direct comparison with experimental results. Work along these directions is in progress and will be reported in the future.

22 ACS Paragon Plus Environment

Page 22 of 51

Page 23 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Acknowledgement This work was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division. Simulations were performed using PNNL Institutional Computing facility. PNNL is a multiprogram National Laboratory operated by Battelle for the U.S. Department of Energy. B.M. thanks Danny Broberg for his assistance in generating Fig. S1 in Supplementary Materials.

Supporting Information Available Figures showing the potential alignment sampling and the phonon DOS of unstable Cr vacancy.

This material is available free of charge via the Internet at http://pubs.acs.

org/.

References (1) ASM Handbook Series, Vol. 13B Corrosion: Materials, and Vol. 13C Corrosion: Environment and Industries; Eds. Cramer, S. D., B. S. Covino, J.; ASM International: Materials Park, OH, USA, 2005/2006. (2) Van Rooyen, G. T. The Potential of Chromium as an Alloying Element. INFACON 6, Proceedings of the 1st International Chromium Steel and Alloys Congress, Cape Town, South Africa, SAIMM, 1992; pp 43–47. (3) Mu, S.; Wysocki, A. L.; Belashchenko, K. D. Effect of Substitutional Doping on the N´eel Temperature of Cr2 O3 . Phys. Rev. B 2013, 87, 054435. (4) Pan, J.; Waghmare, U. V.; Kumar, N.; Ehi-Eromosele, C. O.; Rao, C. N. R. Effect of Nitrogen and Fluorine Co-substitution on the Structure and Magnetic Properties of Cr2 O3 . ChemPhysChem 2015, 16, 1502–1508. 23 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(5) Pati, S. P.; Shimomura, N.; Nozaki, T.; Shibata, T.; Sahashi, M. Nel Temperature of Cr2O3 in Cr2O3/Co Exchange-Coupled System: Effect of Buffer Layer. J. Appl. Phys. 2015, 117, 17D137. (6) Schreiber, D.; Olszta, M.; Saxey, D.; Kruska, K.; Moore, K.; Lozano-Perez, S.; Bruemmer, S. Examinations of Oxidation and Sulfidation of Grain Boundaries in Alloy 600 Exposed to Simulated Pressurized Water Reactor Primary Water. Microsc. Microanal. 2013, 19, 676–687. (7) Schreiber, D.; Olszta, M.; Bruemmer, S. Grain Boundary Depletion and Migration During Selective Oxidation of Cr in a Ni-5Cr Binary Alloy Exposed to High-temperature Hydrogenated Water. Scr. Mater. 2014, 89, 41 – 44. (8) Yu, J.; Rosso, K. M.; Bruemmer, S. M. Charge and Ion Transport in NiO and Aspects of Ni Oxidation from First Principles. J. Phys. Chem. C 2012, 116, 1948–1954. (9) Sushko, M. L.; Alexandrov, V.; Schreiber, D. K.; Rosso, K. M.; Bruemmer, S. M. Multiscale Model of Metal Alloy Oxidation at Grain Boundaries. J. Chem. Phys. 2015, 142, 214114. (10) Alexandrov, V.; Sushko, M. L.; Schreiber, D. K.; Bruemmer, S. M.; Rosso, K. M. Ab Initio Modeling of Bulk and Intragranular Diffusion in Ni Alloys. J. Phys. Chem. Lett. 2015, 6, 1618–1623, PMID: 26263324. (11) Latanision, R.; Staehle, R. Stress Corrosion Cracking of Iron–Nickel–Chromium Alloys.; 1967. (12) King, A.; Johnson, G.; Engelberg, D.; Ludwig, W.; Marrow, J. Observations of Intergranular Stress Corrosion Cracking in a Grain-Mapped Polycrystal. Science 2008, 321, 382–385.

24 ACS Paragon Plus Environment

Page 24 of 51

Page 25 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(13) Huntz, A. M.; Tsai, S. C. Diffusion in Oxide Scales: Application to Cr2 O3 Scales. J. Mater. Sci. Lett. 1994, 13, 821–825. (14) Sabioni, A. C. S.; Huntz, A. M.; Philibert, J.; Lesage, B.; Monty, C. Relation Between the Oxidation Growth Rate of Chromia Scales and Self-Diffusion in Cr2 O3 . J. Mater. Sci. 1992, 27, 4782–4790. (15) Kofstad, P.; Lillerud, K. Chromium Transport Through Cr2 O3 Scales I. On Lattice Diffusion of Chromium. Oxid. Met. 1982, 17, 177–194. (16) Orman, J. A. V.; Crispin, K. L. Diffusion in Oxides. Rev. Mineral. Geochem. 2010, 72, 757–825. (17) Hoshino, K.; Peterson, N. Cation Self-Diffusion in Cr2 O3 . J. Am. Ceram. Soc. 1983, 66, c202–c203. (18) Sabioni, A.; Lesage, B.; Huntz, A.; Pivin, J.; Monty, C. Self-Diffusion in Cr2 O3 I. Chromium Diffusion in Single Crystals. Philos. Mag. A 1992, 66, 333–350. (19) Sabioni, A.; Huntz, A.; Millot, F.; Monty, C. Self-Diffusion in Cr2 O3 II. Oxygen Diffusion in Single Crystals. Philos. Mag. A 1992, 66, 351–360. (20) Sabioni, A.; Huntz, A.; Millot, F.; Monty, C. Self-Diffusion in Cr2 O3 III. Chromium and Oxygen Grain-Boundary Diffusion in Polycrystals. Philos. Mag. A 1992, 66, 361–374. (21) Schmucker, E.; Petitjean, C.; Martinelli, L.; Panteix, P.-J.; Lagha, S. B.; Vilasi, M. Oxidation of Ni-Cr Alloy at Intermediate Oxygen Pressures. I. Diffusion Mechanisms Through the Oxide Layer. Corros. Sci. 2016, 111, 474–485. (22) Tsai, S.; Huntz, A.; Dolin, C. Growth Mechanism of Cr2 O3 Scales: Oxygen and Chromium Diffusion, Oxidation Kinetics and Effect of Yttrium. Mater. Sci. Eng., A 1996, 212, 6–13.

25 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(23) Tuller, H. L.; Bishop, S. R. Point Defects in Oxides: Tailoring Materials Through Defect Engineering. Annu. Rev. Mater. Res. 2011, 41, 369–398. (24) Tsur, Y.; Riess, I. Self-Compensation in Semiconductors. Phys. Rev. B 1999, 60, 8138– 8146. (25) Catlow, C. R. A.; Corish, J.; Hennessy, J.; Mackrodt, W. C. Atomistic Simulation of Defect Structures and Ion Transport in α-Fe2 O3 and α-Cr2 O3 . J. Am. Ceram. Soc. 1988, 71, 42–49. (26) Vaari, J. Molecular Dynamics Simulations of Vacancy Diffusion in Chromium(III) Oxide, Hematite, Magnetite and Chromite. Solid State Ion. 2015, 270, 10 – 17. (27) Lebreau, F.; Islam, M. M.; Diawara, B.; Marcus, P. Structural, Magnetic, Electronic, Defect, and Diffusion Properties of Cr2 O3 : A DFT+U Study. J. Phys. Chem. C 2014, 118, 18133–18145. (28) Lidiard, A. B. The Mott-Littleton method: An Introductory Survey. J. Chem. Soc., Faraday Trans. 2 1989, 85, 341–349. (29) Hohenberg, P.; Kohn, W. Inhomogeneous Electron Gas. Phys. Rev. 1964, 136, B864– B871. (30) Kohn, W.; Sham, L. J. Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. 1965, 140, A1133–A1138. (31) Freysoldt, C.; Grabowski, B.; Hickel, T.; Neugebauer, J.; Kresse, G.; Janotti, A.; Van de Walle, C. G. First-Principles Calculations for Point Defects in Solids. Rev. Mod. Phys. 2014, 86, 253. (32) Mantina, M.; Wang, Y.; Arroyave, R.; Chen, L. Q.; Liu, Z. K.; Wolverton, C. FirstPrinciples Calculation of Self-Diffusion Coefficients. Phys. Rev. Lett. 2008, 100, 215901.

26 ACS Paragon Plus Environment

Page 26 of 51

Page 27 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(33) Youssef, M.; Yildiz, B. Predicting Self-Diffusion in Metal Oxides from First Principles: The case of Oxygen in Tetragonal ZrO2 . Phys. Rev. B 2014, 89, 024105. (34) Lei, Y.; Wang, G. Linking Diffusion Kinetics to Defect Electronic Structure in Metal Oxides: Charge-Dependent Vacancy Diffusion in Alumina. Scr. Mater. 2015, 101, 20 – 23. (35) Liechtenstein, A. I.; Anisimov, V. I.; Zaanen, J. Density-Functional Theory and Strong Interactions: Orbital Ordering in Mott-Hubbard Insulators. Phys. Rev. B 1995, 52, R5467–R5470. (36) Zhang, S. B.; Northrup, J. E. Chemical Potential Dependence of Defect Formation Energies in GaAs: Application to Ga Self-Diffusion. Phys. Rev. Lett. 1991, 67, 2339– 2342. (37) Freysoldt, C.; Neugebauer, J.; Van de Walle, C. G. Fully Ab Initio Finite-Size Corrections for Charged-Defect Supercell Calculations. Phys. Rev. Lett. 2009, 102, 016402. (38) Kumagai, Y.; Oba, F. Electrostatics-Based Finite-Size Corrections for First-Principles Point Defect Calculations. Phys. Rev. B 2014, 89, 195205. (39) Janotti, A.; Van de Walle, C. G. Native Point Defects in ZnO. Phys. Rev. B 2007, 76, 165202. (40) Krukau, A. V.; Vydrov, O. A.; Izmaylov, A. F.; Scuseria, G. E. Influence of the Exchange Screening Parameter on the Performance of Screened Hybrid Functionals. J. Chem. Phys. 2006, 125 . (41) Henkelman, G.; Uberuaga, B. P.; J´onsson, H. A Climbing Image Nudged Elastic Band Method for Finding Saddle Points and Minimum Energy Paths. J. Chem. Phys. 2000, 113, 9901–9904.

27 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 51

(42) Kresse, G.; Hafner, J. Ab initio Molecular Dynamics for Liquid Metals. Phys. Rev. B 1993, 47, 558–561. (43) Kresse, G.; Hafner, J. Ab initio Molecular-Dynamics Simulation of the LiquidMetal˘Amorphous-Semiconductor Transition in Germanium. Phys. Rev. B 1994, 49, 14251–14269. (44) Kresse, G.; Furthm¨ uller, J. Efficient iterative schemes for ab initio Total-Energy Calculations Using a Plane-Wave Basis Set. Phys. Rev. B 1996, 54, 11169–11186. (45) Bl¨ochl, P. E. Projector Augmented-Wave Method. Phys. Rev. B 1994, 50, 17953–17979. (46) Kresse, G.; Joubert, D. From Ultrasoft Pseudopotentials to the Projector AugmentedWave Method. Phys. Rev. B 1999, 59, 1758–1775. (47) Perdew, J. P.; Ruzsinszky, A.; Csonka, G. I.; Vydrov, O. A.; Scuseria, G. E.; Constantin, L. A.; Zhou, X.; Burke, K. Restoring the Density-Gradient Expansion for Exchange in Solids and Surfaces. Phys. Rev. Lett. 2008, 100, 136406. (48) Dudarev, S. L.; Botton, G. A.; Savrasov, S. Y.; Humphreys, C. J.; Sutton, A. P. Electron-Energy-Loss Spectra and the Structural Stability of Nickel Oxide:

An

LSDA+U Study. Phys. Rev. B 1998, 57, 1505–1509. (49) Jain, A.; Ong, S. P.; Hautier, G.; Chen, W.; Richards, W. D.; Dacek, S.; Cholia, S.; Gunter, D.; Skinner, D.; Ceder, G.; Persson, K. A. Commentary: The Materials Project: A Materials Genome Approach to Accelerating Materials Innovation. APL Mater. 2013, 1, 011002. (50) Broberg, D.; Medasani, B.; Zimmerman, N.; Haranczyk, M.; Canning, A.; Asta, M.; Hautier, G. PyCDT: A Python Toolkit for Modeling Point Defects in Semiconductors and Insulators. 2016, arXiv:cond-mat.mtrl-sci/1611.07481. arXiv.org e-Print archive. https://arxiv.org/abs/1611.07481 (accessed Dec 15, 2016). . 28 ACS Paragon Plus Environment

Page 29 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(51) Polo, V.; Kraka, E.; Cremer, D. Electron Correlation and the Self-Interaction Error of Density Functional Theory. Mol. Phys. 2002, 100, 1771–1790. (52) Queisser, H.; Spaeth, J.; Overhof, H. Point Defects in Semiconductors and Insulators: Determination of Atomic and Electronic Structure from Paramagnetic Hyperfine Interactions; Springer Series in Materials Science; Springer Berlin Heidelberg, 2013. (53) Togo, A.; Tanaka, I. First principles phonon calculations in materials science. Scr. Mater. 2015, 108, 1 – 5 . (54) Rohrbach, A.; Hafner, J.; Kresse, G. Ab initio Study of the (0001) Surfaces of Hematite and Chromia: Influence of Strong Electronic Correlations. Phys. Rev. B 2004, 70, 125426. (55) Mosey, N. J.; Liao, P.; Carter, E. A. Rotationally Invariant Ab Initio Evaluation of coulomb and Exchange Parameters for DFT+U Calculations. J. Chem. Phys. 2008, 129 . (56) Shi, S.; Wysocki, A. L.; Belashchenko, K. D. Magnetism of Chromia From FirstPrinciples Calculations. Phys. Rev. B 2009, 79, 104404. (57) Guo, Y.; Clark, S. J.; Robertson, J. Electronic and Magnetic Properties of Ti2 O3 , Cr2 O3 , and Fe2 O3 Calculated by the Screened Exchange Hybrid Density Functional. J. Phys. Condens. Matter 2012, 24, 325504. (58) Wang, L.; Maxisch, T.; Ceder, G. Oxidation Energies of Transition Metal Oxides Within the GGA + U Framework. Phys. Rev. B 2006, 73, 195107. (59) Tahini, H. A.; Tan, X.; Schwingenschlgl, U.; Smith, S. C. Formation and Migration of Oxygen Vacancies in SrCoO3 and Their Effect on Oxygen Evolution Reactions. ACS Catal. 2016, 6, 5565–5570.

29 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(60) Jain, A.; Hautier, G.; Ong, S. P.; Moore, C. J.; Fischer, C. C.; Persson, K. A.; Ceder, G. Formation Enthalpies by Mixing GGA and GGA + U Calculations. Phys. Rev. B 2011, 84, 045115. (61) Lany, S.; Zunger, A. Assessment of Correction Methods for the Band-Gap Problem and for Finite-Size Effects in Supercell Defect Calculations: Case Studies for ZnO and GaAs. Phys. Rev. B 2008, 78, 235104. (62) Krachler, R.; Ipser, H. Triple-Defect Complexes in the B2 Intermetallic Compound NiAl. Phys. Rev. B 2004, 70, 054113. (63) Hine, N. D. M.; Frensch, K.; Foulkes, W. M. C.; Finnis, M. W. Supercell Size Scaling of Density Functional Theory Formation Energies of Charged Defects. Phys. Rev. B 2009, 79, 024112. (64) Chase, M. W. J. NIST-JANAF Thermochemical Tables, 4th Edition; American Institute of Physics: New York, 1998. (65) Crocombette, J.-P.; Torumba, D.; Chartier, A. Charge States of Point Defects in Uranium Oxide Calculated With a Local Hybrid Functional for Correlated Electrons. Phys. Rev. B 2011, 83, 184107. (66) Li, X.; Finnis, M.; He, J.; Behera, R.; Phillpot, S.; Sinnott, S.; Dickey, E. Energetics of Charged Point Defects in Rutile TiO2 by Density Functional Theory. Acta Mater. 2009, 57, 5882 – 5891. (67) Finger, L. W.; Hazen, R. M. Crystal Structure and Isothermal Compression of Fe2 O3 , Cr2 O3 , and V2 O3 to 50 kbars. J. Appl. Phys. 1980, 51, 5362 – 5367. (68) Fang, P. H.; Brower, W. S. Dielectric Constant of Cr2 O3 Crystals. Phys. Rev. 1963, 129, 1561.

30 ACS Paragon Plus Environment

Page 30 of 51

Page 31 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Table 1: Correspondence between the Kroger-Vink (K-V) representation and the representation used in this work for charged vacancies. Defect

Cr vacancy

O vacancy

Charge

K-V Notation

Notation in this study

0

X VCr

0 VCr

-1

VCr

-2 -3

0

−1 VCr

VCr

00

−2 VCr

VCr

000

−3 VCr

0

VOX

VO0

1

VO•

VO1

2

VO••

VO2

Table 2: Lattice parameters a and c, band gap Eg , Cr magnetic moment µ, and dielectric constant  for bulk Cr2 O3 computed with PBEsol functional. Property This study c (˚ A) 13.680 Eg (eV) 2.95 µ (µB /atom) 2.91  (a-axis) 9.77  (c-axis) 11.35

Calculated 27,54 13.850, 13.839 2.8, 2.6 3.1, 3.01

31 ACS Paragon Plus Environment

Experimental 13.566 67 3.4 67 3.8 67 13.3 68 11.9 68

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table 3: Relaxation of the atoms surrounding Cr vacancy. Positive (negative) numbers indicate outward (inward) relaxation. Atom (# of atoms) O (3) O (3) Cr (1) Cr (3)

Unrelaxed Relaxation (in %) ˚ Distance (in A) q=0 q=-1 q=-2 q=-3 1.98 5.90 5.03 5.24 5.53 2.03 5.74 6.74 6.87 7.00 2.68 -6.54 -5.15 -5.80 -6.25 2.91 -3.68 -4.50 -3.91 -3.40

Table 4: Relaxation of the atoms surrounding O vacancy. Positive (negative) numbers indicate outward (inward) relaxation. For q=1, the two atoms in the first three shells have diferent relaxations shown in two middle columns. Atom (# of atoms) Cr(2) Cr(2) O(2) O(2) O(4)

Unrelaxed Distance (in ˚ A) 1.98 2.03 2.63 2.75 2.87

Relaxation (in %) q=0 q=1 q=1* q=2 -1.53 4.10 1.92 9.24 0.16 5.34 3.18 7.27 -0.86 -1.42 -1.84 -2.70 -0.60 -0.90 -1.30 -0.89 -1.23 -1.69

Table 5: Electrostatic corrections to the vacancy formation energies of Cr2 O3 obtained with PyCDT 50 using extended Freysoldt method. 37,38 Defect

VO

VCr

q

Correction (eV) GGA GGA+U 0 0.000 0.000 1 -0.007 0.113 2 0.129 0.534 0 0.000 0.000 -1 0.325 0.348 -2 0.799 0.987 -3 1.436 1.925

32 ACS Paragon Plus Environment

Page 32 of 51

Page 33 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Table 6: Uncorrected and corrected defect transition levels for chromium and oxygen vacancies. The defect transition levels are computed after applying electrostatic corrections, which are given in Tab. 5. The corrections shown in 5th column are computed using the scheme proposed in Ref. 39. Defect VO

VCr

q/q 0 GGA (q/q 0 ) GGA+U (q/q 0 ) corr (q/q 0 ) (q/q 0 ) 1/0 5.954 2.411 -1.041 1.369 2/0 2.907 1.412 -0.440 0.972 2/1 -0.138 0.413 0.161 0.573 -2/-3 0.888 1.746 0.252 1.999 -1/-3 0.612 1.444 0.244 1.688 -1/-2 0.332 1.142 0.238 1.379 0/-3 0.618 1.385 0.225 1.610 0/-2 0.482 1.206 0.213 1.419 0/-1 0.639 1.271 0.186 1.457

Table 7: Bandgap correction to the vacancy formation energies of Cr2 O3 . The corrections are computed using the scheme proposed in Ref. 39. Defect VO

VCr

q Occupancy 0 2 1 1 2 0 0 0 -1 1 -2 2 -3 3

Correction -0.880 0.161 0.000 0.000 0.186 0.425 0.676

Table 8: Defect formation energies (in eV) of Cr2 O3 with Fermi level (EF ) at 0.1 eV above VBM, 0.1 eV below CBM, and in the middle of the bandgap. Defect (EF =)

VO

VCr Cr Frenkel Schottky

q

0 1 2 0 -1 -2 -3 0 0

VBM+0.1 eV Cr rich O rich 2.61 6.01 1.35 4.74 0.87 4.27 4.55 1.93 5.90 3.29 7.18 4.57 9.08 6.46 2.6 2.6 3.97 3.97

This Study CBM-0.1 eV Cr rich O rich 2.61 6.01 4.55 7.94 7.27 10.67 4.55 1.93 2.70 0.08 0.78 -1.83 -0.52 -3.14 2.6 2.6 3.97 3.97

Literature VBM+Eg /2 Cr rich O rich 2.61 6.01 2.95 6.34 4.07 7.47 4.55 1.93 4.30 1.69 3.98 1.37 4.28 1.66 2.6 2.6 3.97 3.97

33 ACS Paragon Plus Environment

Cr rich

O rich 5.12 27

10.07 25 4.84 27

12.27 25 2.49 27

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table 9: Migration barriers for chromium vacancy mediated diffusion. q denotes vacancy charge state. Diffusion Path (q =) V-A V-B V-C

Migration barrier 0 -1 -2 2.66 2.39 2.11 3.43 3.52 2.64 2.81 2.75 2.86

(eV) -3 2.01 2.37 2.28

Table 10: Migration barriers for oxygen vacancy mediated diffusion. q denotes vacancy charge state. Diffusion Path Migration barrier (eV) (q =) 0 1 2 V-A 2.32 1.95 1.22 V-B 3.65 3.09 2.52 V-C 2.47 2.15 1.46 V-D 3.87 3.45 2.68

34 ACS Paragon Plus Environment

Page 34 of 51

Page 35 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 1: Site and element projected density of states for bulk Cr2 O3 . Dashed vertical line marks the position of the Fermi level.

35 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

0 (a) VCr

Page 36 of 51

−1 (b) VCr

−2 (c) VCr

−3 (d) VCr

Figure 2: Relaxed structure surrounding the Cr vacancy site (denoted with magenta sphere) for different vacancy charge states. The black and magenta arrows indicate the displacement of ions away from and towards the Cr vacancy respectively. The arrows are scaled for visibility such that the arrow representing the maximum displacement has a unit length in each of the subfigures.

36 ACS Paragon Plus Environment

Page 37 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(a) VO0

(b) VO+1

(c) VO+2

Figure 3: Relaxed structure surrounding the O vacancy site (denoted with orange sphere) for different vacancy charge states. The black and orange arrows indicate the displacement of ions away from and towards the O vacancy respectively. The arrows are scaled for visibility such that the arrow representing the maximum displacement has a unit length in each of the subfigures.

37 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 38 of 51

0 (a) VCr

−1 (b) VCr

−2 (c) VCr

−3 (d) VCr

Figure 4: Element projected orbital density of states averaged over all sites for Cr vacancies in different charge states. Dashed vertical line marks the position of the Fermi level.

38 ACS Paragon Plus Environment

Page 39 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(a) VO0

(b) VO+1

(c) VO+2

Figure 5: Element projected orbital density of states averaged over all sites for O vacancies. Dashed vertical line marks the position of the Fermi level.

39 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(a) No Correction

(b) With Q correction only

(c) With Q and BG corrections

Figure 6: Vacancy formation energies with electrostatic (Q) and bandgap related (BG) corrections for Cr rich conditions. Green solid line corresponds to Ef (VCr ) and black dashed line refers to Ef (VO ).

(a) No Correction

(b) With Q correction only

(c) With Q and BG corrections

Figure 7: Vacancy formation energies with electrostatic (Q) and bandgap related (BG) corrections for O rich conditions. Green solid line corresponds to Ef (VCr ) and black dashed line refers to Ef (VO ).

40 ACS Paragon Plus Environment

Page 40 of 51

Page 41 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

total

Cr

O d s p 4

5

6

7 8 Energies (eV)

9

10

11

Figure 8: Element projected orbital density of states averaged over all sites for Cr Frenkel defect. Dashed vertical line marks the position of the Fermi level.

41 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 42 of 51

total

Cr

O d s p 4

5

6

7 8 Energies (eV)

9

10

11

Figure 9: Element projected orbital density of states averaged over all sites for a low energy Schottky configuration. The Schottky defect contains two Cr vacancies in neighboring Cr bi-layers aligned along the c-axis and three O vacancies in the intermediate O basal plane. Dashed vertical line marks the position of the Fermi level.

42 ACS Paragon Plus Environment

Page 43 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(a) Cr

(b) O

Figure 10: Migration pathways considered for (a) Cr and (b) O diffusion. The central sphere denoted as V within the polyhedron (pink sphere in (a) and orange sphere in (b)) represents the vacancy site. Atom A is the nearest Cr site in the Cr bilayer basal plane, and atoms B and C are the Cr sites aligned with the vacancy along c-axis in the adjacent Cr bi-layers. In (b) atoms A and B are on the nearest sites on the basal plane, and C and D represent the nearest O sites along the c-axis.

43 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(a) V-A

(b) V-B

(c) V-C

Figure 11: Migration barriers of Cr diffusion along (a) V-A, (b) V-B and (c) V-C pathways, −1 −2 −3 0 as shown in Fig. 10. Red, green, blue, and black lines represent VCr , VCr , VCr , and VCr , respectively.

44 ACS Paragon Plus Environment

Page 44 of 51

Vacancy V-A TS 30 V-B TS V-C TS

VCr0

20

10

0 VCr-2 Density of States

30

20

10

0 VCr-3 30 Density of States

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Density of States

Page 45 of 51

20

10

0

-10

0

10 Frequency

20

30

Figure 12: Phonon density of states (DOS) of the vacancies and the transition states of Cr diffusion. Imaginary frequencies are plotted as negative frequencies. The V-A, V-B, and V-C pathways of Cr diffusion are described in Fig. 10. The phonon DOS of the V-A TS 0 overlaps with that of the V-C TS in the imaginary frequency region for VCr defect. Similalry, phonon DOS of the V-A TS and the V-B TS overlap in the imaginary frequency region for −2 −1 VCr defect. Phonon DOS for VCr defect (not shown here) exhibits imaginary frequncy and −1 the transition states correponding to VCr migration were not analyzed further.

45 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 46 of 51

0 (a) VCr

−1 (b) VCr

−2 (c) VCr

−3 (d) VCr

Figure 13: Site averaged element projected orbital density of states of the intermediate stable site within the Cr diffusion path along V-C for various vacancy charge states. Dashed vertical line marks the position of the Fermi level.

46 ACS Paragon Plus Environment

Page 47 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(b) Band B

(a) Band A

Figure 14: Band decomposed charge distribution of near degenerate defect levels with opposite spins in Cr triple defect with the overall charge -3. The charge associated with each defect level is localized at one of the two Cr vacancies and the neighboring ions. Magenta, cyan and red spheres represent the interstitial Cr ion, Cr ions in regular Cr sublattice, and O ions, respectively.

47 ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 15: Hexagonal bi-pyramid denoting the possible triple defect configurations. Each triple defect can be represented by vacancies at two of the vertices and the central Cr atom (magenta sphere) in the bi-pyramid. The site pairs of (B,B’), (C,C’), and (D,D’) are symmetrically equivalent. In this study, the pairs of the vertices which subtend an angle > 120◦ at the central Cr atom are evaluated.

48 ACS Paragon Plus Environment

Page 48 of 51

Page 49 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(a) V-A

(b) V-B

(c) V-C

(d) V-D

Figure 16: Migration barriers for O diffusion along (a) V-A, (b) V-B, (c) V-C, and (d) V-D pathways, as described in Fig. 10. Red, green, and blue lines represent VO0 , VO+1 , and VO+2 systems respectively.

49 ACS Paragon Plus Environment

Vacancy V-A TS 30 V-B TS V-C TS V-D TS 20

Page 50 of 51

VO0

10

0 VO1 Density of States

30

20

10

0 VO2 30 Density of States

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Density of States

The Journal of Physical Chemistry

20

10

0

-10

0

10 Frequency

20

30

Figure 17: Phonon density of states of the vacancies and the transition states of O diffusion. Imaginary frequencies are plotted as negative frequencies. The V-A, V-B, V-C, and V-D pathways of the O diffusion are described in Fig. 10.

50 ACS Paragon Plus Environment

Page 51 of 51

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Graphical TOC Entry

Vacancy formation energies and migration barriers in Cr2 O3

51 ACS Paragon Plus Environment