UNDERSTANDING OXYGEN VACANCY FORMATION, INTERACTION, TRANSPORT,
AND STRAIN IN SOFC COMPONENTS VIA COMBINED THERMODYNAMICS AND
FIRST PRINCIPLES CALCULATIONS
By
Tridip Das

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
Chemical Engineering-Doctor of Philosophy
2017

ABSTRACT
UNDERSTANDING OXYGEN VACANCY FORMATION, INTERACTION, TRANSPORT,
AND STRAIN IN SOFC COMPONENTS VIA COMBINED THERMODYNAMICS AND
FIRST PRINCIPLES CALCULATIONS
By
Tridip Das
Understanding of the vacancy formation, interaction, increasing its concentration and diffusion,
and controlling its chemical strain will advance the design of mixed ionic and electronic conductor
(MIEC) materials via element doping and strain engineering. This is especially central to improve
the performance of the solid oxide fuel cell (SOFC), an energy conversion device for sustainable
future.
The oxygen vacancy concentration grows exponentially with the temperature at dilute
vacancy concentration but not at higher concentration, or even decreases due to oxygen vacancy
interaction and vacancy ordered phase change. This limits the ionic conductivity. Using density
functional theory (DFT), we provided fundamental understanding on how oxygen vacancy
interaction originates in one of the typical MIEC, La1-xSrxFeO3-δ (LSF). The vacancy interaction
is determined by the interplay of the charge state of multi-valence ion (Fe), aliovalent doping
(La/Sr ratio), the crystal structure, and the oxygen vacancy concentration and/or nonstoichiometry
(δ). It was found excess electrons left due to the formation of a neutral oxygen vacancy get
distributed to Fe directly connected to the vacancy or to the second nearest neighboring Fe, based
on crystal field splitting of Fe 3d orbital in different Fe-O polyhedral coordination. The
progressively larger polaron size and anisotropic shape changes with increasing Sr-content resulted
in increasing oxygen vacancy interactions, as indicated by an increase in the oxygen vacancy
formation energy above a critical δ threshold. This was consistent with experimental results
showing that Sr-rich LSF and highly oxygen deficient compositions are prone to oxygen-vacancy-

ordering-induced phase transformations, while Sr-poor and oxygen-rich LSF compositions are not.
Since oxygen vacancy induced phase transformations, cause a decrease in the mobile oxygen
vacancy site fraction (X), both δ and X were predicted as a function of temperature and oxygen
partial pressure, for multiple LSF compositions and phases using a combined thermodynamics and
DFT approach. A detailed oxygen vacancy migration barrier calculation gave the oxygen ionic
diffusivity and conductivity.
Oxygen vacancy also causes chemical strain, which was treated as a scalar in the literature.
However, in many materials, it should be a tensor, which is anisotropic. We illustrate this effect
on CeO2, in which it explained a puzzling experiment, which shows significant amplification of
measured strain on applied bias in non-stoichiometric Gd doped ceria. The presence of highly
localized 4f valence orbital in Ce causes charge disproportionation on the formation of neutral
oxygen vacancy, producing anisotropic chemical strain in ceria with cubic symmetry.
Understanding of δ and X and anisotropic chemical strain in the lattice has led to the design of
better MIEC via element doping and strain engineering of the lattice.

Copyright by
TRIDIP DAS
2017

Dedicated to the Tax Payers of United States of America
(for funding my research)

v

ACKNOWLEDGEMENTS

I am grateful to my advisor Prof. Yue Qi for introducing me to the world of atomistic modeling.
She taught me density functional theory (DFT) and helped me to learn many more computational
recipes by challenging my abilities. I also want to thank Prof. Jason D. Nicholas for his constant
guidance and motivation. I am grateful to my advisors for their effort and desire to help me
whenever I needed their support, even on the weekends. I want to thank my committee member
Prof. Wei Lai and Prof. Benjamin Levine for their valuable inputs in my thesis work.
I acknowledge the funding from National Science Foundation (Award Number DMR-1410850,
National Science Foundation CAREER Award No. CBET-1254453), and from Department of
Energy (Award Number DE-FE0023315).
I also acknowledge my colleagues, friends, families and high performance computing center
(HPCC) of Michigan State University for their support.

vi

TABLE OF CONTENTS

LIST OF TABLES ........................................................................................................................... ix
LIST OF FIGURES ..........................................................................................................................x
1.

Introduction............................................................................................................................. 1
1.1
Factors Influencing Oxygen Vacancy Formation in Perovskite ..................................... 6
1.2
Oxygen Vacancy Interaction via Vacancy Polaron ...................................................... 11
1.3.
Effect of the oxygen vacancy interaction on Ionic conductivity at SOFC operating
conditions .................................................................................................................................. 13
1.4.
Anisotropic Chemical Strain Effect Due to Oxygen Vacancy Polaron......................... 15

2.

Computational Details: Methodology developed and Used ................................................. 17
2.1.
Magnetic Moment Interpreted Oxidation State Prediction .......................................... 18
2.1.1. Selection of U Parameter .......................................................................................... 19
2.1.2. Determining Oxidation State on Fe .......................................................................... 23
2.1.3. Advantage over Bader Charge.................................................................................. 24
2.2.
Determining Polaron Shape and Size from the Valance Change on Transition Metal 26
2.3.
Calculation Method of Oxygen Vacancy Formation and Migration Barrier Energy .. 27
2.3.1. Oxygen Vacancy Formation Energy at 0 K .............................................................. 27
2.3.2. Predicting Vacancy Ordered Phase Transition ........................................................ 28
2.3.3. Calculating Oxygen Vacancy Migration Barrier ..................................................... 29
2.4.
Thermodynamic Extension of DFT Calculated Energies ............................................. 30
2.4.1. Converting Vacancy Formation Energy to Vacancy Formation Free Energy ......... 30
2.4.2. Connecting Chemical Potential of Oxygen from 0 K to T K..................................... 31
2.5.
Predicting Ionic Conductivity as a Function of Temperature and Phase Change in Air .
....................................................................................................................................... 36
2.5.1. Oxygen Vacancy Concentration in Bulk Lattice ....................................................... 36
2.5.2. Mobility and Diffusivity of Oxygen Vacancy in Bulk Lattice .................................... 37
2.6.
Anisotropic Chemical Strain Due to Oxygen Vacancy ................................................. 37

3. DFT+U as a Reliable Tool to Predict Lattice Structure and Electronic Properties of MIEC
at Low Temperature ...................................................................................................................... 41
3.1
Predicting Ground State Lattice Structure and Magnetic Ordering in La1-xSrxFeO3-δ
Phases ....................................................................................................................................... 42
3.2.
Density of States for LaFeO3-δ and La0.5Sr0.5FeO3-δ ..................................................... 43
3.3.
Study of Different SrFeO3-δ Phases............................................................................... 45
3.3.1 Magnetic Orientations of SrFeO3-δ Phases ............................................................... 47
3.3.2. Density of States for SrFeO3-δ Phases....................................................................... 48
3.3.3. Formation Enthalpy of Brownmillerite SrFeO2.5...................................................... 49
3.3.4. d-Orbital Splitting and Its Impact on the Iron Oxidation State ................................ 49
3.3.5. Solving the Debate on Fe Oxidation State at Square Pyramidal and Octahedral Sites
................................................................................................................................... 50
vii

3.4.
4.

Choices of Crystal Structures at SOFC Operating Temperature ................................. 51

Oxygen Vacancy Interaction via Vacancy Polaron .............................................................. 53
4.1.
Origin of Long Range Charge Transfer due to Oxygen Vacancy Polaron................... 53
4.2.
Oxygen Vacancy Polaron in Different La1-xSrxFeO3-δ Phases...................................... 55
4.2.1. LaFeO3-δ ..................................................................................................................... 56
4.2.2. SrFeO3-δ...................................................................................................................... 58
4.2.3. Rhombohedral and Cubic La0.5Sr0.5FeO3-δ ................................................................ 61

5. Oxygen Vacancy Formation sites and Formation Energy for Interacting and Noninteracting Vacancies at 0 K ......................................................................................................... 63
5.1. Oxygen Vacancy Formation Sites in SrFeO3-δ Phases ..................................................... 63
5.2. Oxygen Vacancy Formation Energies in SrFeO3-δ Phases ............................................... 66
5.3. Computationally Predicted Vacancy Ordered Phase Transition in SrFeO3-δ Phases ...... 68
6.

Oxygen Vacancy Site Fraction as a Function of Temperature and Phase Change (in air) . 72
6.1.
Oxygen Vacancy Site Fractions in Different SrFeO3-δ Phases Predicted with
Temperature in Air .................................................................................................................... 73
6.2.
Oxygen Vacancy Site Fractions in Different LSF Compositions Predicted with
Temperature in Air .................................................................................................................... 76
6.3.
Effect of the Variation of U-Parameter on Oxygen Vacancy Formation ..................... 79

7. Oxygen Ionic Conductivity in La1-xSrxFeO3-δ – Considering Oxygen Vacancy Interaction
and Phase Change ........................................................................................................................ 82
7.1
Oxygen Vacancy Migration in La1-xSrxFeO3-δ .............................................................. 83
7.1.1 Oxygen Vacancy Migration in LaFeO3-δ .................................................................. 84
7.1.2 Oxygen Vacancy Migration in La0.5Sr0.5FeO3-δ ........................................................ 85
7.1.3 Oxygen Vacancy Migration in SrFeO3-δ ................................................................... 88
7.2
Oxygen Ionic Conductivity in La1-xSrxFeO3-δ as a Function of Temperature (in Air) .. 92
8.

Anisotropic Chemical Strain Produced by Non-Stoichiometric Cubic CeO2-δ ..................... 95
8.1
Charge Disproportionation in Non-stoichiometric Ceria ............................................ 95
8.2
PDOS of Perfect and Oxygen Vacancy Adjacent Ce .................................................... 97
8.3
Calculating Anisotropic Chemical Strain ..................................................................... 99

9.

Conclusion and Future Work .............................................................................................. 102
9.1
Thesis Summary and Conclusion ................................................................................ 102
9.2
Future Work ................................................................................................................ 104

APPENDIX.................................................................................................................................. 107
BIBLIOGRAPHY ......................................................................................................................... 110

viii

LIST OF TABLES

Table 2.1. Spin-polarized Bader charge analysis of different LSF phases. Magnetic moment
calculated by VASP and expected formal charge on Fe are also provided for comparison. ....... 25
Table 3.1. DFT predictions of stable LSF structures under ground-state conditions. ................. 41
Table 3.1. Comparison of DFT-calculated and experimentally-measured lattice parameters,
magnetic ordering, and band gaps for several phases of SrFeO3-δ (“See DOS” in the last column
of below-provided table refer to the Figure 3.1 for the density of states plot). N.C. stands for Not
Calculated. Experimental NĂŠel Temperature, TN, is also listed for each phase. .......................... 45
Table 5.1. The vacancy formation energies at different oxygen site calculated for all four
crystalline structures. Each site is noted by its label in Figure 1.1, its position in Fe-O polyhedra,
and the charge of the Fe ions in the Fe-O polyhedra. The bold values indicate the most favorable
oxygen vacancy site....................................................................................................................... 64
Table 8.1. Comparing calculated energy and bond lengths for different Vo-formation cases. .... 95

ix

LIST OF FIGURES

Figure 1.1. Vacancy ordered phases of SrFeO3-δ: (a) cubic, (b) tetragonal, (c) orthorhombic, and
(d) brownmillerite with increasing oxygen non-stoichiometry. The energetically favorable sites for
vacancy creation are highlighted with dashed red circles. Fe-O polyhedra are colored according
to the calculated charge state on Fe, as shown in the legend, with Fe-O octahedra and square
pyramids containing Fe+3.9 or Fe+4 dark brown, Fe-O octahedra containing Fe+3.6 yellowishgreen; Fe+3.3 light purple; Fe+3 light brown, and Fe-O tetrahedra containing Fe+3.2 light orange.
......................................................................................................................................................... 8
Figure 2.1. DFT calculated lattice structures for (a) orthorhombic LaFeO3, (b) rhombohedral
La0.5Sr0.5FeO3, (c) cubic La0.5Sr0.5FeO3 (high temperature phase), and (d) cubic SrFeO3. Atoms
are colored as La – green, Sr - blue, Fe – brown, O – red. The different Fe-O6 octahedra colors
denote polyhedra with different Fe charge states: Fe3+ – light brown, Fe3.5+ – moss green, Fe4+ –
dark brown. ................................................................................................................................... 18
Figure 2.2. The plot of (a) magnetic moment and (b) lattice parameters as a function of Ueff
parameter. The straight lines are experimental data, and blue and red correspond to the data for
SrFeO3 and LaFeO3, respectively. Since LaFeO3 is not cubic, red, green and pink colors in (b)
are corresponding to scaled-lattice (a/√2, b/2, c/√2 respectively) of LaFeO3. ............................ 20
Figure 2.3. Black squares represent the formal charge on Fe resulting from different La/Sr ratios
in LaxSr1–xFeO3 calculated from the magnetic moment on Fe with a GGA+U(=3) (left y-axis). Blue
diamonds show experimentally measured average Fe charge states for varying La/Sr ratios (right
y-axis). The black dotted line represents the predicted linear relationship between the Fe charge
state and the La/Sr ratio. .............................................................................................................. 23
Figure 2.4. Compare different methods to obtain oxygen connection energy in Equation 2.5. .. 35
Figure 2.5. (a) CeO2 unit cell showing O-Ce4 tetrahedral coordination and Ce-O bonding
directions. (b) Charge distribution on Ce around neutral oxygen vacancy formation showing in a
2x2x2 supercell. ............................................................................................................................ 38
Figure 3.1. GGA+U (=3) partial density of states calculations for (a-b) orthorhombic LaFeO3,
(c-d) rhombohedral La0.5Sr0.5FeO3. The gray and black lines represent the Fe 3d and O 2p states
respectively. .................................................................................................................................. 44
Figure 3.2. GGA+U (Ueff = 3) PDOS for four different SrFeO3-δ phases. The gray and black line
represent Fe 3d and O 2p states respectively. The stoichiometric SrFeO 3 has shown metallic
behavior. Whereas, band gaps are observed in oxygen deficient phases. .................................... 46
Figure 3.3. Crystal Field Theory explains d-orbital splitting and electronic distribution for
different Fe-O polyhedra. (a) Square pyramidal Fe-O coordination, (b) octahedral Fe-O
coordination (c) Compressed octahedral Fe-O coordination with Jahn-Teller distortion. ......... 50

x

Figure 4.1. The charge distribution on Fe due to oxygen vacancy formation is described on
vacancy containing lattice plane of a 4x4x4 cubic SrFeO3 supercell with a single oxygen vacancy
at the center of the cell. Figure (a) is the lattice plane without relaxing the atom position (only
electron cloud is allowed to relax after vacancy formation). Figure (b) is the lattice plane
containing atoms after relaxing the atom positions. The light brown atoms are Fe and red atoms
are O. The light blue Sr atoms are above or below the described planes. The oxygen vacancy is
denoted with a black dot. Figure (a) and (b) contains the oxygen vacancy at the center of the (200)
lattice plane. The charge on every Fe atom, rounded to the nearest tenth, is listed next to each Fe
atom............................................................................................................................................... 53
Figure 4.2. Long-range charge transfer areas (denoted by the purple dashed lines in a-d) and
perspective views of the LaxSr1–xFeO3 lattice polarons (e-h) for (a, e) orthorhombic LaFeO3, (b,f)
rhombohedral La0.5Sr0.5FeO3, (c,g) cubic La0.5Sr0.5FeO3, and (d,h) cubic SrFeO3. Note, in (a-d)
different Miller planes were used for each phase to compare lattice orientation with same visual
conformity. In (e-h) the polaron shape was highlighted by only drawing sides to those Fe-O
polyhedra with a Fe atom that received excess electrons due to oxygen vacancy formation. In (ad), the black lines indicate the border of the supercells used for computations, while the arrows
provide an estimate of the polaron size by denoting the distance the distance from the outermost
charge-altered Fe atoms to the oxygen vacancy site. ................................................................... 56
Figure 4.3. Schematic of the d-orbital splitting in the Fe in different polyhedra showing that the
excess electrons generated during oxygen vacancy formation (in red) will be localized to oxygen
vacancy adjacent square pyramidal Fe in orthorhombic LaFeO3 (a), and transferred to second
nearest neighboring octahedral Fe in cubic SrFeO3 (b). ............................................................. 57
Figure 4.4. Oxygen vacancy polaron shape in (a) tetragonal, (b) orthorhombic, (c) brownmillerite
SrFeO3-δ. Fe-O polyhedra colored according to the charge on Fe.27 Here, blue atoms are Sr,
brown are Fe and red are O. ........................................................................................................ 60
Figure 4.5. d-orbital splitting and energy distribution for different Fe-O polyhedra was derived
from crystal field theory. (a) Square planar Fe-O coordination, (b) octahedral FeO coordination
(c) tetrahedral Fe-O coordination. ............................................................................................... 60
Figure 5.1. The vacancy formation energy per oxygen vacancy, calculated from DFT for four
different SrFeO3-δ phases. ............................................................................................................. 66
𝛿

Figure 5.2. Energy required ∆𝐸𝛿 to form the non-stoichiometric structures (𝑆𝑟𝐹𝑒𝑂3−𝛿 + 2 𝑂2 )
from cubic SrFeO3. ....................................................................................................................... 68
Figure 5.3. GGA+U calculated oxygen vacancy formation energies for different LaxSr1–xFeO3
compositions. The Fe-O polyhedra represent the charge on Fe in a perfect lattice. The lines are
guides to the eye. ........................................................................................................................... 69
Figure 6.1. Dotted line represent δ increasing with T at pO2 = 0.21 atm according to Equation 1.6
(dilute vacancy assumption) for cubic SrFeO3-δ (∆𝐸𝑣𝑎𝑐𝑓 from Figure 5.1 for different δ). Solid
line represent δ increasing with T at atmospheric condition according to Equation 2.17
(interacting vacancy assumption) for cubic SrFeO3-δ (where ∆𝐸𝑣𝑎𝑐𝑓 is a function of δ). .......... 72
xi

Figure 6.2. (a) Vacancy formation free energy with nonstoichiometry for all the phases (b)
Nonstoichiometry is plotted as a function of temperature. (c) Vacancy formation site fraction is
plotted as a function of temperature. (the data point corresponding to 3Xcubic = 0.5 is neglected
from fitting as at this vacancy site fraction cubic SrFeO3-δ is easily transformed to brownmillerite
and hence loses its available vacancy sites). ................................................................................ 73
Figure 6.3. (a) The oxygen non-stoichiometry versus temperature in air,(b) the oxygen vacancy
site fraction as a function of temperature in air, and (c) the oxygen vacancy site fraction as a
function of temperature and La/Sr ratio for various LaxSr1–xFeO3 phases. Note, the lines in (a-b)
and are guides to the eye and the lines in (c) are the experimentally determined phase
boundaries.44,65 The oxygen vacancy site fraction in (c) is denoted by the color overlay and the
right-hand color scale. The X values in b) and c) are never zero. They are just smaller than can be
resolved with the linear X scale chosen to highlight the X trends discussed here........................ 77
Figure 6.4. (a) Compare the impact of U parameters on oxygen vacancy formation energy as a
function of oxygen nonstoichiometry; (b) compare the calculated oxygen nonstoichiometry as a
function of temperature with different U parameters with experimental values. ......................... 80
Figure 7.1. Oxygen migration in LaFeO3-δ. The pink atoms in (a) highlight the oxygen migration
path, (b) shows the oxygen migration barrier along this path for two different δ, and (c) shows the
charge change on the neighboring Fe atoms as a function of oxygen atom coordinate on its path
(at δ =0.06). .................................................................................................................................. 85
Figure 7.2. Oxygen migration in LSF55. The pink atoms in (a) highlight the oxygen migration
path in Rhombohedral LSF55 with δ =0.04, and (b) highlight the oxygen migration path in cubic
LSF55 with δ =0.13. Figure (c) shows the oxygen migration barrier in rhombohedral and cubic
LSF55. (d) shows the charge change on the oxygen vacancy neighboring Fe atoms in
rhombohedral LSF55 as a function of oxygen atom coordinate (corresponds to Figure (a)). Figure
(e) shows the charge change on the oxygen vacancy neighboring Fe atoms in cubic LSF55 as a
function of oxygen atom coordinate (corresponds to Figure (b))................................................. 87
Figure 7.3. Oxygen migration in cubic SrFeO3-δ. The pink atoms in Figure (a) highlight the oxygen
migration path, Figure (b) shows the oxygen migration barrier along this path for two different δ,
and Figure (c) shows the charge change on the neighboring Fe atoms as a function of oxygen atom
coordinate on its path. .................................................................................................................. 89
Figure 7.4. Oxygen migration in tetragonal SrFeO3-δ. (Left) non-red smaller atoms highlights the
oxygen migration paths. Paths are defined with Wyckoff positions of oxygen atoms in the perfect
lattice; P1: 4c  16k  8h  16k  4c; P2: 4c  16k  16m  16k  4c; P3: 4c  16k 
4c. (Right) oxygen migration barrier at different coordinates (Îś) of the path. ............................. 90
Figure 7.5. Oxygen migration in orthorhombic SrFeO3-δ. (Left) non-red smaller atoms highlight
the oxygen migration paths. Paths are defined with Wyckoff positions of oxygen atoms in the
perfect lattice; P1: 2b  16r  16r  2b (c-direction); P2: 2b  16r  4j  16r  2b; P3: 2b
 16r  16r  2b (a-direction). (Right) oxygen migration barrier at different coordinates (Μ) of
the path.......................................................................................................................................... 90

xii

Figure 7.6. Oxygen vacancy (a) concentration (b) ion diffusivity, and (c) ion conductivity in
different LSF compositions as a function of temperature in air. In Figure legends ‘cubic SFO’
represent cubic SrFeO3-δ data. ‘other SFO’ represent combined data for tetragonal, orthorhombic,
and brownmillerite phase in Figure 7.6(a) and tetragonal and orthorhombic phases in Figure
7.6(b,c). ‘LSF55(combined)’ represents the LSF55 data with phase change from rhombohedral to
cubic. ‘rhomb LSF55’ in Figure 7.6(c) represent the computed oxygen ionic conductivity data for
rhombohedral LSF55 phase. ......................................................................................................... 91
Figure 8.1. Charge distribution on Ce and displacement of O-atoms around neutral oxygen
vacancy (𝑉𝑂•• ). Red-O moved away from oxygen vacancy site. Green-O came closer to oxygen
vacancy site. Red Ce-O bond signify elongation with respect to bond length in perfect structure
and blue Ce-O bond signify contraction. ...................................................................................... 96
Figure 8.2. Partial density of states plot for Ce-5s, 5p, 5d, 4f and O-2p orbitals in pure CeO2 . 97
Figure 8.3.(a) PDOS for Ce4+-5d and Ce4+-4f orbital for Ce4+ adjacent to 𝑉𝑂•• , after 𝑉𝑂•• formation
in asymmetric case. Orbitals show same shape as Ce4+ in perfect CeO2.(b) PDOS for Ce3+-5d and
Ce3+-4f orbital for Ce3+ adjacent to 𝑉𝑂•• , after 𝑉𝑂•• formation in asymmetric case. Figure shows
localization of excess electron in 4f orbital. .(c) PDOS for Ce3.5+-5d and Ce3.5+-4f orbital for Ce3.5+
adjacent to 𝑉𝑂•• , after 𝑉𝑂•• formation in symmetric case. Figure shows localization of excess
electron in 4f orbital. This scenario is computational artifact as the structure was energetically
unfavorable. .................................................................................................................................. 98
Figure 8.4. The oxygen vacancy formation energy varies with an applied strain in different
directions. The slope of the strain versus energy curve gives the elastic dipole tensor mapped with
applied strain direction. ................................................................................................................ 99
Figure 8.5.(a) Average expansion of single crystal ceria with increasing oxygen nonstoichiometry.
(b) possible orientations of Ce3+ and Ce4+ atoms around oxygen vacancy................................ 101
Figure 9.1. Ruddlesden-Popper type Sr3Fe2O7 lattice with single oxygen vacancy (black atom) at
the center. Blue, brown, and red atoms are Sr, Fe, and O respectively. Sr3Fe2O7 can be a promising
compound for low-temperature SOFC cathode. ......................................................................... 105
Figure 9.2. Crystal structure of 12.5% Gd/Sm/Pr (purple atoms) doped ceria (GDC/SDC/PDC).
Here, Ce and O atoms are colored in yellow and red respectively. Gd/Sm/Pr atoms are randomly
distributed in the lattice. ............................................................................................................. 106

xiii

1.

Introduction

Mixed Ionic Electronic Conductor (MIEC) oxides are used for a variety of electrochemical devices
including gas sensors,1 gas separation membranes,2 memristors,3 solid oxide fuel cell (SOFC)
electrode,4–7 opto-electronics8 etc. Different MIEC oxides with their aforementioned applications
can be summarized by their usability as components of a single device, SOFC. It is a chemical to
electrical energy conversion device suitable for a sustainable future, due to its high energy
conversion efficiency, fuel flexibility and low cost.9 In a SOFC device, oxygen from air gets
incorporated in the surface layer of the cathode as oxygen ion by receiving two electrons from the
cathode (transported through external circuit). Then the oxygen ion transports through cathode,
electrolyte, anode, and to the anode surface, where it reacts with the fuel, producing water and
electricity (if the fuel is hydrogen) or carbon-di-oxide, water, and electricity (if the fuel is natural
gas or syngas). Based on the type of MIEC compound, electronic or oxygen ionic transport prevails
over one another.10,11 For example, the SOFC electrolyte is desired to transport only oxygen ion
through it and it should block the transport of the electron. Therefore, fluorite type compounds:
CeO2-δ,12 Gd1-xCexO2-δ (GDC),13 Pr1-xCexO2-δ (PDC),14 Sm1-xCexO2-δ (SDC),15 or Y1-xZrxO2-δ
(YSZ)16 etcetera with fast oxygen ion conductivity are used as SOFC electrolyte. Whereas, SOFC
electrodes are desired to have good electron and oxygen ion conductivity.4 As a result compounds
of ABO3 perovskite family: La1-xSrxFeO3-δ (LSF),17 La1-xSrxCoO3-δ (LSC),18 La1-xSrxCo1-yFeyO3-δ
(LSCF),7 La1-xSrxMnO3-δ (LSM),19 La1-xSrxCo1-yMnyO3-δ (LSCM),20 Sm1-xSrxCoO3-δ (SSC),21 Ba1xSrxCo1-yFeyO3-δ

(BSCF)2 etcetera serve as the main component for SOFC electrode. Though

conventional cathode was produced in a combination of good electronic conductor, LSM mixed
with GDC, a good ionic conductor22 LSF or LSCF can independently work as a MIEC cathode.23
Unfortunately, the performance of these MIEC compounds is often limited by poor oxygen
ion transport within the material at low temperature (600oC or less).24 To maintain desired oxygen
1

ionic conductivity, it is required to operate the SOFC at high temperature (800-1000oC), which
causes premature aging for different SOFC components.25,26 Further, oxygen ion conductivity does
not exponentially increase with temperature due to oxygen vacancy interaction at higher vacancy
concentration.27–29 Therefore, to improve the performance of these MIEC compounds at low
temperature via doping or lattice strain engineering, we need a thorough understanding of how
different parameters (discussed later) regulate the oxygen ion conductivity in the MIEC. The
oxygen ion conductivity consists of three components, as shown in Equation 1.1.
𝜎 = 𝑐𝑍𝑒𝜇

[1.1]

Where c is the vacancy concentration. Z is the charge on an ion as a multiple of e (the charge of
an electron), and Îź, known as the mobility of the ion. Oxygen vacancy formation and their selfinteraction directly regulate oxygen vacancy concentration and oxygen ion diffusivity. Local
charge redistribution (after oxygen vacancy formation) on (vacancy neighboring) transition metal
and local lattice deformation around the vacancy-site (aka chemical strain field) causes oxygen
vacancy interaction at a higher vacancy concentration, which largely influences the performance
of these MIECs at high temperature. This thesis will investigate the oxygen vacancy formation,
interaction, transport, and its’ chemical strain effect in SOFC components with a focus on two
important members of oxide families: perovskite type, LSF, and fluorite type, CeO2-δ.
Many factors influence the oxygen vacancy formation and migration in perovskite LSF.
These factors include the La/Sr ratio, phase change, charge distribution on multivalent Fe, and the
interactions between the oxygen vacancies in LSF. These factors jointly impact the variation in
oxygen ionic conductivity in LSF at SOFC operating conditions.17,24,30,31 To fully understand the
complicated interplay among oxygen non-stoichiometry, δ, the varying charge on Fe around
oxygen vacancy site, and phase change, several computational challenges have to be solved.17 First,
it was difficult to predict the Fe oxidation state when Fe is present in a lattice with multiple
2

valences. 32–36 With an accurate prediction of the Fe oxidation state in this thesis, the change in Fe
oxidation state can be used to illustrate the excess electron distribution due to oxygen vacancy
formation. The excess electron is also called “polaron”, whose shape and volume were used to
explain the oxygen vacancy interaction. The interacting oxygen vacancy can lead to increasing
vacancy formation energy with vacancy concentration, which means the vacancy cannot be treated
as dilute defects anymore. Therefore, a new computational model is required and was developed
in this thesis, to predict oxygen vacancy concentration dependent nonstoichiometry in both dilute
and non-dilute cases. Using computationally predicted oxygen vacancy concentration and
migration barrier data, the oxygen ionic conductivity can be calculated for LSF as a function of
temperature and phase change, in the air. These resutls will be used to explain the loss of oxygen
conductivity due to vacancy ordered phase transition observed in experiments.37,38
Oxygen vacancy induced chemical strain in the lattice plays a significant role in mechanical
failure at the electrode-electrolyte interface of SOFC. To enhance the mechanical compatibility
between different parts and to arrest the failure at the interface, it is important to understand how
the oxygen vacancy produces this chemical strain. Previous experimental study39 investigated local
lattice distortion and hypothesized anisotropic lattice distortion due to oxygen vacancy in a cubic
CeO2-δ, but couldn’t clearly explain the complicated lattice distortion. Due to the cubic symmetry
of these fluorite structures, earlier computational studies calculated average chemical strain in
place of complete strain tensor.40,41 As a result, they could not capture the anisotropic chemical
strain in the cubic lattice. Therefore, previous experimental and computational studies could not
provide a detailed explanation for the origin of huge chemical strain in GDC thin film on the
application of electric field. This thesis will illustrate the reason behind anisotropic chemical strain
with quantification of it in cubic ceria and explain the complicated lattice distortion around oxygen

3

vacancy. The fundamental reason behind the origin of the chemical strain tensor and its tensorial
relationship with the crystallographic orientation of the lattice will help for better development of
the SOFC electrolyte and phase shift device.42
In summary, this thesis elucidates the role of oxygen vacancy formation, interaction and
transport, which defines the performance of the MIEC compounds (SOFC components) at SOFC
operating conditions, using first principle calculations combined with the thermodynamic
approach. In the rest of Chapter 1, the present challenges and existing literature are discussed for
specific materials and topics in more details. The rest of the thesis will be arranged as the
following:
Chapter 2 introduced newly developed computational methods required for this study used
to solve the defined challenges. For density functional theory (DFT) calculations, Hubbard-U
parameter was selected based on the experimental magnetic moment and lattice length, to correctly
predict the charge on multivalent Fe. Magnetic moment based charge prediction method for
multivalent Fe solved a long-standing debate in literature for charge assignment in a lattice
containing both Fe3+ and Fe4+. This chapter also introduced the method to calculate oxygen
vacancy polaron shape and size from valence change on transition metal. The methods to calculate
oxygen vacancy formation energy and migration energy barrier with varying nonstoichiometry
were also described. A new thermodynamic method developed to calculate oxygen vacancy site
fraction, X for interacting vacancies was introduced. Methods to calculate oxygen vacancy
concentration, ionic diffusivity, mobility and finally ionic conductivity were described. The
derivation for anisotropic chemical strain tensor from energy change due to local deformation and
the long-range elastic strain was also introduced in this chapter.

4

In Chapter 3, the reliability of first-principle calculations as a useful tool to predict ground
state properties of different LSF compositions and phases was established by comparing the
computed lattice structures, electronic properties like magnetic ordering and band gaps with
experiments. This chapter also described all the LSF compositions and phases used in this thesis.
In Chapter 4, oxygen vacancy polaron shape and size were predicted from charge change
on Fe around oxygen vacancy site. The shape and size of this polaron were explained from crystal
field theory with the d-orbital splitting in Fe due to different Fe-O polyhedral coordination. This
predicted polaron shape and size explained the reason behind oxygen vacancy interaction in
different LSF compositions and lattice structures.
Chapter 5 determined the oxygen vacancy formation site and how formation energy
changes as a function of oxygen non-stoichiometry in different LSF phases. The relationship
between the oxygen vacancy formation site and the charge state on Fe was discussed. The critical
vacancy concentration, beyond which oxygen vacancy formation energy increases with
concentration, is correlated with the vacancy polaron size. When oxygen vacancy formation
energy increases with its concentration, the vacancies are considered to be interacting. The
interacting vacancies can lead to phase change, which was also discussed in this chapter.
In Chapter 6, both oxygen vacancy site fraction X and non-stoichiometry δ were predicted
(in air) for different LSF compositions as a function of temperature and phase change. This
calculation utilized the newly developed thermodynamic method for non-dilute (interacting)
vacancies. The agreement of the calculated and experimental values of δ over a broad temperature
range validates our model. Vacancy ordering induced phase change in LSF causes dropping X
with increasing δ, but it is experimentally difficult to separate X and δ. Therefore X was computed,
since it contributes to the mobile oxygen vacancy concentration in LSF.

5

In Chapter 7, oxygen vacancy migration barrier was calculated in different LSF
compositions and phases. The oxygen vacancy migration pathway and how the migration
influences the oxidation state on oxygen vacancy neighboring Fe were analyzed. Calculated
oxygen vacancy concentration along with ionic mobility provided oxygen ion conductivity as a
function of temperature and phase change in LSF compositions. These calculations suggested the
LSF composition and phase with highest ionic conductivity at SOFC operating conditions.
In Chapter 8, anisotropic chemical expansion in cubic ceria was explained with complete
chemical strain tensor along with detailed local lattice distortion. These findings explained the
reason behind anomalous strain effect observed in doped ceria (on application of electric field)
and directional strain produced in the lattice due to the oxygen vacancy.
Finally, based on the knowledge gained and computational methods developed in this
thesis, a future research direction was proposed in Chapter 9, to find new SOFC cathode and to
understand the chemical strain effect in SOFC electrolyte like Gd/Pr/Sm doped ceria.

1.1

Factors Influencing Oxygen Vacancy Formation in Perovskite

Multiple factors like the charge distribution on the multivalent Fe cations, interactions between the
oxygen vacancies, and oxygen nonstoichiometry (δ) influence the oxygen vacancy formation in
perovskite is known for past few decades.43,44 Whereas, the mechanistic details explaining the
complicated interplay between these factors has yet to be fully established.5,6,30,45–47 In order to
promote mobile oxygen vacancy concentrations in strontium ferrite and other perovskite type
materials,28,29,48–51 it is critical to know the location of the oxygen vacancy formation site and the
corresponding formation energy and concentration, which is in turn related to the non-uniformly

6

distributed charge states on the transition metals. Unfortunately, several difficulties encountered
by past literature studies have prevented these relationships from being clearly established.32–36

First, it has been experimentally difficult to definitively identify the oxygen vacancy formation
sites within the lanthanum strontium ferrite lattice. As shown schematically in Figure 1.1 for
SrFeO3-δ (SFO) phases, X-ray44 and neutron32 diffraction studies have shown that with increasing
oxygen nonstoichiometry, the octahedral Fe-O coordination in cubic perovskite SrFeO3 first
transforms into a mixture of square pyramids and octahedra in tetragonal and orthorhombic
strontium ferrite, and finally, transform into a mixture of octahedra and tetrahedra in the
brownmillerite SrFeO2.5 phase. However, due to the unordered arrangement of the as-formed
oxygen vacancies, present experimental techniques have not able to identify their physical
distribution within the strontium ferrite crystal structures. Second, despite the fact that various
calculations have indicated that the oxygen vacancy formation energy varies with concentration in
similar materials,28,29,48,49 this effect has not been accounted for in past studies on SFO. For
instance, experimental studies have routinely obtained activation energies through linear fits of
ionic conductivity versus inverse temperature,38 assuming the activation energy is a constant,
which only holds when the vacancies are dilute. Past modeling studies have either utilized dilute
defect concentration models,52 or ignored the concentration dependence of vacancy formation
energy,24 so they cannot predict the vacancy concentration accurately in materials such as
SFO38,53,54 with interacting point defects. Thus, it is unsurprising that past modeling studies
comparing DFT-predicted oxygen vacancy formation energies55 with experimentally fitted
activation energies measured over large temperature and vacancy concertation spans have shown

7

large discrepancies.24,55,56 Because of this, accurate comparisons of either experimentally derived
activation energies or directly measured oxygen nonstoichiometry44,57–59 with the results of

Figure 1.1. Vacancy ordered phases of SrFeO3-δ: (a) cubic, (b) tetragonal, (c) orthorhombic, and
(d) brownmillerite with increasing oxygen non-stoichiometry. The energetically favorable sites for
vacancy creation are highlighted with dashed red circles. Fe-O polyhedra are colored according
to the calculated charge state on Fe, as shown in the legend, with Fe-O octahedra and square
pyramids containing Fe+3.9 or Fe+4 dark brown, Fe-O octahedra containing Fe+3.6 yellowishgreen; Fe+3.3 light purple; Fe+3 light brown, and Fe-O tetrahedra containing Fe+3.2 light orange.
combined DFT and thermodynamic models24,55 has been hard to achieve. Third, both experiments
and computations have struggled to establish the Fe3+/Fe4+ ratios that exist within strontium ferrite
as a function of oxygen non-stoichiometry. For instance, while it is known that with increasing
oxygen non-stoichiometry the formal charge state on Fe (a.k.a the oxidation state) changes from
Fe4+ in SrFeO3, to a mixture of Fe3+ and Fe4+

57

and finally to Fe3+

in SrFeO2.5,32 it is still being debated if (Fe4+)32,35 or (Fe3+)33,34,36,60 sits in the square pyramidal site
found in tetragonal and orthorhombic SrFeO3-δ.
8

It is well known that the low oxygen content of many “Sr-poor” LSF concentrations (such
as the LSF end member LaFeO3-δ (LFO) which has a large oxygen vacancy formation energy,
𝑓

∆𝐸𝑣𝑎𝑐 , of ~ 4 eV at 0 K)24,61 can be increased through the extrinsic aliovalent Sr substitution of La
via the reaction (in KrĂśger-Vink notation):
′
2𝑆𝑟𝑂 → 2𝑆𝑟𝐿𝑎
+ 2𝑂𝑂𝑥 + 𝑉𝑂

[1.2]

However, under most conditions, the extrinsic aliovalent Sr dopant level, x in La1-xSrxFeO3-δ, does
not independently determine δ because the mixed charge of Fe in La1-xSrxFeO3-δ (which varies
from 3+ in LaFeO3 to 4+ in SrFeO3) can cause δ to vary from 0 ~ 𝑥/2.61,62 In LSF, this internal
redox reaction can be written as:61
1

𝑥
•
2𝐹𝑒𝐹𝑒
+ 𝑂𝑂𝑥 → 2 𝑂2(𝑔) + 𝑉𝑂•• + 2𝐹𝑒𝐹𝑒

[1.3]

In addition, the overall oxygen non-stoichiometry, δ, is not always proportional to the mobile
𝑓

oxygen vacancy site fraction, X. For instance, the LSF end member SFO has a very low ∆𝐸𝑣𝑎𝑐 ≈
0.4 eV and hence at room temperature it possess a 𝛿 of 0.03, but the X is so low that it seldom used
for SOFC applications.20, 21 This is because as 𝛿 increases, the average distance between the
oxygen vacancies shrinks, and beyond a critical 𝛿 (or below a critical vacancy-to-vacancy
distance) oxygen vacancy interactions lead to an increase in the oxygen vacancy formation
energy,28,29 and eventually, as in SFO, oxygen-vacancy interactions induce phase transformations
that decrease X by rearranging oxygen vacancies into a non-mobile, orderly array.20 This
partitioning of 𝛿 between X and the number of non-mobile oxygen vacancies per formula unit
structurally required to form the particular crystal structure, 𝛿 0 , is summarized by the equation:
𝛿 = 𝛿 0 + (3 − 𝛿 0 )𝑋

[1.4]

While many past studies have justified the LSF conductivity maximum occurring at intermediate
LSF compositions as a compromise between doping and oxygen vacancy ordering,63 a detailed
9

understanding of the mechanistic reasons why some intermediate compositions and/or phases are
better than others has yet to emerge.
Unfortunately, past attempts to understand the relationship between the Fe charge state
distribution, the La:Sr ratio, and the oxygen vacancy interactions in LSF have been limited. For
𝑓

instance, Ritzmann et al.,24 calculated that the oxygen vacancy formation energy, ∆𝐸𝑣𝑎𝑐 , increased
in La0.75Sr0.25FeO3-δ, decreased in La0.5Sr0.5FeO3-δ, and remained a constant in LFO, when δ
increased from 1/32 to 1/8. However, their interesting δ-dependent oxygen vacancy formation
energies did not contain an analysis of the oxygen vacancy interactions or the Fe oxidation states
𝑓

in LSF. Further, Mastrikov et al.,55 showed that ∆𝐸𝑣𝑎𝑐 increased with higher La:Sr ratio in LSF
𝑓

(from SFO to La0.5Sr0.5FeO3-δ (LSF55) to LFO) at a fixed δ = 1/8. This increase of ∆𝐸𝑣𝑎𝑐 with
increasing La:Sr ratio was qualitatively related to the oxidation state on Fe but the reasons for the
measured δ peak around a La:Sr ratio of 0.6:0.4 were not discussed. Furthermore, even though the
experimental results of Dann et al.,64 Patrakeev et al.,62 and Fossdal et al.,65 have confirmed that
the La1-xSrxFeO3-δ crystal structure evolves from orthorhombic (0 ≤ x ≤ 0.2) to rhombohedral (0.4
≤ x ≤ 0.7) to cubic (0.8 ≤ x ≤ 1.0) at room temperature, and Fossdal et al.65 observed a rhombohedral
to cubic second-order phase transition when rhombohedral LSF55 was heated from room
temperature to 523Âą50 K (in air), these phase changes were not well captured in previous DFTbased oxygen vacancy calculations. For example, Ritzmann et al.24 and Mastrikov et al.55 used a
cubic supercell for all LSF structures where x ranged from 1/32 to 1/8, even though experiments
have shown that LSF is not cubic at those La:Sr ratios. Recently, Taylor et al. used an orthorhombic
supercell for LFO but unfortunately did not explore other La:Sr ratios.66
Therefore, the aim of the present study was to correctly establish the underlying
mechanistic relationships between the La:Sr ratio, the Fe charge state, the oxygen vacancy polaron
10

size and shape (discussed in Section 1.2), δ and X. Here, three different LSF compositions, LaFeO3δ

(LFO), La0.5Sr0.5FeO3 (LSF55) and SrFeO3-δ (SFO) were studied. SFO and LFO represent the

parent phases65 of LSF and LSF55 was selected for its high oxygen ionic conductivity62,63 and as
a representative structure of the commonly used SOFC cathode material LSF64 (La0.6Sr0.4FeO36,61

δ).

The study objectives were achieved using a combined thermodynamics and DFT+U

approach performed by first, appropriately selecting a suitable U parameter (as described in
Sections 2.1.1), second, using a magnetic moment interpreted charge calculation, instead of the
traditional and in this case insensitive Bader charge analysis, to determine the Fe charge state (as
described in Section 2.1.2), third, validating the DFT simulations with known experimentally
determined physical properties (as described in Chapter 3), fourth, calculating the effect of oxygen
𝑓

vacancy interactions on ∆𝐸𝑣𝑎𝑐 (as described in Chapter 4), and fifth, using a newly-developed
thermodynamic model27 to predict X for various LSF compositions and phases under SOFC
operating conditions (described in Chapter 6).

1.2

Oxygen Vacancy Interaction via Vacancy Polaron

MIEC oxygen vacancy ordering is often discussed in terms of oxygen vacancy polarons, which
Landau67 defined as the excess charge (generated during the formation of a charge neutral oxygen
vacancy) and the surrounding lattice strain field. Since the excess electrons generated during
oxygen vacancy formation both screen the charge on their associated oxygen vacancy site and
induce lattice distortions, they impact the formation of subsequent oxygen vacancies at high δ.
While it is hard to exactly define polaron size and shape because the partial excess charges are
distributed to both the Fe and O atoms around a missing lattice oxygen,68,69 for simplicity, one can
estimate the polaron size and shape based on the spatial deviations in transition metal charge states
11

alone. For example, Castleton et al.,69 defined the oxygen vacancy polaron in ceria as the localized
electrons on two Ce atoms adjacent to an oxygen vacancy site.70 In a recent study Das et al.,27 also
took this approach and found that the strong oxygen vacancy interactions in LSF end member,
SFO was due to a very large oxygen vacancy polaron size. Specifically, it was found that electron
occupancy and crystal-field-splitting-induced differences between the 3d electron orbitals of the
square pyramidally coordinated, oxygen-vacancy-adjacent Fe atoms and the octahedrallycoordinated, oxygen-vacancy-second-nearest-neighboring Fe atoms caused the two electrons left
behind during oxygen vacancy formation in cubic SFO to spread to the second nearest neighbor
Fe (instead of the Fe directly connected to the vacancy) resulting in extended lattice distortions
that promoted strong oxygen vacancy interactions.27 However, the distribution of the excess
electrons generated by oxygen vacancy formation in La1-xSrxFeO3-δ has not been examined. Further,
a detailed understanding of the relationship between the Fe charge state distribution associated
with polaron and oxygen vacancy interactions in LSF has been missing from the literature. This is
partially due to the difficulty in determining the actual Fe charge distribution within LSF, for both
experimentalists and theorists.24,43,61,62,71 For instance, while it is recognized that La1-xSrxFeO3-δ
has a mixture of both Fe3+ and Fe4+ ions, it is still debated where Fe3+ or Fe4+ sits in vacancycontaining lattices.
This new “long-range charge transfer” mechanism results in large oxygen vacancy polaron
sizes that cause strong oxygen vacancy interactions and an increase in the oxygen vacancy
formation energy with increasing oxygen vacancy site fraction. The objective of the present work
was to examine the impact of lanthanum doping on this SFO oxygen vacancy formation behavior
and the oxygen vacancy migration behavior of a variety of LSF compositions and crystal structures.

12

Further, use this charge transfer concept to predict oxygen vacancy interactions in Mn or Cocontaining perovskites, La1-xSrxMnO3-δ or La1-xSrxCoO3-δ.

1.3.

Effect of the oxygen vacancy interaction on Ionic conductivity at SOFC operating
conditions

Performance of the MIEC materials is mostly limited due to an exponential loss in oxygen ionic
conductivity at a lower temperature.72 To achieve the ultimate goal of designing MIEC materials
with high ionic conductivity operating at a relatively lower temperature, we have to understand
the factors (mentioned earlier) controlling the ionic conductivity in it. The specific aim of this
study is to predict the oxygen ionic conductivity of the SOFC material from computationally
predicted oxygen vacancy concentration and diffusivity. The ionic conductivity, 𝜎 can be
predicted with Equation 1.1. Where, c is nothing but δ/Vu, (Vu: volume per formula unit) if there
is no vacancy ordering in the lattice. In case of oxygen vacancy ordering, c should be related to X
with below equation:
𝑐 = (3 − 𝛿 0 )𝑋/𝑉𝑢

[1.5]

The above equation holds for both vacancies ordered and unordered lattices. For lattice without
vacancy ordering, 𝛿 0 = 0. Therefore, in that case δ = 3X, which satisfies previous relation.
Further, it is required to develop a thermodynamic model which can predict oxygen
vacancy site fraction as a function of temperature, pressure, and oxygen vacancy concentration
and then use that data to predict ionic conductivity. The model should be robust enough that it will
efficiently predict the oxygen site fraction for both interacting and non-interacting vacancies. The
thermodynamic model prior to this work used to predict vacancy site fraction is shown in Equation

13

𝑓

1.6.52 Where, ∆𝐺𝑣𝑎𝑐 is the vacancy formation free energy at the SOFC operating temperature (T)
and partial pressure of oxygen (p). 𝑘𝐵 denotes Boltzmann constant.

𝑋

= exp (−
1−𝑋

𝑓

∆𝐺𝑣𝑎𝑐 (𝑇,𝑝)
𝑘𝐵 𝑇

)

[1.6]

Hence, it is required to develop a thermodynamic model as shown in Section 2.5, which can predict
oxygen vacancy site fraction as a function of temperature, pressure and concentration dependent
formation free energy.
Einstein relation can relate the mobility of an ion to diffusivity (D), is given as,
𝑍𝑒𝐷

𝜇=𝑘

[1.7]

𝐵𝑇

Furthermore, the diffusivity can be calculated as:
𝐸

𝐷 = 𝐷0 exp (− 𝑘 𝑎𝑇)
𝐵

[1.8]

Where, 𝐷0 is the pre-factor of diffusion. 𝐸𝑎 is the activation energy for diffusion can be obtained
𝑓

from vacancy formation free energy, ∆𝐺𝑣𝑎𝑐 and migration energy, ∆𝐻𝑚𝑖𝑔 . ∆𝐻𝑚𝑖𝑔 can be
calculated via nudged elastic band method (NEB). If this method is established to predict ionic
conductivity at SOFC operating conditions, then we can directly compare the DFT predicted
results with experimental measurements. Hence, this prediction will elucidate a method to directly
predict ionic conductivity for interacting vacancies of newly designed materials without any
experimental inputs.

14

1.4.

Anisotropic Chemical Strain Effect Due to Oxygen Vacancy Polaron

Pure and doped Ceria are amongst the most studied nonstoichiometric oxides, due to its application
as a commercial catalyst.73,74 The catalytic capability of ceria owes to the reversible transition
between Ce4+ and Ce3+. The nonstoichiometric Gd-doped ceria (GDC) though mainly used for
catalytic application but the reversible transition of Ce4+ and Ce3+ enables its’ application for
electrostriction.75 The application of electric field on a nonstoichiometric GDC thin film orients
the local chemical strain fields (with polaron) producing observable mechanical stress, is a good
example of the electro-chemo-mechanical coupling concept.76 The change of chemical strain field
on the application of electricity implies the anisotropic nature of chemical strain in ceria, which is
the prime focus of this study.
Previous experimental and computational studies were devoted to measure/calculate the
average chemical strain effect in pure or doped ceria. For example, using molecular dynamics
simulation Marrocchelli et al.40 predicted average chemical expansion coefficient, 𝛼𝐶 , by creating
different amount of vacancies (𝛿) and predicting average strain (𝜀𝐶 ) at those 𝛿. Then from the
slope of 𝜀𝐶 with 𝛿, 𝛼𝐶 (= 𝜀𝐶 /𝛿) was calculated. Further, using DFT Shenoy et al.41 calculated
isotropic 𝜀𝐶 , assuming isotropic chemical strain field in ceria, due to its cubic symmetry. Shenoy
et al. used only principle components of the elastic dipole tensor to calculate average chemical
strain, neglecting anisotropic effect. Whereas, the orientation of chemical strain dipole on
application of electric field suggests that, due to non-uniform bond distortion around oxygen
vacancy site, the oxygen vacancy induced chemical strain is anisotropic in nature, which helped
to orient the strain field on application of electricity.
The objective is to understand how the oxygen vacancy polaron causes asymmetric lattice
distortion locally in cubic ceria, thus anisotropic chemical strain is produced in nonstoichiometric

15

ceria. To understand the oxygen vacancy induced anisotropic chemical strain due to non-uniform
bond distortion around oxygen vacancy site, a complete strain tensor, 𝜺𝑪 should be used. The
complete strain tensor can be calculated from elastic dipole tensor, G and elastic stiffness tensor,
ℂ as:42
𝜺𝑪 =

−𝛿(ℂ−1 𝑮)

[1.9]

𝑉𝑈

Where δ is oxygen nonstoichiometry, and 𝑉𝑈 is the volume per formula unit. The detailed
derivation of the above equation is provided in Section 2.6. With this method, the reason behind
anisotropic chemical strain tenor was quantified and its tensorial relationship with the
crystallographic orientation of the lattice will explain the giant strain effect in GDC under electrical
bias.42

16

2.

Computational Details: Methodology developed and Used

The plane wave based ab initio simulation package, VASP (Vienna Ab initio Simulation
Package)77 was used for all the ground state (0 K, zero pressure) calculations. Projectoraugmented-wave (PAW)78,79 potentials with valence configurations of 5s25p65d16s2 for La,
4s24p65s2 for Sr, 3d74s1 for Fe, and 2s22p4 for O were used to describe the valence electrons. The
generalized gradient approximation (GGA) functional along with Perdew, Burke, and Ernzerhof
(PBE)80,81 parameters were used to describe the exchange-correlation potentials of the constituting
elements. As GGA predicted results suffer from self-interaction error in strongly correlated
transition metal complexes,82 Fe was treated with the GGA+U method83 with a Ueff = 3. A detailed
discussion of the Fe atom U parameter selection process is provided in Sections 2.1.1. Spinpolarized calculations were performed and the magnetic moment on each Fe was calculated by
spherical integration.84 The accuracy of the electronic calculations was within 1 ÎźeV. A single
point energy calculation was performed to test the convergence of the cut-off energy and k-points
for each input lattice. It was concluded from the k-mesh calculations that a k-spacing of 0.2 Å-1
was sufficient for the required accuracy of 1 meV/atom when the cutoff energy was set to 500 eV.
Ionic relaxations were performed until the Hellmann-Feynman force on each atom reached the
order of 10 meV/ Å.85
First, the calculation methods were tested for their ability to accurately predict the
experimentally observed LSF crystal structures, lattice parameters, and low-temperature Fe
magnetic ordering arrangements. According to experiments, the representative LFO, SFO, and
LSF55 compositions exhibit orthorhombic, cubic, and rhombohedral structures, respectively, at
room temperature.44,46,86 In addition to the different crystalline structures shown in Figure 1, they
can also exhibit magnetic ordering, as indicated by MĂśssbauer results showing G-type Fe

17

antiferromagnetic (AFMG) ordering for 0< x< 0.3 and paramagnetic Fe ordering for x > 0.4, at
room temperature.87

The ground state energy for all three oxygen-stoichiometric (i.e., δ = 0) LSF compositions
was calculated with six different initial structures (three ground state crystal structures, as shown
in Figure 2.1, multiplied with two magnetic ordering arrangements), by allowing the cell shape
and atom positions to relax. Among these six input structures, the one with the minimum energy
per formula unit was considered as the (ground state) crystal structure for that composition. The
calculated crystal structures, lattice parameters, bond lengths, low-temperature magnetic ordering
arrangements, and band gaps for stoichiometric SFO, LSF55, and LFO were benchmarked with
experimental or computational studies whenever available in the literature.

Figure 2.1. DFT calculated lattice structures for (a) orthorhombic LaFeO3, (b) rhombohedral
La0.5Sr0.5FeO3, (c) cubic La0.5Sr0.5FeO3 (high temperature phase), and (d) cubic SrFeO3. Atoms
are colored as La – green, Sr - blue, Fe – brown, O – red. The different Fe-O6 octahedra colors
denote polyhedra with different Fe charge states: Fe3+ – light brown, Fe3.5+ – moss green, Fe4+ –
dark brown.
2.1.

Magnetic Moment Interpreted Oxidation State Prediction

Though Bader charge analysis is a widely used and efficient method to define formal charge on
atoms,88,89 it is not universally applicable for all compounds.48,53 For example, Bader charge
analysis on Fe in LaFeO3, La0.75Sr0.25FeO3-δ and La0.5Sr0.5FeO3 show a charge of 1.70 e, not able
to determine the charge change on Fe due to change in La/Sr-ratio. Further, after oxygen vacancy
18

formation the charge on vacancy adjacent Fe was not significantly different for changing La/Srratios.24 Therefore, a new way to define charge on Fe was required. As a result, we propose a
magnetic moment interpreted charge for a magnetic element, like Fe. As magnetic moment is
directly related to a number of unpaired electrons, so from magnetic moment on Fe formal charge
can be predicted and can be compared with experimental measurement. Only DFT study was not
able to predict the magnetic moment on Fe correctly (refer to U = 0 in Figure 2.2(a)), and HubbardU was necessary for proper treatment of the localized 3d-electrons of Fe.

2.1.1. Selection of U Parameter

DFT with a Hubbard-U approach (DFT+U)90 was used to address the on-site electronic interactions
of the 3d-orbitals in the multivalent Fe atoms. This method has been shown to accurately predict
the enthalpy of formation of oxides, the open circuit voltage in Li-ion batteries, phase diagrams,
chemical reaction barriers, etc. for many transition-metal-containing compounds.30,91–93 In the
DFT+U framework, spherically averaged Coulombic and exchange interactions for the electrons,
denoted by U and J respectively, enter into DFT calculations as overhead terms.94 Following the
approach of Dudarev et al.,90 these U and J parameters can be combined as Ueff = U – J, or used
simply as a U parameter. The Hubbard-U parameter can be selected either directly employing ab
initio methods like constrained-DFT94,95 or based on empirically fitting experimental properties
such as the band gap, magnetic moment, lattice parameters, bulk modulus, etc.3,90,92,96,97 In some
DFT+U studies, different transition metal charge states were treated with different U parameters.
For example, Mosey et al.,95 treated Fe2+ and Fe3+ with Ueff = 3.7 and 4.3 respectively, giving
reasonable lattice parameter, magnetic moment, band gap and other physical property agreement
for Fe2O3 and FeO. However, the energies calculated with different U parameters should not be
19

compared (since the electrons are more localized with increasing U), making it impossible to use
that approach to perform comparative studies of different LSF phases containing mixed Fe
oxidation states.
In the present study, because the local formal charge state on Fe was allowed to vary from
2+ ~ 4+ (due to the changing La/Sr ratio and the changing oxygen vacancy content), a single U
parameter was required to treat all the multi-valent Fe atoms. GGA+U calculations were performed
for (orthorhombic, FM) LFO and (cubic, FM) SFO with the U parameter over a range of 0.2 to 4.0
eV. The U-parameter was selected based on a careful comparison with the experimentally
observed magnetic moments (as shown in Figure 2.2(a)) and lattice parameters (as shown in Figure
2.2(b)) of the LSF end members. The reason behind the selection of U to match both LFO and
SFO was that the charge on Fe in LSF varies between 3+ and 4+, and the charge on Fe in LSF can
be further interpreted as a linear mixing of 3+ and 4+ (as shown in Figure 2.3).98

Figure 2.2. The plot of (a) magnetic moment and (b) lattice parameters as a function of Ueff
parameter. The straight lines are experimental data, and blue and red correspond to the data for
SrFeO3 and LaFeO3, respectively. Since LaFeO3 is not cubic, red, green and pink colors in (b)
are corresponding to scaled-lattice (a/√2, b/2, c/√2 respectively) of LaFeO3.

20

Figure 2.2 (cont’d)

Figure 2.2(a) shows a comparison between the experimental and computed Fe magnetic
moments with varying U-parameter (from 0.2 to 4) for LFO and SFO, represented by lines and
solid dots, respectively. SFO showed a steady increase in the magnetic moment of Fe with
increasing U parameter but the magnetic moment of Fe in LFO fluctuated (Figure 2.2(a)), possibly
due to the existence of additional minima as a result of different possible electronic distributions.99
The experimental magnetic moment for Fe in LFO was reported by Koehler et al.100 as 𝜇𝐹𝑒 3+ =
4.6 ± 0.2 μB (Figure 2.2(a), dashed line) and it is a widely acceptable value in the literature.101–103
In contrast, difficulty in the experimental sample preparation of stoichiometric cubic SFO104,105
has resulted in a large variability in the experimentally measured Fe magnetic moments in
SFO.32,105,106 Nevertheless, most of the experimental data indicates that Fe has a high spin state in
SFO.104,107–109 In a recent comprehensive study, Reehuis et al.105 measured the magnetic moment
on Fe in stoichiometric cubic SFO as 𝜇𝐹𝑒 4+ = 2.96 μB (Figure 2.2(a), dotted line) with neutron
diffraction, at 2 K.
21

Figure 2.2(b) compares the computed lattice parameters of LFO and SFO with
experimental results. Note, to allow a comparison of the computed and experimental lattice
parameters in a single plot, all the LFO lattice parameters were divided by √2 before plotting in
Figure 2.2(b). As done previously for SFO,27 based on a comparison of the computed versus
experimental lattice parameters and magnetic moments for both LFO and SFO in Figure 2.2, a Ueff
= 3 was selected. With this value of U, the lattice parameters showed a maximum deviation of
1.75 % from their experimentally-measured room temperature values. The selection of Ueff = 3
also resulted in a reasonable agreement between the computed and experimentally measured
magnetic moments. The computed magnetic moment was 𝜇𝐹𝑒 3+ = 4.23 μB for Fe3+ in LFO, and
was 𝜇𝐹𝑒 4+ = 3.61 μB for the Fe4+ present in SFO, both of which were similar to the experimental
values.100,101,105 In some earlier computational studies with multivalent Fe atoms, the focus was on
Fe3+ containing compounds, for which a Ueff = 4 110 or 5.4 111 was employed. However, only a few
studies exist where Fe multi-valence states were treated with a single U parameter. For example,
Muñoz-García et al.,49 treated 2+, 3+, and 4+ charge states on Fe with U eff = 4 and Ritzmann et
al.,24 used Ueff = 5.4, where the Fe charge state varied from 3+ to 3.5+. One consequence of using
such a high U-parameter is that the vacancy formation energy might be underestimated.27 For
example, it was reported in our previous study that Ueff = 4 would lead to a lower oxygen vacancy
formation energy in SFO, resulting in a ~200˚C difference in the temperature required to generate
the oxygen non-stoichiometry observed in cubic SrFeO3-δ in the air.27 Thus, when Ritzmann et
al.,24 used Ueff = 5.4 in their calculations, the underestimated vacancy formation energy led to an
unphysical (negative) oxygen vacancy formation free energy at 700oC for LSF55. As will be
shown in Chapter 3, a Ueff =3 also gave the correct band gap and phases for various SFO
compositions and phases.

22

Figure 2.3. Black squares represent the formal charge on Fe resulting from different La/Sr ratios
in LaxSr1–xFeO3 calculated from the magnetic moment on Fe with a GGA+U(=3) (left y-axis). Blue
diamonds show experimentally measured average Fe charge states for varying La/Sr ratios (right
y-axis). The black dotted line represents the predicted linear relationship between the Fe charge
state and the La/Sr ratio.
2.1.2. Determining Oxidation State on Fe
DFT with a Hubbard-U approach (DFT+U)90 was used to address the on-site electronic interactions
of the 3d-orbitals in the multivalent Fe atoms. This method has been shown to accurately predict
the enthalpy of formation of oxides, the open circuit voltage in Li-ion batteries, phase diagrams,
chemical reaction barriers, etc. for many transition-metal-containing materials.30,91–93 In the
DFT+U framework, spherically averaged Coulombic and exchange interactions for the electrons,
denoted by U and J respectively, enter into DFT calculations as overhead terms.94 Following the
approach of Dudarev et al.,90 these U and J parameters can be combined as Ueff = U – J, or used
simply as a U parameter. The Hubbard-U parameter can be selected either directly employing ab
initio methods like constrained-DFT94,95 or based on empirically fitting experimental properties
such as the band gap, magnetic moment, lattice parameters, bulk modulus, etc.3,90,92,96,97 In some
DFT+U studies, different transition metal charge states were treated with different U parameters.

23

For example, Mosey et al.,95 treated Fe2+ and Fe3+ with Ueff = 3.7 and 4.3 respectively, giving
reasonable lattice parameter, magnetic moment, band gap and other physical property agreement
for Fe2O3 and FeO. However, the energies calculated with different U parameters should not be
compared (since the electrons are more localized with increasing U), making it impossible to use
that approach to perform comparative studies of different LSF phases containing mixed Fe
oxidation states.
In the present study, because the local formal charge state on Fe was allowed to vary from 2+ ~
4+ (due to the changing La/Sr ratio and the changing oxygen vacancy content), a single U
parameter was required to treat all the multi-valent Fe atoms. GGA+U calculations were
performed for (orthorhombic, FM) LFO and (cubic, FM) SFO with the U parameter over a range
of 0.2 to 4.0 eV. The U-parameter was selected based on a careful comparison with the
experimentally observed magnetic moments (as shown in Figure 2.2(a)) and lattice parameters (as
shown in Figure 2.2(b)) of the LSF end members. The reason behind the selection of U to match
both LFO and SFO was that the charge on Fe in LSF varies between 3+ and 4+, and the charge
on Fe in LSF can be further interpreted as a linear mixing of 3+ and 4+ (as shown in Figure 2.3).98

2.1.3. Advantage over Bader Charge
The calculated total Bader charge analysis, as shown in Table 2.1, indicates charge on Fe in LFO,
LSF55 and SFO were 2.00+, 2.03+ and 2.04+ respectively, whereas the expected formal charge
on Fe in these structures are 3.0+, 3.5+, and 4.0+ respectively. This signifies that the calculated
total Bader charge on Fe was insensitive to the local electronic environment in LSF, and therefore
could not be used to determine the formal charge on Fe. Fortunately, the total unpaired electrons
obtained via spin polarized Bader charge analysis of Fe show a significant variation from 3.88 to

24

4.64 e from SFO to LFO. Hence, the magnetic moment, which is a manifestation of the unpaired
electrons in an element, was used to interpret the charge state on Fe. Accordingly, the magnetic
moment calculated by VASP was used to interpret the formal charge state on Fe. Note, neither the
spherical integration of charge used in a VASP calculation nor the Voronoi integration used in the
Bader method gave the exact formal charge on Fe in LSF. Therefore, a magnetic moment of 4.23
ÎźB for Fe was assigned to represent Fe3+ and 3.61 ÎźB for Fe4+ in agreement with the LSF end
member Fe charge states and magnetic moment discussed in Section 2.1.
Table 2.1. Spin-polarized Bader charge analysis of different LSF phases. Magnetic moment
calculated by VASP and expected formal charge on Fe are also provided for comparison.

It is possible that after neutral oxygen vacancy generation, some oxygen-vacancy-neighboring Fe3+
atoms will adopt a charge closer to Fe2+. Although Fe2+ was not considered in the U-parameter
selection, its charge was calibrated with the magnetic moment of Fe2+ in (planar) SrFeO2. Since
Fe2+ has the same number of unpaired electrons as Fe4+, it gives rise to a magnetic moment of 3.68
μB which is close to the 𝜇𝐹𝑒 4+ = 3.61 μB (and comparable to experimentally reported SrFeO2
magnetic moment at 10 K).112
25

Figure 2.3 shows the computed magnetic moment and the formal charge on Fe compared to the
experimentally measured average charge on Fe in LSF, which also shows a linear relationship.
Hence, in the current work, an intermediate magnetic moment was interpreted as a linear
interpolation between 3+ and 4+ charge states.98 In LSF55, Fe is expected to carry a formal charge
of 3.5+ and from the linear relation of charge and magnetic moment, Fe3.5+ corresponds to a
magnetic moment of 3.92 ÎźB, which was very close to the VASP calculated magnetic moment of
3.94 ÎźB for Fe3.5+ in rhombohedral LSF55. Overall, the selection of Ueff = 3 and the magnetic
moment interpreted charge jointly captured the Fe charge state ranging from Fe2+ to Fe4+ in LaxSr1–
xFeO3-

2.2.

.
Determining Polaron Shape and Size from the Valance Change on Transition Metal

In 1933 Lev Landau,67 proposed the concept of polaron as a quasiparticle. It was described as an
electron moving in a dielectric crystal where the atoms move from their equilibrium positions to
effectively screen the charge of an electron. Following the concept we have defined our oxygen
vacancy polaron, where two excess electron left due to the oxygen vacancy formation gets
distributed to the nearest neighboring Fe atoms, causing charge change and lattice distortion. For
accuracy of the calculation if the magnetic moment interpreted charge on a Fe changed by Âą 0.1e
after oxygen vacancy formation, then it was considered as the “charge change” on Fe. Further,
from charge change on Fe atoms neighboring to oxygen vacancy site, we found the extent of charge
transfer and shape and size of the polaron. If the excess electron generated due to a neutral oxygen
vacancy formation transferred to second nearest neighboring Fe, then it was named as “long range
charge transfer”.

26

2.3.

Calculation Method of Oxygen Vacancy Formation and Migration Barrier Energy

Oxygen vacancy formation and migration energies were calculated using DFT+U method. The 0
K calculated energy with temperature and pressure correction can be converted to experimentally
measurable quantity like oxygen vacancy nonstoichiometry, diffusivity etc. The vacancy
formation and migration energies were calculated with varying size of supercell for each LSF
phases. If the calculated energy didn’t change with lattice size then it was considered as dilute
vacancy scenario, otherwise interacting vacancy scenario.
2.3.1. Oxygen Vacancy Formation Energy at 0 K

1

For the defect reaction,113 𝑂𝑂𝑥 ⇌ 𝑉𝑂•• + 2𝑒 ′ + 2 𝑂2, the neutral oxygen vacancy formation energy
𝑓

(∆𝐸𝑣𝑎𝑐 ) per vacancy at 0 K was calculated using the relationship:85
1

𝑓

𝐷𝐹𝑇
𝐷𝐹𝑇
∆𝐸𝑣𝑎𝑐 = 𝐸𝑑𝑒𝑓𝑒𝑐𝑡𝑖𝑣𝑒
− 𝐸𝑝𝑒𝑟𝑓𝑒𝑐𝑡
+ 2 𝐸𝑂𝐷𝐹𝑇
2

[2.1]

𝐷𝐹𝑇
𝐷𝐹𝑇
where 𝐸𝑑𝑒𝑓𝑒𝑐𝑡𝑖𝑣𝑒
is the energy of a structure with one oxygen vacancy, 𝐸𝑝𝑒𝑟𝑓𝑒𝑐𝑡
is the energy of

the perfect crystal structure, and 𝐸𝑂𝐷𝐹𝑇
is the energy of an oxygen molecule114 in its gas phase.
2
𝐸𝑂𝐷𝐹𝑇
was computed as -9.86 eV for an isolated oxygen molecule placed inside a 202020 Å3
2
simulation cell, with a predicted bond length of 1.23 Å. The predicted energy and bond length are
𝐷𝐹𝑇
in well agreement with values reported in the literature.85,115,116 𝐸𝑑𝑒𝑓𝑒𝑐𝑡𝑖𝑣𝑒
was calculated in

supercells of varying size containing a single oxygen vacancy, mimicking a broad range of oxygen
nonstoichiometry. In SrFeO3-δ, increasing the absolute oxygen non-stoichiometry (𝛿) causes the
oxygen vacancy site fraction (𝑋) to increase until oxygen vacancy ordering induces a phase
transformation that reduces the total number of lattice oxygen sites and resets 𝑋 to zero; a fact
which can be summarized with the Equation 1.4. Where 𝛿 0 is the value of 𝛿, when 𝑋 = 0 for each
27

perfect crystal structure (i.e. 𝛿 0 equals 0, 0.125, 0.25, and 0.5 in cubic SrFeO3, tetragonal
SrFeO2.875, orthorhombic SrFeO2.75, and brownmillerite SrFeO2.5, respectively). Hence, a single
oxygen vacancy created in a supercell with n formula units of cubic strontium ferrite (SrnFenO3n)
corresponded to a δ = 1/n, while a single oxygen vacancy created in a supercell with n formula
units of brownmillerite strontium ferrite (SrnFenO2.5n) corresponded to a δ=0.5+1/n.
It is important to note that as shown in Equation 1.4, the same 𝛿 can be achieved by having a large
𝑋 in a structure with a low 𝛿 0 or by having a small 𝑋 in a structure with a higher 𝛿 0 . For crystal
structures containing multiple oxygen vacancy sites (Figure 1.1(b-d)), the oxygen vacancy
formation energy was calculated for all possible sites in order to identify the site with the lowest
oxygen vacancy formation energy. The largest supercell used for the cubic SFO, rhombohedral
LSF55, cubic LSF55 and orthorhombic LFO calculations was 4x4x4, 3x3x3, 3x3x3, and 2x2x2,
which represents an average distances of 15.52 Å, 16.55 Å, 16.58 Å, 12.78 Å between the
vacancies for δ = 0.02, 0.02, 0.02, and 0.03 respectively. The present study focuses on the variation
of non-stoichiometry at a uniform distribution of the oxygen vacancies. DFT+U minimized lattice
parameters were used in vacancy formation calculations, since it is estimated that the thermal
expansion will only induce lattice parameter changes less than 3% between 0 and 1000 K.106
Further, it is important to note that the material undergoes a continuous phase change in the
measured temperature range.

2.3.2. Predicting Vacancy Ordered Phase Transition
It is known that DFT+U calculation does a poor job in predicting phase transition in materials,
especially when metal to insulator phase transition happens but our selection of U parameter in
this study was able to capture the energetically favored phases when oxygen vacancy ordered phase
28

transitions occurred at a certain oxygen non-stoichiometry. In order to determine energetically
favorable arrangements, the excess energy (ΔEδ) of various oxygen deficient structures (cubic,
tetragonal, orthorhombic, and brownmillerite) with respect to the perfect cubic SrFeO3 was defined
and compared via the relation:
𝛿

𝐷𝐹𝑇
𝐷𝐹𝑇
∆𝐸 𝛿 = 𝐸𝑆𝑟𝐹𝑒𝑂
+ 2 𝐸𝑂𝐷𝐹𝑇
− 𝐸𝑆𝑟𝐹𝑒𝑂
2
3
3−𝛿

[2.2]

𝐷𝐹𝑇
where 𝐸𝑆𝑟𝐹𝑒𝑂
is the energy of the oxygen deficient arrangement of interest, 𝐸𝑂𝐷𝐹𝑇
is the energy
2
3−𝛿
𝐷𝐹𝑇
of an isolated oxygen molecule, and 𝐸𝑆𝑟𝐹𝑒𝑂
is the energy of the perfect cubic SrFeO3.
3

2.3.3. Calculating Oxygen Vacancy Migration Barrier
The oxygen migration barriers between all the possible oxygen vacancy sites in cubic,
tetragonal, and orthorhombic SFO, rhombohedral LSF55, and orthorhombic LFO were calculated
and compared at 0 K. Then, the impact of a change in the magnetic moment of the Fe atoms along
the oxygen vacancy migration path was determined. The migration energy barrier can be converted
to diffusivity using thermodynamic data as shown in Section 2.5.2. The migration energy can be
calculated using Climb Image Nudge Elastic Band (CINEB) method117 which will also provide
the transition state structure. In the NEB method, one only gets the minimum energy path
(MEP).118,119 The maximum value of the potential energy along the path will give the activation
energy (barrier) for the diffusion process. Using CINEB it not only gives the MEP but also finds
the transition state structure at the saddle point.120

29

2.4.

Thermodynamic Extension of DFT Calculated Energies

DFT calculations are performed at 0 K and zero pressure. Therefore, to compare DFT calculated
properties with experimentally observed results a thermodynamic extension to the DFT calculated
properties is necessary.

2.4.1. Converting Vacancy Formation Energy to Vacancy Formation Free Energy
To facilitate comparisons between computational predictions and experimentally measurable
properties (such as the oxygen vacancy concentration at finite temperatures), DFT computed
oxygen vacancy formation energies were converted into temperature (T) and oxygen partial
𝑓

pressure (𝑝𝑂2 ) dependent oxygen vacancy formation free energies (∆𝐺𝑣𝑎𝑐 ) via the relation:
1

𝑓

∆𝐺𝑣𝑎𝑐 (𝑇, 𝑝𝑂2 ) = 𝐸𝑑𝑒𝑓𝑒𝑐𝑡𝑖𝑣𝑒 (𝑇, 𝑝𝑂2 ) − 𝐸𝑝𝑒𝑟𝑓𝑒𝑐𝑡 (𝑇, 𝑝𝑂2 ) + 2 𝜇𝑂2 (𝑇, 𝑝𝑂2 )

[2.3a]

where 𝐸𝑑𝑒𝑓𝑒𝑐𝑡𝑖𝑣𝑒 (𝑇, 𝑝𝑂2 ), 𝐸𝑝𝑒𝑟𝑓𝑒𝑐𝑡 (𝑇, 𝑝𝑂2 ) and 𝜇𝑂2 (𝑇, 𝑝) are the energy of a structure with a single
oxygen vacancy, the energy of the perfect structure, and the chemical potential of oxygen gas,
respectively, at a particular T and 𝑝𝑂2 . Here, 𝜇𝑂2 (𝑇, 𝑝) is defined via the relation:
𝜇𝑂2 (𝑇, 𝑝𝑂2 ) = 𝐸𝑂𝐷𝐹𝑇
+ ∆𝜇𝑂02 (𝑇) + 𝑘𝑇𝑙𝑛𝑝𝑂2
2

[2.3b]

where 𝐸𝑂𝐷𝐹𝑇
is the energy of an oxygen molecule, and the last two terms describe the change in the
2
chemical potential of oxygen from 0 K to T , and under the partial pressure of 𝑝𝑂2 . These two terms
are collectively referred to as the thermodynamic connection energy to the DFT calculated energy
for isolated molecular oxygen. Equation 2.1 is the 0 K version of Equation 2.3a. Equations 2.3a
and 2.3b can then be combined with Equation 2.1 to yield:
𝑓

𝑓

1

∆𝐺𝑣𝑎𝑐 = ∆𝐸𝑣𝑎𝑐 + 2 (∆𝜇𝑂02 (𝑇) + 𝑘𝑇𝑙𝑛𝑝𝑂2 )

30

[2.4]

where the free energy contribution due to the phonon vibration of the extra oxygen atom in the
perfect lattice is lumped into ∆𝜇𝑂02 (𝑇).121
Defining the chemical potential of an individual atomic species in a crystalline compound
is inherently difficult. In this case, the chemical potential of oxygen in crystal lattice was set by
the temperature and pressure of the gas phase O2, since it is in chemical equilibrium with lattice
oxygen.122 The energy of oxygen dimer, 𝐸𝑂𝐷𝐹𝑇
, and the subsequent oxygen vacancy formation
2
energies are sensitive to the exact choice of exchange-correlation functional, pseudopotential
etcetera in DFT calculations. This problem is mitigated by various methods of estimating ∆𝜇𝑂02 (𝑇)
in the literature.121–124 In this paper, the method of Zhang et al.121 was modified to describe oxygen
vacancy formation in materials containing Fe-O polyhedra.
Using the derivation provided in Section 2.4.2, ∆𝜇𝑂02 (𝑇) was broken into a thermodynamic
connection energy from 0 K to room temperature (∆𝜇𝑂02 (𝑇𝑟 )), and a thermodynamic connection
energy from room temperature (𝑇𝑟 = 298 K) to a temperature of interest, T, using the expression:
𝑇

𝑇

𝑟

𝑟

∆𝜇𝑂02 (𝑇) = ∆𝜇𝑂02 (𝑇𝑟 ) + ∫𝑇 𝐶𝑝 𝑑𝑇 − 𝑇 ∫𝑇

𝐶𝑝
⁄ 𝑑𝑇 − (𝑇 − 𝑇𝑟 )𝑆𝑟
𝑇

[2.5]

where Cp is the specific heat of oxygen gas at temperature T, and Sr is the entropy of oxygen gas
at the standard state.125 It is important to note that the last three terms in Equation 2.5 are easily
calculated from the data available in any thermochemistry handbook.125,126

2.4.2. Connecting Chemical Potential of Oxygen from 0 K to T K
Assuming the oxygen atom in the lattice site is in chemical equilibrium with the gas phase oxygen
molecule we can write Equation 2.6.
31

1

1

1

2

2

2

𝜇𝑂 (𝑇, 𝑝) = 𝜇𝑂2 (𝑇, 𝑝) = 𝐸𝑂𝐷𝐹𝑇
+ ∆𝜇𝑂2 (𝑇, 𝑝)
2

[2.6]

where ∆𝜇𝑂2 (𝑇, 𝑝) is the thermodynamic connection energy between DFT (0 K) calculated energy
of an isolated oxygen molecule with free energy of oxygen molecule at any temperature, T and
partial pressure, p. Equation 2.7 brakes ∆𝜇𝑂2 (𝑇, 𝑝) into three parts. The first part, ∆𝐻𝑂2 |𝑇0𝐾 is the
enthalpy of connection energy. All the enthalpy related terms in this derivation are highlighted in
red color. The second part, 𝑇∆𝑆𝑂2 |𝑇0𝐾 is related to entropy correction. All the entropy related terms
𝑝

are highlighted in green color. The third part, 𝑘𝑇𝑙𝑛 (𝑝 ) is the pressure correction term and
0

highlighted in blue color.
𝑝
∆𝜇𝑂2 (𝑇, 𝑝) = ∆𝐻𝑂2 |𝑇0𝐾 − 𝑇∆𝑆𝑂2 |𝑇0𝐾 + 𝑘𝑇𝑙𝑛 ( )
𝑝0

[2.7]

Equation 2.8 substitutes the enthalpy and entropy related terms into their integral form.
𝑇

𝑇

∆𝜇𝑂2 (𝑇, 𝑝0 ) = ∫ 𝐶𝑝 𝑑𝑇 − 𝑇 ∫
0𝐾

0𝐾

𝐶𝑝
𝑝
𝑑𝑇 + 𝑘𝑇𝑙𝑛 ( )
𝑇
𝑝0

[2.8]

The enthalpy contribution to the connection energy is obtained by integrating the heat capacity,
Cp over the temperature range 0 K to T (Equation 2.8) and it is separated into two parts as shown
in Equation 2.9 for the ease of calculation. Similarly, the entropy term of the connection energy
was also calculated in two parts as shown in Equation 2.10.
∆𝜇𝑂2 (𝑇, 𝑝0 ) =

𝑇𝑟 =298𝐾
∆ℎ𝑂2 |0𝐾

+

𝑇

𝑇

∍

𝐶𝑝 𝑑𝑇 − 𝑇 ∫

𝑇𝑟 =298𝐾

𝑇 =298𝐾

𝑟
Here, ∆ℎ𝑂2 |0𝐾

0𝐾

𝐶𝑝
𝑝
𝑑𝑇 + 𝑘𝑇𝑙𝑛 ( )
𝑇
𝑝0

[2.9]

is the connection enthalpy for oxygen molecule between 0 K and standard

state (298 K) which can be obtained in multiple ways as explained in next section.
𝑇𝑟 =298𝐾

𝑇

∆𝜇𝑂2 (𝑇, 𝑝0 ) =

𝑇𝑟 =298𝐾
∆ℎ𝑂2 |0𝐾

+

∍
𝑇𝑟 =298𝐾

𝐶𝑝 𝑑𝑇 − 𝑇

∍
0𝐾

32

𝐶𝑝
𝑑𝑇 − 𝑇
𝑇

𝑇

∍
𝑇𝑟 =298𝐾

𝐶𝑝
𝑝
𝑑𝑇 + 𝑘𝑇𝑙𝑛 ( )
𝑇
𝑝0

[2.10]

Rearrange the terms in Equation 2.10, it becomes
𝑇𝑟 =298𝐾

∆𝜇𝑂2 (𝑇, 𝑝) =

𝑇𝑟 =298𝐾
∆ℎ𝑂2 |0𝐾

−𝑇

∍
0𝐾

𝐶𝑝
𝑑𝑇 +
𝑇

𝑇

𝑇

∍

𝐶𝑝 𝑑𝑇 − 𝑇

𝑇𝑟 =298𝐾

∍
𝑇𝑟 =298𝐾

𝐶𝑝
𝑝
𝑑𝑇 + 𝑘𝑇𝑙𝑛 ( )
𝑇
𝑝0

[2.11]
𝑇 =298𝐾 𝐶𝑝

Subtracting and adding 𝑇𝑟 ∫0𝐾𝑟

𝑇

𝑑𝑇 in Equation 2.11 helps us to relate terms with chemical

potential of oxygen, as shown in later steps.
𝑇𝑟 =298𝐾

∆𝜇𝑂2 (𝑇, 𝑝) =

𝑇𝑟 =298𝐾
∆ℎ𝑂2 |0𝐾
𝑇

−𝑇

∍
𝑇𝑟 =298𝐾

𝑇 =298𝐾 𝐶𝑝

Now, ∫0𝐾𝑟

𝑇

− 𝑇𝑟

∍
0𝐾

𝐶𝑝
𝑑𝑇 − (𝑇 − 𝑇𝑟 )
𝑇

𝑇𝑟 =298𝐾

∍
0𝐾

𝐶𝑝
𝑑𝑇 +
𝑇

𝑇

∍

𝐶𝑝 𝑑𝑇

𝑇𝑟 =298𝐾

𝐶𝑝
𝑝
𝑑𝑇 + 𝑘𝑇𝑙𝑛 ( )
𝑇
𝑝0

[2.12]

𝑑𝑇 is the standard entropy and can be obtained from thermodynamic handbook

data and represented as 𝑆𝑟 in Equation 2.13.
𝑇𝑟 =298𝐾

∆𝜇𝑂2 (𝑇, 𝑝) =

𝑇𝑟 =298𝐾
∆ℎ𝑂2 |0𝐾

− 𝑇𝑟

0𝐾

𝑇

−𝑇

∍
𝑇𝑟 =298𝐾

∆𝜇𝑂2 (𝑇, 𝑝) =

− 𝑇𝑟

𝑇

∍

𝐶𝑝 𝑑𝑇

𝑇𝑟 =298𝐾

𝐶𝑝
𝑝
𝑑𝑇 + 𝑘𝑇𝑙𝑛 ( )
𝑇
𝑝0

𝑇𝑟 =298𝐾
𝑇𝑟 =298𝐾
∆ℎ𝑂2 |0𝐾

∍

𝐶𝑝
𝑑𝑇 − (𝑇 − 𝑇𝑟 ) 𝑆𝑟 +
𝑇

∍
0𝐾

𝐶𝑝
𝑑𝑇 +
𝑇

𝑝
− (𝑇 − 𝑇𝑟 ) 𝑆𝑟 + 𝑘𝑇𝑙𝑛 ( )
𝑝0

[2.13]

𝑇

𝑇

∍

𝐶𝑝 𝑑𝑇 − 𝑇

𝑇𝑟 =298𝐾

∍
𝑇𝑟 =298𝐾

𝐶𝑝
𝑑𝑇
𝑇

[2.14]
𝑇 =298𝐾

𝑟
Defining the enthalpy and entropy of connection between 0 K and 298 K (∆ℎ𝑂2 |0𝐾

𝑇 =298𝐾 𝐶𝑝

𝑇𝑟 ∫0𝐾𝑟

𝑇

−

𝑑𝑇) to be the change in chemical potential of oxygen ∆𝜇𝑂02 (𝑇𝑟 ) between 0 K and

standard state, leads to:
𝑇

𝑇

𝐶𝑝

𝑟

𝑟

𝑇

∆𝜇𝑂02 (𝑇) = ∆𝜇𝑂02 (𝑇𝑟 ) + ∫𝑇 =298𝐾 𝐶𝑝 𝑑𝑇 − 𝑇 ∫𝑇 =298𝐾
33

𝑑𝑇 − (𝑇 − 𝑇𝑟 ) 𝑆𝑟

[2.15]

The 0 K to 𝑇𝑟 connection energy for oxygen molecule, ∆𝜇𝑂02 (𝑇𝑟 ) can be calculated in
various ways. For instance, Reuter et al.124 used enthalpy and entropy data for oxygen at low
temperatures from a thermochemical handbook127 to estimate ∆𝜇𝑂02 (𝑇𝑟 ) as 𝑇 → 0 for their RuO2
surface calculations. In contrast, Lee et al.123 and Wang et al.92 computed ∆𝜇𝑂02 (𝑇𝑟 ) by utilizing
metal oxide formation reactions. They first took the connection energy of the reaction as the
difference between DFT calculated formation enthalpy of oxides and their corresponding enthalpy
of formation at standard state, then subtracted the connection energies of metal and oxides obtained
by integrating heat capacity and entropy data obtained from thermodynamics handbooks, in a
temperature range of 0 K to 𝑇𝑟 . We followed the method of Zhang et al.121 They constructed metal
oxide thermodynamic cycles to calculate ∆𝜇𝑂02 (𝑇𝑟 ) for various oxides and showed that it varies
with the forming elements since different metal-oxygen bonds are formed.121 Since Fe-O
polyhedra play the lead role in determining the oxygen vacancy concentration in the SrFeO3-δ
system, a ∆𝜇𝑂02 (𝑇𝑟 ) value of -0.56 eV for iron oxide from Zhang et al. was used to model SrFeO3δ.

The advantage of Zhang et al.’s metal oxide formation based approach121 is that it implicitly

accounts for the vibrational contribution of the oxygen atom in the lattice structure of the metal
oxide, and is therefore compatible with the use of Equation 2.4. Various methods exist in the
literature to calculate the oxygen connection energy and the chemical potential of oxygen.

34

Figure 2.4. Compare different methods to obtain oxygen connection energy in Equation 2.5.
Here, the computed oxygen connection energy, 0.5∆𝜇𝑂2 (𝑇, 𝑝 = 1𝑎𝑡𝑚), was compared in
a broad temperature range determined using the present model against other methods available in
the literature, as shown in Figure 2.4. Reuter et al.124 used the CRC handbook values for oxygen
connection energy to connect chemical potential of oxygen from 0 K to thermodynamic standard
state and they also considered the contribution of phonons. Lee et al.123 used an average oxide
formation energy shift for oxygen connection energy in their study. All three methods give the
same trend of ∆𝜇𝑂02 (𝑇, 𝑝) decreasing with temperature and similar ∆𝜇𝑂02 (𝑇, 𝑝) values (which
differ by less than 0.1 eV). This suggests that the free energy contribution due to the phonon
vibration of the extra oxygen atom in the perfect lattice is not large enough to alter the conclusions
𝑓

𝑓

presented in this paper. The 0.1eV difference in ∆𝐺𝑣𝑎𝑐 will impact more when ∆𝐺𝑣𝑎𝑐 is small. At
most, the predicted δ(T) relationship shown in Figure 6.4(b) could shift by ~200K at low δ region
and by ~100K at high δ region.

35

2.5.

Predicting Ionic Conductivity as a Function of Temperature and Phase Change in
Air

To calculate ionic conductivity, 𝜎 using Equation 1.1, it is required to compute oxygen vacancy
concentration and ionic mobility first, as Z and e are constant for oxygen vacancies.

2.5.1. Oxygen Vacancy Concentration in Bulk Lattice
Oxygen vacancy concentration is directly related to oxygen vacancy concentration as shown in
Equation 1.5. If the vacancy formation energy is independent of the oxygen non-stoichiometry (as
in case of materials with a dilute concentration of oxygen vacancies), the site fraction of oxygen
vacancies (X) at the SOFC operating conditions can be calculated via the expression:52,53
𝑋

= exp (−
1−𝑋

𝑓

∆𝐺𝑣𝑎𝑐 (𝑇,𝑝)
𝑅𝑇

)

[2.16]

𝑓

where R is the ideal gas constant and ∆𝐺𝑣𝑎𝑐 is the vacancy formation free energy calculated via
𝑓

Equation 2.4. However, if ∆𝐺𝑣𝑎𝑐 varies with X (as would be the case for interacting, non-dilute
oxygen vacancies), the right hand side of the Equation 2.16 has to be modified to include the
𝑓

𝑓

change in ∆𝐺𝑣𝑎𝑐 with X. In some cases (for example, cubic SrFeO3- ), the ∆𝐸𝑣𝑎𝑐 can change as
much as 1 eV with X. When this change is more dominant than the configurational entropy
contribution (which is ~0.1 eV at temperature > 600oC) due to ordered oxygen vacancies induced
by periodic boundary conditions, Equation 2.16 becomes:
𝑋

= exp (−
1−𝑋

𝑓

∆𝐺𝑣𝑎𝑐 (𝑇,𝑝,𝑋)
𝑅𝑇

)

[2.17]

(to account for the complex defect configurational entropy contribution, Equation 2.5 should be
modified accordingly following the method discussed by Xi et al.)128 Equation 2.17 needs to be
solved to find the correct operating temperature, T for a predetermined pressure, 𝑝𝑂2 with
36

calculated vacancy formation energy for a particular non-stoichiometry. Here, a numerical method
with an incremental search technique, presented in the Appendix, was developed to solve Equation
2.17 as 𝑋 was varied in each phase, T was varied from 0-2000 K and 𝑝𝑂2 was held constant at 0.21
atm.

2.5.2. Mobility and Diffusivity of Oxygen Vacancy in Bulk Lattice
Calculated oxygen ion migration barrier with CINEB method (described in Section 2.3.3)
will be used to calculate ionic diffusivity in solids. As the ultimate goal is to calculate ionic
conductivity using Equation 1.1, so for the calculation of diffusivity in Equation 1.8 only migration
barrier will be used in place of activation energy.129 Because in Equation 1.1 vacancy concentration
already incorporated the free energy of vacancy formation. The diffusivity prefactor D 0 can be
calculated from jump distance (d), jump frequency (ν), and jump correlation factor (f) using below
equation:130
1

𝐷𝑜 = 6 𝑓𝑑 2 𝜈

[2.18]

In a lattice with multiple vacancy migration pathways, d represent the average jump distance. The
jump frequency, ν was assumed as 1013 Hz following the solid state diffusion.131 Jump correlation
factor, f was assumed as 1.132 Once diffusivity is calculated at a certain temperature, the calculation
of ionic mobility is straight forward from Equation 1.7.

2.6.

Anisotropic Chemical Strain Due to Oxygen Vacancy
A 2x2x2 ceria cubic-supercell (Figure 2.5) was used for chemical strain calculation with a

single oxygen vacancy (VO•• ) formation (δ = 0.03125). Such a low δ was selected to avoid any
vacancy-vacancy interaction due to periodic boundary condition.27 In this 2x2x2 CeO2-δ supercell
37

the average distance between oxygen vacancies is 11.0 Å. Fully spin-polarized calculations with
generalized gradient approximation with Hubbard-U (GGA+U) was performed on a perfect and
oxygen vacancy containing (defective) lattices. Ueff = 4.5 was used to treat Ce 4f for highly
localized nature of f-orbital, following the rotationally invariant approach proposed by Dudarev et
al.90 and calculated for ceria by Farbis et al.133 This Ueff value for Ce 4f has provided satisfactory
results, as reported in previous literature.134,135 The partial density of states analyses for Ce atoms
adjacent to oxygen vacancy site can confirm the electron localization after the formation of neutral
oxygen vacancy. In earlier studies, it was observed that oxygen vacancy produces complicated
local deformation due to the charge disproportionation on Ce (4Ce4+ + 2e-  2Ce4+ + 2Ce3+).136,137

Figure 2.5. (a) CeO2 unit cell showing O-Ce4 tetrahedral coordination and Ce-O bonding
directions. (b) Charge distribution on Ce around neutral oxygen vacancy formation showing in a
2x2x2 supercell.

Esch et al. has shown experimentally, the charge localization on Ce (surface and subsurface) after VO•• formation.138 Nolan et al.139 observed charge localization on f-orbital forming
Ce3+ but their two Ce3+ were not equidistant from VO•• . Though the charge disproportionation was
observed in previous studies, they didn’t use the local lattice deformation to calculate the chemical
38

strain tensor for the defect containing structure.139 In stoichiometric CeO2 fluorite structure, CeO8 forms cubic coordination but O-Ce4 forms a tetrahedral coordination (Figure 2.5(b)). As each
O is connected with the four Ce4+ atoms, on formation of VO•• , two electrons should be distributed
to the four adjacent Ce atoms but due to the presence of highly localized 4f orbital (as discussed
above) two electrons are localized to two Ce-atoms converting them to Ce3+, whereas, other two
remained Ce4+. The charge localization is shown in Equation 2.19 with KrĂśger-Vink notation.
1

𝑥
′
2𝐶𝑒𝐶𝑒
+ 𝑂𝑂𝑥 → 2 𝑂2(𝑔) + 𝑉𝑂•• + 2𝐶𝑒𝐶𝑒

[2.19]

Following the computational study of Tang et al.140 and Castleton et al.,141 calculations were
performed with two different scenarios. In the first case, atom positions were relaxed after VO••
formation keeping the cell volume constant, where four Ce atoms equally moved away from the
VO•• site producing uniform distortion. In the second case, after VO•• formation two Ce atoms were
perturbed then atom positions were relaxed keeping the cell volume constant. In second case two
Ce atoms moved further than other two, causing asymmetric distribution of Ce atoms around VO•• .
The anisotropic chemical strain coefficient tensor was calculated following the work of
James et al.42 The chemical strain tensor was calculated by minimizing the total energy, ∆𝐸𝑡𝑜𝑡𝑎𝑙
which consists of energy due to local deformation from oxygen vacancy formation, ∆𝐸𝑠ℎ𝑜𝑟𝑡 and
the energy due to long-range elastic strain, ∆𝐸𝑙𝑜𝑛𝑔 . Short range energy change, ∆𝐸𝑠ℎ𝑜𝑟𝑡 was
originated due to local strain dipole, 𝐺 from the VO•• . The strain dipole can be calculated, taking
𝑓

𝑓

𝑓

first order Taylor expansion of VO•• formation energy with strain applied, 𝐸V•• ,𝜀 . 𝐸V•• ,𝜀 = 𝐸V•• ,𝜀=0 +
O

O

O

𝑓

𝐺: 𝜀, where 𝐸V•• ,𝜀=0 is the VO•• formation energy without any applied strain. Taking derivative of
O

𝑓

the preceding equation with respect to 𝜀, we obtain 𝐺 = 𝑑𝐸V•• ,𝜀 ⁄𝑑𝜀 . Energy for the short and long
O

𝛿

1

rangee interaction can be obtained as ∆𝐸𝑠ℎ𝑜𝑟𝑡 = 𝑉 𝐺: 𝜀, and ∆𝐸𝑙𝑜𝑛𝑔 = 2 (ℂ: 𝜀): 𝜀 respectively.
𝑈

Where, 𝑉𝑈 is the volume per formula unit of the perfect CeO2 and ℂ is the elastic stiffness tensor
of the perfect lattice. Chemical strain tensor, 𝜀𝐶 and chemical expansion coefficient tensor, αC can

39

be obtained by minimizing the total energy with respect to 𝜀𝐶 from Equation 2.20. The derivation
is provided below:
𝛿

1

∆𝐸𝑡𝑜𝑡𝑎𝑙 = ∆𝐸𝑠ℎ𝑜𝑟𝑡 + ∆𝐸𝑙𝑜𝑛𝑔 = 𝑉 𝐺: 𝜀 + 2 (ℂ: 𝜀): 𝜀

[2.20]

𝑈

𝑑∆𝐸𝑡𝑜𝑡𝑎𝑙

=

𝑑𝜺
𝑑∆𝐸𝑠ℎ𝑜𝑟𝑡
𝑑𝜺
𝑑∆𝐸𝑙𝑜𝑛𝑔
𝑑𝜺
𝑑∆𝐸𝑙𝑜𝑛𝑔
𝑑𝜺

=

1 𝑑
2 𝑑𝜺
1

𝑑∆𝐸𝑠ℎ𝑜𝑟𝑡
𝑑𝜺
𝑑

+

𝑑∆𝐸𝑙𝑜𝑛𝑔

=0

𝑑𝜺

𝛿

[2.21]

𝛿

= 𝑑𝜺 (𝑉 𝑮: 𝜺) = 𝑉 𝑮
𝑈

[2.22]

𝑈

1 𝑑(ℂ:𝜺)

{(ℂ: 𝜺): 𝜺} = {
2
𝑑𝜺

𝑑𝜺

1

𝑑𝜺

1

𝑑𝜺

} : 𝜺 + 2 (ℂ: 𝜺): 𝑑𝜺
1

1

= 2 {ℂ: 𝑑𝜺} : 𝜺 + 2 (ℂ: 𝜺): 𝑑𝜺 = 2 (ℂ: 𝜺) + 2 (ℂ: 𝜺) = ℂ: 𝜺
𝑑∆𝐸𝑡𝑜𝑡𝑎𝑙
𝑑𝜺

𝛿

= 𝑉 𝑮 + ℂ: 𝜺 = 𝟎
𝑈

𝜀𝐶 =
𝛼𝐶 =

−𝛿 (ℂ−1 𝐺)

𝛿

=

[2.24]
[2.25]
[2.26]

𝑉𝑈
𝜀𝐶

[2.23]

−(ℂ−1 𝐺)
𝑉𝑈

40

[2.27]

3.

DFT+U as a Reliable Tool to Predict Lattice Structure and Electronic Properties of
MIEC at Low Temperature

This Section studies the perfect LSF lattices with DFT+U (0 K) calculations. Here, we have
summarized properties like magnetic orientation, band structure, lattice parameters etc from 0 K
calculations. The predicted properties can be compared with low-temperature experimental
observations.
Table 3.1. DFT predictions of stable LSF structures under ground-state conditions.

The bold phase data under each column signify the stable structure with the preferred magnetic
orientation.
The tilt angle column denotes the average octahedral tilt among Fe-O octahedra in the lattices.
* The last row contains experimental data for comparison.
Figure 2.1 shows the energetically favored crystal structures for LFO, SFO, and LSF55. In
order to determine the ground state lattice structure for each LSF composition, three possible
crystal structures each with two possible Fe magnetic ordering arrangements (ferromagnetic, FM,
and anti-ferromagnetic, AFMG), were computed. From energy minimization, the stable, lowestenergy structures, presented in Table 3.1, were determined and illustrated in Figure 2.1. In Table
3.1, ΔE denotes the excess energy of a structure above the minimum energy of each composition,
the ‘μB’ column denotes the magnetic moment on Fe in the structure, and the “tilt angle” denotes
the Fe-O6 octahedral tilt angle (which varies with the lattice structure and La/Sr ratio).

41

3.1

Predicting Ground State Lattice Structure and Magnetic Ordering in La1-xSrxFeO3δ Phases

In accordance with experiments,100 the orthorhombic (Pbnm) crystal structure with AFMG
ordering was predicted as the stable LFO structure. The calculated lattice of 5.672 Å x 7.923 Å x
5.582 Å also compared well to the experimentally reported lattice parameters.86,100 The Fe-O6
octahedra in LFO showed a Z-out distortion (the b-axis in Figure 2.1(a)) with four 2.040 Å Fe-O
bonds in the XY plane (the ac plane in Figure 2.1(a)) and two 2.057 Å bonds in the Z-direction (baxis in Figure 2.1(a)). These Fe-O bond lengths are comparable to the experimentally measured
average bond length of 2.006 Å.86 The computed octahedral tilt was 154o (the bond angle of FeO-Fe) which was comparable to experimentally reported octahedral tilt of 157o.86
For stoichiometric SFO, the predicted phase is the experimentally reported lowtemperature cubic structure (space group: 𝑃𝑚3̅𝑚) with FM Fe ordering and a lattice parameter of
3.881 Å, which agreed well with experiments32 (note, FM ordering was used to simplify the
experimentally-reported helical Fe magnetic orientation arrangement in the 2K neutron diffraction
SFO study of Reehuis et al.).105 The predicted magnetic moment on Fe was 3.61 ÎźB, which is
similar to a literature DFT+U value of 3.7 ÎźB calculated with U=5.4.101 All the Fe-O bond lengths
in SFO were equal to 1.940 Å due to octahedral (Oh) symmetry (Figure 2.1(d)) and were
comparable to the experimentally reported bond length of 1.928 Å (measured at room
temperature).145 A detailed study on other SFO phases are provided in Section 3.3.
For LSF55, the cubic, orthorhombic, rhombohedral structures were computed. The unit
cell for rhombohedral LSF55 was depicted by a hexagonal lattice following the work of Kharton
et al.46 In light of the 0 to ~600K single phase behavior of LSF55 observed by X-ray
diffraction65,143,145 and the insensitive oxygen vacancy formation energy on the site ordering of
La/Sr simulated by Ritzmann et al.,24 a La/Sr ordered LSF55 phase was assumed. Thus, the 6a
42

Wyckoff positions were assumed to be equally shared by La and Sr and an alternate layer of LaO and Sr-O was assumed, separated by the Fe-containing layer, as shown in Figure 3.1(b). The
present calculations correctly predicted that a hexagonal (rhombohedral) lattice structure with an
R32 space group is energetically favorable for LSF55. The calculated lattice length of the
hexagonal unit cell was a = 5.564 Å and c = 13.451 Å, and compared well with the experimental
values (a = 5.5150 Å, c = 13.4287 Å).46 In any octahedron, three of the calculated Fe-O bond
lengths were 1.986 Å and the other three were 1.971 Å. Along the Fe-O-Fe-O direction, the Fe-OFe bond angle between two Fe-O6 octahedra alternated between 166˚ and 161˚ and the Fe-O bond
lengths also alternated between 1.986 and 1.971 Å. These calculated bond lengths are comparable
to experimentally reported the average Fe-O bond length of 1.954 Å (at room temperature).145
Although it was predicted that FM order had a lower energy than AFM, the energy difference was
only ~0.06 eV/f.u., and hence not in disagreement with the AFM magnetic orientation measured
in LSF55 with MĂśssbauer spectroscopy at 78 K.98
3.2.

Density of States for LaFeO3-δ and La0.5Sr0.5FeO3-δ

Figure 3.1 shows the band gap for each energetically favored lattice structure. The partial density
of states (PDOS) was calculated for both FM and AFMG ordering in LFO as shown in Figures
3.1(a) and 3.1(b). The computed band gap was 0.5 and 2.0 eV for FM and AFMG ordering,
respectively. The band gap predicted in AFMG LFO was comparable to the experimentally
observed optical band gap of 2.1 eV measured at room temperature.146 Cubic SFO shows a metallic
behavior,105 as correctly predicted in Figure 3.2(a). The detailed study on measured properties for
the different SrFeO3- phases are provided in Section 3.3. The PDOS was calculated for both FM
and AFMG magnetic ordering in LSF55 and reported in Figure 3.1(c) and 3.1(d). Though there

43

are few DOS observed close to Fermi energy with an FM orientation (Figure 3.1(c)), the AFMG
ordered structure has an energy gap of 1.6 eV, comparable to earlier computational studies.24
Overall, the DFT settings successfully predicted that orthorhombic AFM, rhombohedral
(hexagonal) (FM or AFM), and cubic (FM) were stable structures for LFO, LSF55, and SFO,
respectively, at low temperature.

Figure 3.1. GGA+U (=3) partial density of states calculations for (a-b) orthorhombic LaFeO3,
(c-d) rhombohedral La0.5Sr0.5FeO3. The gray and black lines represent the Fe 3d and O 2p states
respectively.

44

3.3.

Study of Different SrFeO3-δ Phases

Table 3.1. Comparison of DFT-calculated and experimentally-measured lattice parameters,
magnetic ordering, and band gaps for several phases of SrFeO3-δ (“See DOS” in the last
column of below-provided table refer to the Figure 3.1 for the density of states plot). N.C.
stands for Not Calculated. Experimental NĂŠel Temperature, TN, is also listed for each
phase.

* Calculated from neutron powder diffraction data
Calculated from powder X-ray diffraction data

#

Figure 1.1 shows the calculated perfect crystal structures (at X=0), the charge state of Fe, and the
various oxygen sites for each of the perfect SrFeO3-δ phases observed from 0-1000oC in the air:
cubic SrFeO3, tetragonal SrFeO2.875, orthorhombic SrFeO2.75, and brownmillerite SrFeO2.5. Table
3.1 lists the computed lattice parameter, magnetic ordering, and band gap information on these
phases with the DFT settings discussed in Section 2.1.1. As discussed, Ueff = 3 led to the best
agreement of magnetic moment and lattice parameters with the experimental values for both
LaFeO3 and SrFeO3 bulk structures compared to other U-parameters over the range of 0.2 to 4.0.

45

The computed lattice parameters for all four phases of SrFeO3-δ were also in good agreement with
experimental results.

Figure 3.2. GGA+U (Ueff = 3) PDOS for four different SrFeO3-δ phases. The gray and black line
represent Fe 3d and O 2p states respectively. The stoichiometric SrFeO3 has shown metallic
behavior. Whereas, band gaps are observed in oxygen deficient phases.

46

3.3.1

Magnetic Orientations of SrFeO3-δ Phases

SrFeO3-δ displays complicated magnetic ordering. Reehuis et al.105 detected incommensurate
magnetic structure in cubic SrFeO3 (helical) and SrFeO2.87 with neutron diffraction at 2 K. They
also reported commensurate G-type antiferromagnetic (AFMG) ordering for SrFeO2.75 in the same
study. Greaves et al.150 reported AFMG ordering on brownmillerite SrFeO2.5 by analyzing neutron
powder diffraction data collected at 4.2 K. For cubic SrFeO3, 0 K DFT+U (Ueff = 3) calculated
energies for different AFM (G, A, C) and ferromagnetic (FM) configurations indicate that FM has
the lowest energy (AFMG and FM energies are compared in Table 3.1). Due to computational
complexity, FM ordering for Fe in cubic SrFeO3 was chosen instead of helical ordering. Like the
helical magnetic structure, FM ordering also maintains the metallicity of cubic SrFeO3, as later
shown with the density of states calculations. Further, these low temperature (0 K) calculations
indicate that while the FM structure has lower energy in cubic and tetragonal phases, orthorhombic
and brownmillerite phases show lower energy in AFMG structures.105,150 In addition, the
calculations correctly predict the experimentally-observed AFM structures in SrFeO2.75 and
SrFeO2.5. Although Fe ions in SrFeO3-δ phases display AFM orientation at low temperature,151 the
energy difference is small (less than 0.26 eV). Furthermore, these materials are well above their
Neel temperatures at the SOFC operating temperatures of interest.105,152,153 Therefore, in the
present work it was not considered necessary to account for antiferromagnetic ordering123 when
performing oxygen vacancy formation computations at SOFC operating temperatures, and only
ferromagnetic ordering was utilized. As orthorhombic and brownmillerite phases prefer AFMG
orientation at low temperature, we have calculated oxygen vacancy formation in these two phases.
Our calculations showed that vacancy formation energies in AFMG orthorhombic and
𝑓

brownmillerite phases were ~0.4 and ~0.2 eV higher than their ∆𝐸𝑣𝑎𝑐 (1.3 and 3.0 eV respectively)

47

in FM phases. Further converting these oxygen vacancy formation energy to Gibbs free energy
using Equation 4 and relating corresponding oxygen nonstoichiometry with temperature applying
Equation 7 showed that AFMG orthorhombic and brownmillerite phases overestimated the
required temperature for the oxygen nonstoichiometry by 280oC and 140oC respectively at 𝑝𝑂2 =
0.21 atm. Therefore, it was reasonable to perform oxygen vacancy calculations on FM SrFeO3-δ
phases.

3.3.2. Density of States for SrFeO3-δ Phases
It is known that GGA calculations underestimate the band gap in materials, and GGA+U
sometimes can predict band gaps closer to experimental values.154 Figure 3.2 shows the GGA+U
(Ueff = 3) predicted partial density of states (PDOS) plots of Fe 3d and O 2p orbitals in the phases
of SrFeO3-δ. In order to compare with low-temperature experimental results, the PDOS for the four
low-temperature phases were first plotted in Figure 3.2 (a-d). Figure 3.2(a) shows that no band gap
exists in the bulk cubic FM SrFeO3 structure, consistent with its observed metallic behavior.105
Figure 3.2(d) shows Brownmillerite AFMG SrFeO2.5 has a clear band gap of 1.7 eV, which is
comparable to the experimental band gap of 1.5 Âą 0.2 eV.147 Figure 3.2(b) shows that tetragonal
FM SrFeO2.875 has no distinct band gap. However, the density of states is relatively low within the
region of the Fermi level to 0.3 eV above, which may be related to the observed small band gap in
this phase. Figure 3.2(c) shows orthorhombic AFMG SrFeO2.75 has a band gap of ~1.7 eV but
small polaron states exist in the band gap region. No experiments on the exact composition of these
two phases have been reported. However, experimental bandgap values of 1.5 Âą 0.2 and 1.0 Âą 0.2
eV were observed for non-stoichiometric SrFeO3-δ i.e., SrFeO2.52 and SrFeO2.82 respectively.147
GGA+U (Ueff = 3) gives reasonable values of band gaps in SrFeO3-δ, especially since it captures
the metallic to insulator transition from cubic SrFeO3 to Brownmillerite SrFeO2.5. Since the oxygen
48

vacancy concentrations will be computed using FM phases, Figure 3.2(e) and 3.2(f) present the
PDOS for these two structures. These figures show that a low DOS exists around and below the
Fermi level in orthorhombic FM SrFeO2.75 and Brownmillerite FM SrFeO2.5, consistent with the
high electronic conductivity in these phases at SOFC operating conditions.
3.3.3. Formation Enthalpy of Brownmillerite SrFeO2.5
Before calculating the oxygen vacancy formation energy, the enthalpy of formation for
brownmillerite (FM and AFMG both) phases were calculated and compared with experimental
data to further calibrate the choice of DFT+U parameters. The formation enthalpy of
brownmillerite SrFeO2.5 from its constituent simple oxides was calculated as 0.58 eV and 0.86 eV
for FM and AMFG phases, respectively using Equation 3.1:
1

𝑆𝑟𝑂 + 2 𝐹𝑒2 𝑂3 = 𝑆𝑟𝐹𝑒𝑂2.5

[3.1]

The enthalpy of formation of the FM brownmillerite phase is comparable with previous
computational results obtained using the GGA method (0.54 eV)155 and with experimental
estimation at room temperature (0.48 Âą 0.02 eV).156

3.3.4. d-Orbital Splitting and Its Impact on the Iron Oxidation State
Outside a crystal lattice, the half-filled outer d-orbital of Fe3+ (with a [Ar] 3d5 electron
configuration) makes it energetically more stable than Fe4+ (with a [Ar] 3d4 electron configuration).
However, in strontium ferrite, crystal field Fe-O d-orbital splitting dramatically alters the local
charge state of iron. As clearly demonstrated in the calculations here and determined from the
literature of Fe magnetic moment experiments,32,33 the outer Fe 3d electrons in strontium ferrite
adopt a high spin configuration. Thus, Fe4+ in SP coordination has eg2b2g1a1g1 configuration that

49

requires significant amounts of energy (based on the high energy of the SP b1g site in Figure 3.3(a))
to make it a Fe3+ with an eg2b2g1a1g1b1g1 electron configuration. In contrast, Figure 3.3(b) and past
experiments105 show that the Oh coordinated Fe4+ in cubic SrFeO3 with an electronic configuration
of t2g3 eg1 requires less energy to become Fe3+ (based on the lower energy of the Figure 3.3(b) eg
orbital compared to the Figure 3.3(a) b1g orbital). Even in distorted Oh of tetragonal strontium
ferrite (as shown in Figure 3.3(c)), Fe4+ with a b2g1eg2b1g1 electron configuration that can relatively
easily be reduced compared to the Fe4+ in SP. This is the underlying reasons for the various Fe
charge distributions found within different SrFeO3-δ phases.

Figure 3.3. Crystal Field Theory explains d-orbital splitting and electronic distribution for
different Fe-O polyhedra. (a) Square pyramidal Fe-O coordination, (b) octahedral Fe-O
coordination (c) Compressed octahedral Fe-O coordination with Jahn-Teller distortion.
3.3.5. Solving the Debate on Fe Oxidation State at Square Pyramidal and Octahedral
Sites
The Fe charge distribution reported here also resolves the debate regarding the charge state of Fe
ions in strontium ferrite phases containing mixed iron oxidation states, such as tetragonal
SrFeO2.875 and orthorhombic SrFeO2.75. In particular, there has been much debate over whether SP
50

coordinated iron is Fe4+

32,35

or Fe3+.

33,34,36,60

For instance, Hodges et al.32 concluded from their

MĂśssbauer study and empirical bond strength analysis that Fe4+ and Fe3+ ions occupy the SP and
Oh sites respectively. Adler35 produced similar conclusions based on his MĂśssbauer analysis. On
the contrary, Schmidt et al.33 assigned the Fe inside the SP to be Fe+3 and the Fe inside Oh to be
Fe+4, based on the magnetic structure obtained from MĂśssbauer measurements. Unfortunately,
Schmidt et al.’s conclusion may be due to the incorrect assumption that the observed Fe-O
octahedra compression resulted from Fe+4 induced Jahn-Teller distortions.33 This assumption is at
odds with 1) literature MĂśssbauer studies on cubic SrFeO3 showing that octahedrally coordinated
Fe+4 has no Jahn-Teller distortion,142,157 and 2) results here indicating there is no Jahn-Teller
distortion for Fe-O Oh containing Fe+4 in either SrFeO3 or SrFeO2.875. In fact, Jahn-Teller
distortions are only observed for octahedra containing Fe with oxidation states ranging from +3 to
+3.6 based on the present work. In contrast to the results reported here, DFT calculations by Vidya
et al.34 and Ravindran et al.36,60 have also claimed SP iron to be Fe+3. However, those modeling
results are likely to be incorrect since those studies did not apply U corrections to their GGA-based
simulations, and U parameters are generally needed to correct the orbital energies of partially
occupied localized d-states in Fe-O polyhedra.91
3.4.

Choices of Crystal Structures at SOFC Operating Temperature

Since high temperature oxygen vacancy nonstoichiometry is important for SOFC applications,
oxygen vacancy formation was studied in the high-temperature LSF phases experimentally
determined to be stable.65,158 In terms of high-temperature crystal structure, although LFO
transforms from orthorhombic to rhombohedral at 1278Âą5 K (in air), rhombohedral LFO was not
considered since the phase transition temperature is above typical SOFC operating temperatures.
On the other hand, LSF55 shows a second order phase transition from rhombohedral to hightemperature cubic at around 523Âą50 K (in air).65 For this reason, oxygen vacancy formation
51

calculations were performed in both cubic and rhombohedral LSF55 phases. The rhombohedral
LSF55 structure is shown in Figure 2.1(c). SFO undergoes multiple vacancy ordering induced
phase transitions as it progresses from cubic to tetragonal to orthorhombic to brownmillerite in a
temperature range of 273 to 1273 K.44 The oxygen vacancy formation calculations in these four
SFO phases were performed by Das et al.,27 and some of those computational results are presented
in this paper for comparison. As shown in Figure 2.1, in all these structures, the perovskite (ABO3)
B-site (transition metal, here Fe) was coordinated to six oxygen atoms, forming an octahedral (Oh)
Fe-O6 coordination.159 All the oxygen atoms were shared by two octahedra at their corners. After
oxygen vacancy formation, the octahedral Fe adjacent to oxygen vacancy formation site went into
Fe-O5 square-pyramidal (SP) coordination.
In term of magnetic ordering, the NĂŠel temperature for LFO is 750 K100, and it is rapidly
reduced with increasing Sr content. In fact, for LSF55 it is ~210 K,87 and for SFO it is 133 K.105
These NĂŠel temperatures are much below the temperature range used for SOFC operation (and
hence the oxygen vacancy site fraction calculations). Thus, Fe in all the LSF phases lose their
magnetic ordering and become paramagnetic (PM) at the SOFC operating temperature (>1000 K).
However, due to the computational complexity of modeling PM ordering, FM ordering was
assumed on Fe,27,55 for all oxygen vacancy formation calculations.

52

4.

Oxygen Vacancy Interaction via Vacancy Polaron
This study first time explains the oxygen vacancy interaction in compounds of LSF family

via oxygen vacancy polaron. The oxygen vacancy polaron plays a significant role in controlling
oxygen vacancy concentration due to vacancy interaction at a higher concentration as discussed in
Chapter 5.
4.1.Origin of Long Range Charge Transfer due to Oxygen Vacancy Polaron

Figure 4.1. The charge distribution on Fe due to oxygen vacancy formation is described on
vacancy containing lattice plane of a 4x4x4 cubic SrFeO3 supercell with a single oxygen vacancy
at the center of the cell. Figure (a) is the lattice plane without relaxing the atom position (only
electron cloud is allowed to relax after vacancy formation). Figure (b) is the lattice plane
containing atoms after relaxing the atom positions. The light brown atoms are Fe and red atoms
are O. The light blue Sr atoms are above or below the described planes. The oxygen vacancy is
denoted with a black dot. Figure (a) and (b) contains the oxygen vacancy at the center of the (200)
lattice plane. The charge on every Fe atom, rounded to the nearest tenth, is listed next to each Fe
atom.

The strong long-range interactions between oxygen vacancies in cubic SrFeO3explained by the unique charge transfer mechanism that occurs in this material, which is also due
to the different Fe-O d-orbital splitting in Oh and SP. For instance, Figure 4.1 shows how oxygen

53

vacancy formation affects the charge of neighboring Fe in cubic SrFeO2.98. Here, one oxygen
vacancy is created at the center of a 4x4x4 cubic SrFeO3 supercell, in which the distance between
the oxygen vacancies is 15.52 Å (i.e. larger than the interaction distance of 11 Å identified in
Section 5.2). In this (200) oriented cross-section, the oxygen vacancies and the center of the OhOh octahedra sit on the plane, while the Sr atoms sit above and below the plane. Figure 4.1(a)
shows the unrelaxed crystal structure after an oxygen vacancy is created. Here, the Fe sits on the
basal plane of the SP polyhedra, and oxygen vacancy formation has reduced the charge on the two
iron atoms bordering the oxygen vacancy from 4+ to 3.7+. In contrast, Figure 4.1(b) shows that
after atomic relaxation, the Fe relaxes into the SP structure by moving 0.10 Å from the SP basal
plane, and the O in the SP polyhedra have moved toward the SP Fe. Together, these changes
decrease Fe-O bond lengths in the SP polyhedra and increase those in the neighboring Oh
polyhedra. Further, after the atomic relaxation, the two iron atoms in SP coordination become
Fe4.2+. Due to the different Fe-O d-orbital splitting in Oh and SP, the SP iron adjacent to a newly
created oxygen vacancy in a cubic or tetragonal SrFeO3-δ structure (not shown) remains 4+ by
maintaining its unpaired electron and passes the electrons resulting from oxygen loss to nearby
octahedrally coordinated iron instead of reducing its own oxidation state. The transfer of electrons
to the second nearest neighboring Fe atoms, instead of localization on the Fe atoms connected to
the oxygen vacancy (i.e. transfer to the first nearest neighbor Fe atom) is a previously unknown
long-range charge transfer mechanism. This comparison also indicates that structural relaxations
are a necessary step for observing Fe4+ in SP polyhedra. The charge transfer discussed here was
calculated at 0 K, although electrons will be more delocalized, as this material remains
paramagnetic at SOFC operating temperature, it is anticipated that the d-orbital splitting for Oh
and SP coordination will qualitatively remain same (as shown in Figure 3.3) with increasing

54

temperature and/or a with change in oxygen partial pressure; only the gap between energy levels
will change which is not going to impact the overall electronic distribution. Thus, the long-range
charge transfer phenomena originated by different d-orbital splitting in Oh and SP will still present.
Surrounding each oxygen vacancy, the eight Fe atoms connected to the SP polyhedra via the SP
basal plane oxygen are reduced from 4+ to 3.8+. Therefore, one can think of the vacancy polaron
as having a “pan-cake” shape that encloses these reduced Fe atoms, where the purple dotted line
in Figure 4.1(b) represents the cross-section of the pancake. As shown in Figure 4.1(b), the distance
from the oxygen vacancy to the reduced Fe is ~4.3 Å. Due to the reduced charge on these Fe,
vacancy formation around them becomes more difficult. This results in the increased oxygen
vacancy formation energy when the distance between the vacancies is less than twice of 4.3 Å (i.e.
when the vacancy polarons overlap). Thus, the long-range charge transfer on Fe is responsible for
the strong vacancy interaction (indicated by the rapid increase in ΔEδ when δ > 0.1 in Figure 5.2)
occurring when the distance among oxygen vacancies is within 8 Å, as is true for all SrFeO 3polaron, one would
expect that they may prefer to be aligned along the short-axis of the pancake, which will be subject
to future studies, and may illuminate the oxygen vacancy ordering and phase transformation
𝑓

mechanisms. This may lead to some low energy configurations, so ∆𝐸𝑣𝑎𝑐 for cubic SrFeO3- at
𝑓

lowered, but the increasing trend of ∆𝐸𝑣𝑎𝑐
long-range charge transfer induced repulsion interaction.
4.2.Oxygen Vacancy Polaron in Different La1-xSrxFeO3-δ Phases
In this section, the vacancy formation sites and the oxidation state changes on Fe upon oxygen
vacancy generation in various LSF phases were analyzed to determine the shape and size of the

55

oxygen polarons shown in Figure 4.2. The Fe oxidation state changes on were explained by the
electron occupancy and crystal-field-theory schematics shown in Figure 4.3.

Figure 4.2. Long-range charge transfer areas (denoted by the purple dashed lines in a-d) and
perspective views of the LaxSr1–xFeO3 lattice polarons (e-h) for (a, e) orthorhombic LaFeO3, (b,f)
rhombohedral La0.5Sr0.5FeO3, (c,g) cubic La0.5Sr0.5FeO3, and (d,h) cubic SrFeO3. Note, in (a-d)
different Miller planes were used for each phase to compare lattice orientation with same visual
conformity. In (e-h) the polaron shape was highlighted by only drawing sides to those Fe-O
polyhedra with a Fe atom that received excess electrons due to oxygen vacancy formation. In (ad), the black lines indicate the border of the supercells used for computations, while the arrows
provide an estimate of the polaron size by denoting the distance the distance from the outermost
charge-altered Fe atoms to the oxygen vacancy site.
4.2.1. LaFeO3-δ
In LFO, there are two different oxygen vacancy formation sites, one of which is between the long
(2.057 Å) and short (2.040 Å) Fe-O bonds while the other one is between two short (2.040 Å) FeO bonds (Figure 2.1(a)). The oxygen vacancy formation calculations showed that the two-different
oxygen vacancy formation sites have similar oxygen vacancy formation energies, differing by less
than 0.01eV (i.e. within the energy resolution of the present study). The magnetic moment

56

interpreted charge discussed in Section 2.1, showed that after oxygen vacancy formation, the two
Fe atoms connected to the oxygen vacancy site changed their charge states from 3+ to 2.3+, while
all the other Fe atoms in the lattice maintained their charge at 3+ (as in perfect LFO).

Figure 4.3. Schematic of the d-orbital splitting in the Fe in different polyhedra showing that the
excess electrons generated during oxygen vacancy formation (in red) will be localized to oxygen
vacancy adjacent square pyramidal Fe in orthorhombic LaFeO3 (a), and transferred to second
nearest neighboring octahedral Fe in cubic SrFeO3 (b).
This indicates that after oxygen vacancy formation, the excess electrons generated by
neutral oxygen vacancy formation remained localized to the oxygen vacancy adjacent area and
formed a relatively small polaron, as shown in Figure 4.2(a). This charge localization can be
simply explained by the electronic distribution on the oxygen-vacancy-adjacent Fe d-orbitals.
According to crystal field theory, the Fe d-orbitals in the octahedral coordination found in the
perfect LFO lattice split into t2g3eg2, which can be further split into b2g1eg2b1g1a1g1 due to JahnTeller distortions, as shown in Figure 4.3(a). All the Fe atoms in perfect LFO are in the 3+ highspin state with each Fe containing five electrons in its 3d orbital ([Ar] 3d5), indicated by the blue
arrows in Figure 4.3(a). As a geometric consequence of losing a neighboring oxygen during
oxygen vacancy formation, the two first-nearest-neighboring Fe atoms become SP coordinated,
creating a different d orbital splitting situation of eg2b2g1a1g1b1g1. To minimize the system energy,
the two electrons produced during neutral oxygen vacancy formation are distributed to the lowest
57

energy level of d-orbitals of each Fe in the SP coordination. Hence, in LFO the excess (i.e. polaron)
electrons remain localized and adjacent to the oxygen vacancy site. A schematic of the polaron
(shown in yellow in Figure 4.2(a) and with a perspective view of oxygen vacancy adjacent lattice
area in Figure 4.2(e)), show the resulting small LaFeO3-δ polaron which is elongated along the Z
direction with a length of ~4.0 Å.

4.2.2. SrFeO3-δ
In cubic SFO, as discussed in Section 4.1 the electrons produced by a neutral oxygen vacancy
formation do not remain localized on the first-nearest-neighboring Fe connected to the oxygen
vacancy. The Fe in perfect SFO has an [Ar] 3d4 4+ charge state, i.e., 4 electrons in the 3d orbitals,
as indicated by the blue arrows in Figure 4.3(b). After oxygen vacancy formation, the excess
electrons generated due to oxygen vacancy formation (i.e. the red arrows in Figure 4.3(b)) transfer
to a1g level in the Oh-Fe atoms instead of a higher energy b1g of the nearest-neighbor, SPcoordinated Fe ions. Further, since there are eight second-nearest neighboring Oh-Fe present with
the same a1g energy level around each oxygen vacancy site, the excess electrons are distributed
between these eight Fe atoms to produce the “pancake” shaped polaron shown in Figures 4.1(d)
and 4.1(h) with an estimated polaron diameter and height of ~7.8 and ~3.9 Å, respectively.
Figure 4.4(a) shows that the oxygen vacancy polaronic distribution in tetragonal SFO
structure resembles the polaron shape in cubic SFO (Figure 4.2(h)). Here, the formation of an
oxygen vacancy pushes the excess electrons to the second nearest neighboring octahedral (Oh) Fe
rather than localizing it to the oxygen vacancy adjacent square-pyramidal (SP) Fe. The charge on
the oxygen-vacancy-adjacent SP-Fe remains ~4+ (a result obtained from a magnetic moment
analysis) before and after oxygen vacancy formation. In contrast, the charge on the second nearest
neighboring Oh-Fe atoms is reduced from 3.6+ to 3.3+ due to a gain in partial electrons. Due to
58

symmetry, the two electrons are shared between the eight second nearest neighbor Oh-Fe atoms
and the shape of the polaron is denoted by the dashed purple line in Figure 4.4(a).

Unlike tetragonal SFO, the excess electrons left behind by neutral oxygen vacancy
formation remain localized to the oxygen-vacancy-adjacent Fe atoms in orthorhombic and
brownmillerite SFO. This occurs because, in orthorhombic SFO, oxygen vacancy generation
occurs next to Fe4+ ions originally in square pyramidal coordination and transform to a square
planar (SP) coordination with oxygen vacancy formation (as shown in Figure 4.4(b)). Since a
Figure 4.5(a) SP dx2-y2 electron would have a relatively high-energy state, after oxygen vacancy
formation the five electrons in SP- Fe3+ adopt a low spin configuration after oxygen vacancy
formation (with one electron each in the z2 and xy states and three electrons shared between the yz
and zx states). The orthorhombic to brownmillerite phase transition is brought about by the need
to accommodate the electrons generated by oxygen vacancy formation. Because the secondnearest neighbor Fe atoms in orthorhombic SFO are octahedrally coordinated and have 5 electrons
(each one occupying one of the electron states in Figure 4.5(b)) it is not energetically favorable
for the electrons formed by oxygen vacancy formation to be transferred to the second nearest
neighbors. Further, since localization of these extra electrons on the nearest SP-Fe would
significantly increase the energy of the system, the FeO4 SP coordination undergoes a structural
transformation to tetrahedral (Td) coordination (with a comparatively lower energy state for 5
unpaired d-orbital electrons, as shown in Figure 4.5(c)), in the process causing the orthorhombic
to brownmillerite phase transition.

59

Figure 4.4. Oxygen vacancy polaron shape in (a) tetragonal, (b) orthorhombic, (c) brownmillerite
SrFeO3-δ. Fe-O polyhedra colored according to the charge on Fe.27 Here, blue atoms are Sr,
brown are Fe and red are O.

Figure 4.5. d-orbital splitting and energy distribution for different Fe-O polyhedra was derived
from crystal field theory. (a) Square planar Fe-O coordination, (b) octahedral FeO coordination
(c) tetrahedral Fe-O coordination.
For dilute oxygen vacancy formation (δ < 0.1), the vacancy formation energies for cubic
SFO, rhombohedral LSF55, and orthorhombic LFO are 0.4, 1.6, and 3.5 eV respectively. This ~3
eV variation shows that even in the dilute range the La/Sr ratio dramatically impacts the oxygen
vacancy concentration. Furthermore, it is significant to note that maximum oxygen vacancy site
fraction that can be achieved in SFO, LSF, and LFO within SOFC operating conditions is ~ 10%
that occurs at 600oC in orthorhombic SFO.160
60

4.2.3. Rhombohedral and Cubic La0.5Sr0.5FeO3-δ
In both stoichiometric rhombohedral and stoichiometric cubic LSF55, the computationally
predicted charge state on Fe is 3.5+. This indicates that, on average, formally there are four and a
half electrons in the Fe d-orbitals. One could imagine, as a combination of LFO and SFO, the
excess electrons due to oxygen vacancy formation could, therefore, be distributed to both the first
and second nearest neighboring Fe atoms and this is indeed the case, as shown in Figures 4.2(f)
and 4.2(g).
In cubic LSF55, after oxygen vacancy formation all the Oh-Fe connected to the SP-Fe via
two 180˚ Fe-O-Fe bonds bend to 160o. The charge on first and second nearest neighboring Fe
atoms were 3.35+, suggesting that the electrons left behind by neutral oxygen vacancy formation
were distributed to the first (two SP-Fe) and second (eight Oh-Fe) nearest neighboring Fe atoms
to the oxygen vacancy. Thus, as shown in Figures 4.2(c) and 4.2(g), in cubic LSF55 the oxygen
vacancy polaron exhibit a “pancake” shaped polaron with an average diameter and height of ~8.3
and ~4.0 Å, respectively.
In rhombohedral LSF55, the polaron size is smaller than that in cubic LSF55, due to larger
Fe-O polyhedra tilt angle in the rhombohedral structure. The oxygen vacancy formation
calculations showed that oxygen vacancy formation in the Sr-O3 plane was energetically favorable,
compared to oxygen vacancy formation in the La-O3 plane, by ~0.1 eV. After oxygen vacancy
formation, the two SP Fe charge changed to 3.6+ from 3.5+. In contrast, the four Oh-Fe connected
with the tilted Fe-O-Fe bond shown in Figure 4.2(c) gained electrons and became 3.3+ from 3.5+.
This difference is related to the Fe-O-Fe bond angle. Note there are four Oh-Fe connecting the
four O atoms at the square bottom of the SP-Fe polyhedra, via a Fe-O-Fe bridge. In rhombohedral
LSF55, two of the Fe-O-Fe bond angles are ~170o (compared to 180˚ in the cubic structure) and

61

the other two Fe-O-Fe bond angles are ~153o (i.e they are more bent). Since a more bent bond
angle results in more Jahn-Teller distortions, the split eg energy level in the more tilted Oh will be
at a lower energy state and the excess electrons generated by oxygen vacancy formation will
preferably transfer to the more titled Oh-Fe atoms, as shown in Figures 4.2(c) and 4.2(g). In
rhombohedral LSF55, these excess electrons produced a quarter of a “pancake” like polaron shape
(Figure 4.2(f)). The polaron size was approximately one-fourth of that in cubic LSF55, with an
average polaron radius and height of 4.3 and 4.0 Å, respectively.
As shown in Figure 4.2, the size of the polaron is the smallest in LFO and it increases with
increasing Sr doping LSF. The small and more localized polaron in LFO can screen the charge of
the oxygen vacancy more effectively, thus the oxygen vacancies will remain non-interacting and
will not impact each other’s polaronic strain field until they are very close. However, the larger
and more delocalized polarons can not screen the vacancy charge effectively, and thus the
vacancies may start to interact, even at a large distance. Therefore, in the next section, the oxygen
vacancy formation energy as a function of vacancy concentration was studied and correlated with
the polaron size.

62

5.

Oxygen Vacancy Formation sites and Formation Energy for Interacting and Noninteracting Vacancies at 0 K

In most of the LSF phases, multiple oxygen vacancy formation sites are available, observed from
the symmetry. For vacancy formation calculation it is important to find all the possible vacancy
formation sites and after that energetically favorable vacancy formation site which will give
oxygen vacancy formation energy for that structure.

5.1. Oxygen Vacancy Formation Sites in SrFeO3-δ Phases
The dashed red circles in Figure 1.1 identify the preferred oxygen vacancy generation sites, as
determined by comparing the Table 5.1 oxygen vacancy formation energies for the various oxygen
sites within each crystal structure. Figure 1.1(a) shows that in a perfect cubic SrFeO3, all of the
oxygen atoms are situated at the corners of symmetrical Fe-O octahedra, 1.94 Å from the central
Fe ion, and shared between two neighboring octahedra (here the designation Oh is used to denote
an octahedron, while the oxygen site at the corner of two neighboring octahedra is denoted as OhOh). The high order of geometric symmetry (point group: Oh) and the single Fe4+ oxidation state
ensures that only one type of oxygen site (denoted O1) is present in cubic SrFeO3. Also, according
to MĂśssbauer spectroscopic study, the charge on Fe in cubic SrFeO3 does not show a charge
disproportionation reaction, like what happens in CaFeO3 (where Fe4+  Fe3+ + Fe5+).161 When
an oxygen atom leaves the Oh-Oh site in cubic SrFeO3, an oxygen vacancy is formed and the two
neighboring Oh-Oh polyhedra collapse into two square pyramidal (SP) polyhedra. Curiously, for
reasons described in Section 3.3.4, the two iron atoms bordering this oxygen vacancy (the ones
now in SP coordination) remain as Fe4+.
In contrast, Figure 1.1(b) shows that in perfect tetragonal SrFeO2.875 both octahedral (Oh)
and square pyramidal (SP) polyhedra are present. Charge interpretation based on the magnetic
63

moment of iron indicates that the SP coordinated iron is Fe4+ while the Oh coordinated iron is
either Fe3.6+ (denoted as Oh1) or Fe3.9+ (denoted as Oh2). The Oh1 octahedra are deformed and
tilted (producing a Fe-O-Fe bond angle of 168.67o between two Oh1 octahedra) while there is no
tilt between two Oh2 octahedra (producing a Fe-O-Fe bond angle of 180˚, just as in the cubic
SrFeO3 structure). As shown in Table 5.1, the most favorable oxygen vacancy formation site in
tetragonal SrFeO2.875 (denoted as the O6 or Oh2-Oh2 site) is at the corner shared by two Oh2
octahedra.
Table 5.1. The vacancy formation energies at different oxygen site calculated for all four
crystalline structures. Each site is noted by its label in Figure 1.1, its position in Fe-O
polyhedra, and the charge of the Fe ions in the Fe-O polyhedra. The bold values indicate
the most favorable oxygen vacancy site.
Compound

Label of O in
Figure 1.1

Position in
Fe-O Polyhedra

The charge state
of Fe

ΔEvac (eV)

SrFeO3
(cubic)

O1

Oh-Oh

+4 --- +4

0.71

O2
O3
O4
O5
O6
O7
O8
O9
O10
O11
O12

Oh1-Oh2
Oh1-SP
Oh1-Oh1
SP-SP
Oh2-Oh2
Oh-SP
Oh-Oh
SP-SP
Oh-Oh
Oh-Td
Td-Td

+3.6 --- +3.9
+3.6 --- +4
+3.6 --- +3.6
+4 --- +4
+3.9 --- +3.9
+3.3 --- +3.9
+3.3 --- +3.3
+3.9 --- +3.9
+3 --- +3
+3 --- +3.2
+3.2 --- +3.2

1.06
1.05
1.14
1.15
0.92
2.07
2.12
1.25
3.44
3.28
3.00

SrFeO2.875
(tetragonal)

SrFeO2.75
(orthorhombic)
SrFeO2.5
(brownmillerite)

Figure 1.1(c) shows that in the perfect orthorhombic SrFeO2.75 structure only a single type of Oh
polyhedron (all containing a Fe3.3+ ion) and a single type of SP polyhedron (all containing a Fe3.9+
ion) are present. These results clearly show that this structure is not made up of iron atoms that are
distinctly 3+ or 4+ in charge, as has been postulated using the concept of formal charge,32 but
instead made up of iron atoms with fractional oxidation states interpreted from the calculated
64

magnetic moment. This results from internal redox reactions with the neighboring oxygen atoms
that result in informal charges on the oxygen atoms.162 Comparing Figure 1.1(b) and 1.1(c), the
orthorhombic strontium ferrite can be considered as tetragonal strontium ferrite which has lost one
oxygen from the Oh2-Oh2 site and where the two neighboring Oh2-Oh2 polyhedra have collapsed
into two SP polyhedra. The two iron atoms in these SP polyhedra remain ~4+, while the iron atoms
in the neighboring Oh sites reduce their charge to 3.3+ with the transformation from tetragonal to
orthorhombic SrFeO2.75. In the orthorhombic structure, all of the Oh polyhedra are deformed and
tilted (producing a Fe-O-Fe bond angle of 165.22o between two Oh polyhedra). As shown in Table
5.1, the most favorable oxygen vacancy formation site in orthorhombic SrFeO2.75 is at the shared
corner of two SP polyhedra.
Figure 1.1(d) shows that in the brownmillerite SrFeO2.5 structure both Oh and tetrahedral
(Td) Fe-O coordination is present. Although some literature sources have claimed that all Fe ions
in SrFeO2.5 carry a formal charge state of Fe3+ in SrFeO2.5,32 the analysis here indicates that Td Fe
B

and a charge of 3.2+, while the Oh Fe has a magnetic moment

+. This agrees with experimentally measured MĂśssbauer data
indicating that the magnetic moment of Oh Fe is slightly higher than that on Td Fe.148 As shown
in Table 5.1, the most favorable oxygen vacancy formation site in brownmillerite SrFeO2.5 is at
the shared corner of two Td polyhedra (at O12).
The charge disproportionation of Fe in oxygen deficient SrFeO3-δ phases (Figure 1.1) is
accompanied by gradual loss of occupied states near the Fermi level (Figure 3.2), causing a
reduction in the electronic conductivity. This phenomenon has been reported by Hirai et al. during
charge disproportion induced phase transformation in tetragonal SrFeO3-δ at 70 K.163 As
demonstrated in Table 5.1, in each of the strontium ferrite structures modeled here, oxygen atoms
65

bonded to the highest charged Fe atoms in polyhedra are always the most favorable locations for
oxygen vacancy formation.
5.2. Oxygen Vacancy Formation Energies in SrFeO3-δ Phases

Figure 5.1. The vacancy formation energy per oxygen vacancy, calculated from DFT for four
different SrFeO3-δ phases.
𝑓

Figure 5.1 shows the oxygen vacancy formation energy (∆𝐸𝑣𝑎𝑐 ) calculated using Equation 2.1 (on
a per vacancy basis) for cubic, tetragonal, orthorhombic, and brownmillerite strontium ferrite as a
function of oxygen nonstoichiometry (plotted in terms of the absolute nonthese calculations, only the site with the lowest oxygen vacancy formation energy in each phase
identified in Table 5.1 was utilized. It is important to note that the trends observed here should
also hold at high temperature since Equation 2.4 indicates that oxygen’s temperature and pressure
effects will shift the formation energy in each phase by the same amount. Figure 5.1 reveals that
the 0 K oxygen vacancy formation energy increases with each phase transition. This is because,
as discussed in Section 4.1, oxygen vacancy formation is always the easiest next to the iron with
66

higher oxidation state, which decreases in quantity with each new structure. It is important to note
that the vacancy formation energy in cubic SrFeO3-δ is a very low 0.4 eV, which explains why the
material has such a high oxygen vacancy concentration at room temperature in air and why it has
been difficult to produce stoichiometric cubic SrFeO3.32,44 Figure 5.1 also shows that the oxygen
vacancy formation energy in cubic SrFeO3-

<0.05;

increases to ~0.5 eV for 0.05< δ <0.1; and quickly increases with concentration when δ >~0.1.
The increase in oxygen vacancy formation energy with increasing oxygen non-stoichiometry is
due to oxygen vacancy interactions. A similar trend due to oxygen vacancy interactions was also
reported by Kuklja et al.29 in their computational work on cubic and hexagonal
BaxSr1−xCoyFe1−yO3−δ. None of the other SrFeO3-

phases (tetragonal, orthorhombic and

brownmillerite SrFeO3energy on the oxygen vacancy concentration. The concentration dependent oxygen vacancy
formation energy in cubic SrFeO3-

,

when the average distance between oxygen vacancies are ~11 Å, and the interaction becomes
stronger when δ > 0.1 as the average oxygen vacancy distance decreases to ~8 Å. This average
distance between neighboring oxygen vacancies is more than 4 times larger than the Fe-O bond
distance (1.94 Å).

67

5.3. Computationally Predicted Vacancy Ordered Phase Transition in SrFeO3-δ Phases

𝛿

Figure 5.2. Energy required (∆𝐸 𝛿 ) to form the non-stoichiometric structures (𝑆𝑟𝐹𝑒𝑂3−𝛿 + 2 𝑂2 )
from cubic SrFeO3.
Experimentally, with increasing T and the increasing δ accompanying it, strontium ferrite displays
several phase transitions as it progresses from cubic SrFeO3 to tetragonal SrFeO2.875 (with
𝛿 0 =0.125), to orthorhombic SrFeO2.75 (with 𝛿 0 =0.25), and finally to brownmillerite SrFeO2.5 (with
𝛿 0 =0.5). However, it has not been clear if these phase changes are induced by T or by δ. Thankfully,
DFT calculations such as those in Figure 5.2 where the X for any structure can be arbitrarily
increased as long as δ is > 𝛿 𝑜 for the structure of interest (including those not observed in nature
because of a phase transformation to a lower energy crystal structure) allow the effects of oxygen
non-stoichiometry alone to be examined. Specifically, Figure 5.2 shows that at 0 K cubic strontium
ferrite has the lowest energy from 0< δ <0.125, tetragonal strontium ferrite has the lowest energy
from 0.125<𝛿<0.25, tetragonal and orthorhombic strontium ferrite have similar energies from
0.25<𝛿<0.50, and brownmillerite has the lowest energy for 𝛿>0.50. This suggests that oxygen
68

non-stoichiometry instead of temperature is the dominant force driving the strontium ferrite phase
transformations observed from 0-1400 oC in air. Further, these DFT calculations indicate that the
enthalpy gain in the cubic to tetragonal or orthorhombic to brownmillerite phase transitions are
larger than the configurational entropy contributions.

Figure 5.3. GGA+U calculated oxygen vacancy formation energies for different LaxSr1–xFeO3
compositions. The Fe-O polyhedra represent the charge on Fe in a perfect lattice. The lines are
guides to the eye.
In Figure 5.3, the DFT+U calculated oxygen vacancy formation energies for various LSF
phases are plotted versus increasing oxygen non-stoichiometry per formula unit (i.e δ). Vacancy
formation energies of cubic SFO replotted from Figure 5.1 for comparison purpose. The vacancy
formation energies in LFO remained constant at 3.57 Âą 0.04 eV and did not show any significant
interaction even at high oxygen non-stoichiometry (i.e. at δ = 0.25). This is consistent with the
small and localized vacancy LFO polaron that allows the oxygen vacancies to be very close
without interacting (and hence treated as dilute).

69

Figure 5.3 also shows that rhombohedral LSF55, which as shown in Figure 4.2 has tilted
Fe-O polyhedra and a small polaron size compared to the untilted Fe-O polyhedra and large
polaron sizes found in cubic LSF55, exhibits a lower oxygen vacancy formation energy and less
oxygen vacancy interactions than cubic LSF55. The vacancy formation energy in rhombohedral
LSF55 remained constant at ~1.7 eV for δ ≤ 0.25 until the polarons began to overlap at an average
distance between oxygen vacancies of ~7 Å, larger than the polaron size due to its irregular shape.
In contrast, cubic LSF55 showed an increasing trend even for δ <0.1, due to its much larger polaron
size. Note, in cubic LSF55 it was not computationally possible to find a dilute limit for oxygen
vacancy formation. This is because the atomic position relaxation of the cubic LSF55 supercell
with very low oxygen non-stoichiometry (δ = 0.03) resulted in increased octahedral tilt around
each oxygen vacancy site producing a locally rhombohedral LSF55 structure. This resulted in the
calculated cubic LSF55 oxygen vacancy formation energy that was less than that found in
rhombohedral LSF55 (1.4 eV vs.1.7 eV) at very low oxygen non-stoichiometry. This suggests that
cubic LSF55 may only exist at higher δ and/or temperature values far from the rhombohedral to
cubic LSF phase boundary.
Overall, the oxygen vacancy formation energies in rhombohedral and cubic LSF55 suggest
that smaller polaron size leads to weaker interaction among oxygen vacancies. Further, the
observed trends with La:Sr ratio indicate that progressively larger polaron sizes result in increasing
oxygen vacancy formation energies with increasing Sr or δ. This behavior is consistent with
experimental results showing that Sr-rich LSF and highly oxygen deficient compositions are more
prone to oxygen-vacancy-ordering-induced phase transformations, while Sr-poor and oxygen rich
LSF compositions are not.63

70

In addition, the crystal-field-theory based polaron size estimation technique presented here
can be applied to other perovskite structures. For example, LaMnO3 has a structure similar to
LaFeO3, which shows almost no oxygen vacancy interactions. However, Mn3+ has the same
number of d-electrons as Fe4+ (i.e. both have an electron configuration of [Ar] 3d4). Therefore, the
excess electrons generated by oxygen vacancy formation should be distributed to the Ohcoordinated Fe, and thus the long-range charge transfer observed in SrFeO3 should also occur in
LaMnO3 and result in a large oxygen vacancy polaron size that leads to strong vacancy interactions
et al.28

71

6.

Oxygen Vacancy Site Fraction as a Function of Temperature and Phase Change (in
air)

Figure 6.1. Dotted line represent δ increasing with T at pO2 = 0.21 atm according to Equation
𝑓

1.6 (dilute vacancy assumption) for cubic SrFeO3-δ (∆𝐸𝑣𝑎𝑐 from Figure 5.1 for different δ). Solid
line represent δ increasing with T at atmospheric condition according to Equation 2.17
𝑓

(interacting vacancy assumption) for cubic SrFeO3-δ (where ∆𝐸𝑣𝑎𝑐 is a function of δ).
It is important to correctly account for the concertation dependence of the oxygen vacancy
formation energy when calculating the equilibrium oxygen non-stoichiometry. Specifically, Figure
6.1 shows that using the dilute model and assuming concentration independent oxygen vacancy
formation energy values for cubic SrFeO3-δ (the dashed lines) to predict oxygen non-stoichiometry
values (i.e. following Equation 1.6 or 2.16) leads to results that are wildly different from
experimentally observed values and sometimes unphysical. For instance, assuming ΔEvac = 0.4 eV
for cubic SrFeO3- (the smallest value calculated in Figure 5.1) results in a δ ≤ 0.02 at 0 oC that
quickly rises to ~3 at 300oC after applying the thermodynamic connection energy of T at 𝑝𝑂2 =
0.21 atm. This result is unphysical because a δ=3 would indicate that cubic strontium ferrite would
72

be completely reduced to its constituent metals at 300oC, despite the fact that strontium ferrite is
stable to well above 800oC.44,158 Similarly, assuming a concentration independent cubic strontium
ferrite oxygen vacancy formation energy of ΔEvac = 1.45 eV (the largest value calculated for cubic
SrFeO3-

in Figure 5.1) and the necessary thermodynamic connection energy predict that cubic

strontium ferrite is stable at high temperature, and that the predicted oxygen non-stoichiometry at
25oC in air is many orders of magnitude below experimental value.44 It is important to note that
when oxygen vacancies interact, their formation energy increases with the site fraction X
𝑓

(especially in cubic phase). Therefore ∆𝐺𝑣𝑎𝑐 (𝑇, 𝑝, 𝑋) depends on temperature, oxygen partial
pressure, and the vacancy site fraction, as demonstrated in Figure 6.2 and in previous studies.24,164
6.1. Oxygen Vacancy Site Fractions in Different SrFeO3-δ Phases Predicted with
Temperature in Air

Figure 6.2. (a) Vacancy formation free energy with nonstoichiometry for all the phases (b)
Nonstoichiometry is plotted as a function of temperature. (c) Vacancy formation site fraction is
plotted as a function of temperature. (the data point corresponding to 3Xcubic = 0.5 is neglected
from fitting as at this vacancy site fraction cubic SrFeO3-δ is easily transformed to brownmillerite
and hence loses its available vacancy sites).

73

Figure 6.2 (cont’d)

It is therefore important to apply the non-dilute model and solve Equation 2.17 to calculate the
oxygen vacancy site fraction X. The solid line in Figure 6.1 shows that much more realistic oxygen

74

𝑓

non-stoichiometry values can be obtained using the full range of concentration-dependent ∆𝐸𝑣𝑎𝑐
values obtained in Figure 5.1 using Equation 2.17.
Figure 6.2(a) shows the full range of temperature and oxygen partial pressure corrected
𝑓

concentration-dependent free energy of oxygen vacancy formation, ∆𝐺𝑣𝑎𝑐 (𝑇, 𝑝, 𝑋), calculated via
Equation 2.4 for the various strontium ferrite phases as a function of oxygen non-stoichiometry.
𝑓

𝑓

Although only modest increases in ∆𝐺𝑣𝑎𝑐 occur with  for  < 0.5, the ∆𝐺𝑣𝑎𝑐 jumps dramatically
once the brownmillerite phase forms.
Based on the solution of Equation 2.17, a unique relationship between X and T can be

Figure 6.2(b) compares the experimentally measured and our predicted oxygen non-stoichiometry
of strontium ferrites (a trend line is added as the moving average) from 300-1300oC in the air. The
computationally predicted oxygen non-stoichiometry vs. temperature is similar to the
experimentally measured values of both Takeda et al.44 and Vashuk et al.57 Such agreement
demonstrates that the DFT+U predicted oxygen vacancy formation sites and activation energies
reported here are correct and it is important to account for the concentration dependent vacancy
formation energy in predicting vacancy concentrations. That said, some deviation between the
computationally-predicted and experimentally-measured oxygen non-stoichiometry exists. This
deviation may result from multiple phases can coexist within several oxygen non-stoichiometry
windows.44
Lastly, Figure 6.2(c) shows the strontium ferrite oxygen vacancy site fraction as a function
of temperature. Interestingly, the oxygen vacancy site fraction (which together with the charge and
oxygen vacancy diffusivity determines the oxygen vacancy concentration important to so many

75

MIEC applications) peaks around 600oC then declines at higher temperatures due to the much
larger oxygen vacancy formation energy in the brownmillerite structure compared to the cubic,
tetragonal, and orthorhombic structures (Figure 5.1). As will be explored in a future publication,165
since the cubic structure shows the lowest oxygen vacancy formation energy but suffers from
strong vacancy interactions caused by the long-range electron transfer, it is reasonable to suggest
that dopants (such as La on the Sr site as occurs in the commonly used SOFC material
La0.6Sr0.4FeO3-δ)6 that can localize electrons around an oxygen vacancy may reduce the interaction
between oxygen vacancies and eliminate vacancy ordering, allowing for cubic SrFeO3- with high
ionic conductivity. This possibility will be explored in a future paper.
6.2. Oxygen Vacancy Site Fractions in Different LSF Compositions Predicted with
Temperature in Air

Figure 6.3 shows the oxygen vacancy site fraction as a function of temperature and partial pressure
of oxygen (pO2 = 0.21 atm) for three different LSF compositions calculated using the non-dilute
oxygen vacancy site fraction calculation method of Equation 6. This method was developed and
X as a function of temperature in SFO in the air, and (considering
phase changes in the material) the predicted trends compared well with the experimental TGA
data.
Figure 6.3(a) compares the calculated oxygen vacancy non-stoichiometry, δ, in
rhombohedral and cubic LSF55 structures with experimentally measured data11 and shows
reasonably good agreement at SOFC operating temperatures (600 ˚C ~1200 ˚C ) The experimental
data in Figure 6.3(a)
decreases with increasing temperature. This is caused by rhombohedral to cubic phase

76

transformation at high temperature65 and the stronger oxygen vacancy interactions in cubic LSF55
than in rhombohedral LSF55 at low temperature.

Figure 6.3. (a) The oxygen non-stoichiometry versus temperature in air,(b) the oxygen vacancy
site fraction as a function of temperature in air, and (c) the oxygen vacancy site fraction as a
function of temperature and La/Sr ratio for various LaxSr1–xFeO3 phases. Note, the lines in (a-b)
and are guides to the eye and the lines in (c) are the experimentally determined phase
boundaries.44,65 The oxygen vacancy site fraction in (c) is denoted by the color overlay and the
right-hand color scale. The X values in b) and c) are never zero. They are just smaller than can
be resolved with the linear X scale chosen to highlight the X trends discussed here.

77

Figure 6.3 (cont’d)

As discussed previously, rhombohedral LSF55 generally has a lower oxygen vacancy
formation energy, and therefore Figure 6.3(a)

rhombohedral LSF than for

cubic LSF at similar temperatures. The rhombohedral to cubic LSF phase transformation
temperature increases with higher La:Sr ratio, as shown in the experimentally overlaid phase
diagrams of Figure 6.3(c). For example, the phase transition temperature occurs at ~450oC for
LSF55 but is increased to 800˚C for La0.6Sr0.4FeO3-δ (LSF64). Hence, it is beneficial to adjust the
La:Sr ratio so that a highly oxygen deficient rhombohedral LSF, instead of cubic LSF, is formed.
This likely contributes to the high performance of La0.6Sr0.4FeO3-δ (LSF64) at traditional 800oC
SOFC operation, as LSF64 retains the rhombohedral phase at 800oC according to the experimental
phase diagram,
Considering the various phase transitions occurring in LSF, the oxygen vacancy site fraction, X,
was computed for the specific phase occurring at that temperature according to the experimental
phase diagrams from Fossdal et al.,65 and Takeda et al.44 The computed X values for LFO, SFO,
78

and LSF55 are plotted in Figure 6.3(b). Figure 6.3(b) shows that due to the oxygen vacancy
formation energies are shown in Figure 5.3 and those calculated for other phases, increasing La:Sr
ratio tends to result in a lower X at a given temperature. In LFO, there were almost no oxygen
vacancy sites till 1000oC. In contrast, SFO is highly oxygen deficient SFO but experiences a
reduction in X with increasing temperature due to oxygen vacancy ordering once brownmillerite
is formed. Figure 6.3(b) also shows that the X in orthorhombic SrFeO3- , is larger than in LSF55
at ~600oC. Additional studies on the oxygen ion conductivity are performed to determine if this
results in a higher (~600oC) SFO ionic conductivity than LSF55 and reported in Chapter 7.31,160
6.3.

Effect of the Variation of U-Parameter on Oxygen Vacancy Formation

The qualitative trends in oxygen vacancy formation energy with vacancy concentration do not
depend on the selection of U parameter, and that Ueff = 3 gives the best match to the measured
oxygen nonstoichiometry. This is shown in Figure 6.4(b) where the oxygen vacancy formation
energies for four different phases at different oxygen nonstoichiometries using Ueff = 2, 3, and 4
are shown. The solid black line shows the trend of oxygen vacancy formation energy with
increasing oxygen nonstoichiometry at Ueff = 3 as presented in Figure 6.2(b). The dotted and
dashed line showed results corresponding to Ueff = 2 and 4 respectively. It is important to note that
oxygen vacancy formation energies calculated with Ueff = 2 overestimate the results for cubic and
tetragonal phase by ~ 0.3 eV, whereas, calculations made with Ueff = 4 underestimates the results
by ~ 0.3 eV. Further, the effect of a change in U parameter (from 2 to 3, or 3 to 4) on the calculated
oxygen vacancy formation energies for orthorhombic and brownmillerite phases were insignificant
(< 0.1 eV).

79

Figure 6.4. (a) Compare the impact of U parameters on oxygen vacancy formation energy as a
function of oxygen nonstoichiometry; (b) compare the calculated oxygen nonstoichiometry as a
function of temperature with different U parameters with experimental values.
The oxygen nonstoichiometry and temperature relation, δ(T) were also calculated with Ueff = 2 and
Ueff = 4 and plotted in Figure 6.4(b) and compared with the results of Ueff = 3. At lower oxygen
nonstoichiometries (δ = ~0 – 0.3) they respectively overestimated and underestimated the
80

temperature, required to generate the experimentally reported δ, by ~200oC. Overall, Ueff = 3 led
to predicted δ(T) relationship more consistent with experimental data for SrFeO3-δ phases, mainly
due to it is suited for both Fe+3 and Fe+4.

81

7.

Oxygen Ionic Conductivity in La1-xSrxFeO3-δ – Considering Oxygen Vacancy
Interaction and Phase Change

Oxygen vacancy concentration and diffusion in any perovskite is central to its ion conduction
capability. The product of oxygen vacancy concentration, diffusivity, and charge on the oxygen
ion gives the total ionic conductivity (Equation 1.1). Chapter 6 showed the oxygen vacancy
concentration as a function of La/Sr ratio and temperature. Present chapter will focus on how
diffusivity and ionic conductivity vary with La/Sr ratio and try to locate the composition with the
maximum ionic conductivity at SOFC operating conditions. Patrakeev et al.62 have shown LSF55
possesses the highest oxygen ionic conductivity within LSF family at 950oC in the air in their
experimental study. Whereas, most of the previous experimental studies to improve ionic
conductivity or oxygen surface exchange in LSF, were performed with a different composition,
La0.6Sr0.4FeO3-δ (LSF64).6,47,166–169 With this computational study on ionic conductivity we have
explained the reason for the selection of LSF64 over LSF55.

Oxygen vacancy formation energy and oxygen migration barrier are the major controlling
factors for oxygen vacancy concentration and diffusion respectively, while temperature and
pressure kept constant. Computational study on the oxygen vacancy formation energy shows, it
varies by ~3 eV with change in La/Sr ratio at the A site of ABO3-perovskite, LSF (Figure 5.3).
Also, significant variation (~2.5 eV) was observed in oxygen vacancy formation energy due to
phase change in LSF end member SFO (Figure 5.1). Whereas, previous experimental and
computational studies on oxygen vacancy migration in LSF predicted that A-site doping has
negligible effect on oxygen migration barrier in LSF.18,20,43,55 This was confirmed from the
computational work of Mastrikov et al.55 They mentioned in their paper, “In (La,Sr)(Mn,Fe,Co)O3-

82

δ

perovskites, the oxygen vacancy diffusion coefficients were experimentally found to be almost

independent of the cation composition, with a typical migration barrier of ≈ 0.8 eV”. However,
these earlier studies on oxygen migration barrier in LSF are not complete for a broader range of
La/Sr ratio, phase change and/or varying non-stoichiometry in LSF.

In this Chapter we have computationally evaluated oxygen vacancy migration barrier in
different LSF compositions and phases, then predicted oxygen ion diffusivity, mobility and
ultimately oxygen ionic conductivity to understand the effect of oxygen vacancy concentration
and diffusion barrier separately, on the oxygen ion conduction. Previously calculated oxygen
vacancy site fraction data is used to calculate oxygen vacancy concentration as a function of
temperature, in the air.

7.1

Oxygen Vacancy Migration in La1-xSrxFeO3-δ
It was shown by Chroneos et al.170 that oxygen ion/vacancy moves via a vacancy hopping

mechanism in a perovskite, like LSF. Using the CINEB calculation method described in Section
2.3.3 oxygen ion migration barriers were calculated in different LSF phases. From the oxygen
vacancy site fraction calculations shown in Figure 6.2 and 6.3 within SOFC operating temperature
(in air), it can be noticed the oxygen vacancy concentration remains close to dilute in most phases,
except high temperature cubic SFO and cubic LSF55.
Therefore, oxygen vacancy migration barrier was calculated at different representative
concentrations in all phases including: orthorhombic LFO, rhombohedral LSF55, cubic LSF55,
cubic SFO, tetragonal SFO, and orthorhombic SFO. Due to the extremely high oxygen vacancy

83

formation energy of ~3 eV in brownmillerite SFO, this material is of less interest for MIEC
applications, and therefore its migration barrier calculations were not performed.

7.1.1

Oxygen Vacancy Migration in LaFeO3-δ
Figure 7.1(a) shows the oxygen migration pathway in LFO, highlighted with pink atoms.

Lattice distortion shown in the figure corresponds to the transition state lattice distortion.
Calculated migration barrier at two different δ (= 0.06 and 0.25), shown in Figure 7.1(b), possess
the same migration energy barrier of 0.89 eV. This confirms that even at δ=0.25, the vacancies do
not interact in LFO, which agrees with our dilute vacancy concept in LFO, described in Section
4.2.1. This calculated oxygen vacancy migration barrier in LFO is comparable to experimentally
observed oxygen migration barrier of 0.77 eV obtained via oxygen tracer diffusion43 and 1.1 eV
obtained in conductivity relaxation experiment.171 Figure 7.1(c) shows how the magnetic moment
changes on the oxygen-vacancy-neighboring Fe atoms during oxygen vacancy migration. This
calculation shows that magnetic moment interpreted oxidation state of the second nearest
neighboring Fe remains constant at 3+ during oxygen vacancy migration.27 This finding is
consistent with the small polaron volume calculated earlier for this material. Further, it can be
noticed that when an oxygen atom moved from Îś = 0 to Îś = 1, the Fe2 atom reduced its charge
from 3+ to 2+, as it changes from the 2NN to the 1NN of the oxygen vacancy. The magnetic
moment on Fe1, oxygen vacancy adjacent Fe, maintained same magnetic moment (corresponding
to Fe2+) before and after migration, while at the transition state it was reduced due to gain in
electron but the oxidation state was beyond our scale of measurement, defined in Section 2.1.2.

84

Figure 7.1. Oxygen migration in LaFeO3-δ. The pink atoms in (a) highlight the oxygen migration
path, (b) shows the oxygen migration barrier along this path for two different δ, and (c) shows the
charge change on the neighboring Fe atoms as a function of oxygen atom coordinate on its path
(at δ =0.06).

7.1.2

Oxygen Vacancy Migration in La0.5Sr0.5FeO3-δ

Figure 7.2(a) and (b) show the oxygen migration path in rhombohedral and cubic LSF55
respectively. The calculated oxygen migration barriers are shown in Figure 7.2(c), for
rhombohedral and cubic LSF55. The calculated migration barrier in rhombohedral LSF55 was
0.45 eV for both dilute and interacting vacancies (at δ = 0.04 and 0.17). The calculated migration
barrier in cubic LSF55 was 0.7 eV at δ = 0.12. Due to the lattice distortion in cubic LSF55, lattice
locally transformed towards “rhombohedral LSF55” distortion with a significant reduction in
structural energy at ζ = 0.25/0.75 which can be called as “trapped state” in vacancy migration path.
85

These results show, the effect of phase change on the oxygen vacancy migration barrier is not as
substantial (differ by ~0.25eV), as it was on oxygen vacancy formation energy (differ by ~0.5eV
Figure 5.3).

In Figure 7.2(d) and (e), change of the magnetic moment on oxygen vacancy neighboring
Fe atoms are shown in rhombohedral and cubic LSF55, respectively. Oxygen vacancy adjacent
Fe1 maintained the same oxidation state before and after oxygen vacancy migration in both,
rhombohedral and cubic LSF55. Whereas, the oxidation state on Fe4 (second nearest neighbor to
oxygen vacancy) was reduced due to partial charge transfer after oxygen vacancy migration in
rhombohedral LSF55. After the migration of the oxygen vacancy from Fe3 adjacent site to Fe2
adjacent site, the oxidation state on Fe2 and Fe3 interchanged in cubic LSF55. However, the
magnetic moments Fe2 and Fe3 did not simply interchange in rhombohedral LSF55, due to the
anisotropic shape of polaron (Figure 4.2(b) and (f)). Due to symmetric polaron shape in cubic
LSF55, the magnetic moment on Fe4 did not change, as it remains as the second nearest neighbor
to the oxygen vacancy. A sudden dip in calculated magnetic moment or increase in Fe oxidation
state can be observed at ζ = 0.25/0.75, due to the “trapped state” in the migration path. The larger
change of the magnetic moment on Fe is in coordination with the higher migration energy barrier
in cubic LSF55.

86

Figure 7.2. Oxygen migration in LSF55. The pink atoms in (a) highlight the oxygen migration
path in Rhombohedral LSF55 with δ =0.04, and (b) highlight the oxygen migration path in cubic
LSF55 with δ =0.13. Figure (c) shows the oxygen migration barrier in rhombohedral and cubic
LSF55. (d) shows the charge change on the oxygen vacancy neighboring Fe atoms in
rhombohedral LSF55 as a function of oxygen atom coordinate (corresponds to Figure (a)). Figure
(e) shows the charge change on the oxygen vacancy neighboring Fe atoms in cubic LSF55 as a
function of oxygen atom coordinate (corresponds to Figure (b)).

87

7.1.3

Oxygen Vacancy Migration in SrFeO3-δ

Figure 7.3(a) shows the oxygen migration path (pink atoms) in cubic SFO. The lattice
distortion shown in the figure corresponds to the transition state lattice distortion. As shown in
Figure 7.3(b), the calculated oxygen migration barriers for cubic SFO were 0.58, 0.62 and 1.07
eV at δ = 0.04, 0.12, and 0.5 respectively. These results show that the effect of oxygen
nonstoichiometry on the oxygen vacancy migration barrier is insignificant (<0.05 eV) below a
critical nonstoichiometry, above which oxygen vacancy strongly interacts (at δ = 0.5). From the
oxygen vacancy site fraction calculations shown in Figure 6.2 and 6.3 within SOFC operating
temperature (in air), it can be noticed the oxygen vacancy concentration is less than ~0.125 in
cubic SFO until it transforms to tetragonal phase. Below this concentration, the migrantion energy
barrier is ~0.6eV and does not change with δ.

Figure 7.3(c) shows the magnetic moment change on neighboring Fe atoms during the
migration. The magnetic moment on vacancy adjacent Fe1 goes through a minimum during
oxygen vacancy migration, i.e., Fe oxidation state was increased at transition state during oxygen
vacancy migration, because excess electrons were pushed to the second nearest neighboring Fe4.
Fe1 maintained same oxidation state, ~4.1+ before and after oxygen vacancy migration due to its
square pyramidal coordination with neighboring oxygen atoms. Initially (Îś = 0), the oxygen
vacancy was between Fe1 and Fe3. Therefore, the oxidation state on square pyramidal-Fe1 and
Fe3 were same, (~4.1+) corresponding to vacancy adjacent or nearest neighboring Fe and magnetic
moment on Fe2 and Fe4 were same, (~3.8+) corresponding to oxygen vacancy second nearest
neighboring Fe (Figure 7.3(c)). After oxygen vacancy moved to the site between Fe1 and Fe2 (Îś =
1), oxidation state on Fe2 increased to 4.1+ due to sqare pyramidal coordination and second nearest

88

neighboring Fe3 and Fe4 have oxidation state of 3.8+ due to octahedral coordination. So, the
oxidation state on Fe2 and Fe3 interchanged due to oxygen vacancy migration. Same as cubic
LSF55, in cubic SFO (which has a large polaron volume Figure 4.2(d)) the charge on the second
nearest neighboring Fe4 remained constant during oxygen migration (Figure 7.3(c)). Overall,
Figure 7.3(c) also shows that the extent of magnetic moment/Fe oxidation change on the oxygenvacancy-neighboring Fe

B)

was not as large as it was for LFO

B),

due to the

distributed nature of the polaron.

Figure 7.3. Oxygen migration in cubic SrFeO3-δ. The pink atoms in Figure (a) highlight the oxygen
migration path, Figure (b) shows the oxygen migration barrier along this path for two different δ,
and Figure (c) shows the charge change on the neighboring Fe atoms as a function of oxygen
atom coordinate on its path.

89

Figure 7.4. Oxygen migration in tetragonal SrFeO3-δ. (Left) non-red smaller atoms highlights the
oxygen migration paths. Paths are defined with Wyckoff positions of oxygen atoms in the perfect
lattice; P1: 4c  16k  8h  16k  4c; P2: 4c  16k  16m  16k  4c; P3: 4c  16k 
4c. (Right) oxygen migration barrier at different coordinates (Îś) of the path.

Figure 7.5. Oxygen migration in orthorhombic SrFeO3-δ. (Left) non-red smaller atoms highlight
the oxygen migration paths. Paths are defined with Wyckoff positions of oxygen atoms in the
perfect lattice; P1: 2b  16r  16r  2b (c-direction); P2: 2b  16r  4j  16r  2b; P3: 2b
 16r  16r  2b (a-direction). (Right) oxygen migration barrier at different coordinates (Μ) of
the path.

Due to the presence of multiple oxygen vacancy migration paths in tetragonal and orthorhombic
strontium ferrite (multiple migration paths are designated in Figure 7.4-7.5 by the Wyckoff
positions of the oxygen atoms in perfect lattice along each path), the minimum energy path for
oxygen ion migration was used to estimate the vacancy migration barrier. Figure 7.4 and 7.5 show
that for tetragonal and orthorhombic SFO, the migration barriers are 0.71 and 1.17 eV, respectively.

90

Figure 7.4 shows in tetragonal SFO, three different migration paths: P1, P2, and P3 have almost
the same oxygen vacancy migration barrier of ~0.7 eV. Therefore, migration barrier in tetragonal
SFO is isotropic in nature, independent of direction of migration. Whereas, Figure 7.5 shows
oxygen vacancy migration barrier in orthorhombic SFO is anisotropic in nature.

Figure 7.6. Oxygen vacancy (a) concentration (b) ion diffusivity, and (c) ion conductivity in
different LSF compositions as a function of temperature in air. In Figure legends ‘cubic SFO’
represent cubic SrFeO3-δ data. ‘other SFO’ represent combined data for tetragonal, orthorhombic,
and brownmillerite phase in Figure 7.6(a) and tetragonal and orthorhombic phases in Figure
7.6(b,c). ‘LSF55(combined)’ represents the LSF55 data with phase change from rhombohedral to
cubic. ‘rhomb LSF55’ in Figure 7.6(c) represent the computed oxygen ionic conductivity data for
rhombohedral LSF55 phase.

91

Figure 7.6 (cont’d)

7.2

Oxygen Ionic Conductivity in La1-xSrxFeO3-δ as a Function of Temperature (in Air)

Using oxygen vacancy site fractions calculated in Chapter 6, oxygen vacancy concentration was
calculated for different LSF compositions, with Equation 1.5. Oxygen ion diffusion coefficient or
diffusivity in different LSF compositions and phases were calculated from the migration barrier
(discussed in Section 7.1) using Equation 1.8. The overall oxygen ion conductivity was calculated
using Equation 1.1. The data was plotted as a function of temperature in air in Figure 7.6. As the
phase diagram (Figure 6.3(c)) have shown several phase transformations. This was included in the
calculation for “LSF55 (combined)”, in which the rhombohedral to cubic phase transformation
occurs at 250Âą50oC.65 To illustrate the impact of phase change on ionic conductivity, individual
phases, such as “cubic SFO” and “rhomb LSF55” were calculated (ignoring the phase change) for
comparison. “other SFO” means tetragonal, orthorhombic and brownmillerite SFO phases were
considered in their specific phase stable temperature range. The ionic conductivity for LFO were
too low to show in Figure 7.6(c).

92

Figure 7.6(a) shows oxygen vacancy concentration in different phases of LSF. Due to very
high oxygen vacancy formation energy, ~3.5 eV in LFO (Figure 5.3), vacancy concentration in it
was relatively low. Oxygen vacancy concentration data for LSF55 shows a kink around 250oC due
to rhombohedral to cubic phase transition in bulk structure.65 Cubic SFO remains almost constant
due to low oxygen vacancy formation energy in this material. Whereas, oxygen vacancy
concentration in other (tetragonal, orthorhombic, brownmillerite) SFO phases fluctuates with
temperature, due to oxygen vacancy ordered phase transition which causes loss in mobile oxygen
vacancy concentration from vacancy site ordering. Once brownmillerite SFO formed at high
temperature, the vacancy concentration in it also drastically reduced due to high oxygen vacancy
formation energy, ~3 eV. As a result, oxygen vacancy concentration in brownmillerite SFO
becomes comparable to vacancy concentration in LFO. Maximum oxygen vacancy concentration
can be observed in orthorhombic SFO at ~ 650oC in the air.
The oxygen diffusivity in LFO is lower than the other LSF compositions. Rhombohedral
LSF55 shows the highest oxygen ionic diffusivity in LSF family at low temperature but it drops
after 250oC, due to rhombohedral to cubic phase transition and oxygen vacancy migration barrier
in cubic LSF55 was higher than rhombohedral LSF55 (Section 7.1.2). Diffusivity of oxygen ion
in tetragonal and orthorhombic SFO was lower than cubic SFO due to higher migration barrier in
previous phases as reported in Section 7.1.3.

Figure 7.6(c) shows the ionic conductivity comparison in different LSF compositions. It clearly
captured LSF55 possesses the highest ionic conductivity at or above 800oC in the air. The dotted
line represent hypothetical rhombohedral LSF55 phase (as rhombohedral LSF55 transform to
cubic at ~250oC in air). The calculated ionic conductivity for LSF55 at 950oC was 0.3977 S/cm
was comparable to experimentally measured ionic conductivity (0.5 S/cm) at the same
93

temperature.62 As LSF64 has better stability for rhombohedral phase, this may be the reason for
the broader application of LSF64 as the best MIEC candidate in LSF family.166 If the rhombohedral
LSF55 can be stabilized at a higher temperature, it will be even more desirable than LSF64, due
to its higher oxygen vacancy concentration. On the other-hand, SFO will lose its ionic conductivity
at high temperature due to continuous phase change, starting from the tetragonal phase. Though it
was previously shown in Figure 6.3(b-c) orthorhombic SFO possess highest oxygen vacancy
concentration within the SOFC operating window (0 – 1000oC) in the air but due to high oxygen
ion migration barrier, it is not best candidate as oxygen ion conductor. Due to low oxygen vacancy
formation energy, cubic SFO possess the highest oxygen ionic conductivity at low temperature
(below 300oC). Therefore, stabilizing cubic SFO (avoiding phase transformation) to higher
temperature will lead to higher oxygen ionic conductivity.

94

8.

Anisotropic Chemical Strain Produced by Non-Stoichiometric Cubic CeO2-δ

It can be noticed from Figure 2.5(a) each oxygen atom/ oxygen vacancy site (VO•• ) in Ceria is
coordinated with four Cerium atoms. As discussed in Section 2.6 on neutral oxygen vacancy
formation two electrons are generated, it can be distributed to four adjacent Ce-atoms (symmetric
distribution) or can remain localized to two Ce atoms (asymmetric distribution).
8.1

Charge Disproportionation in Non-stoichiometric Ceria
The asymmetric scenario provided energetically favorable minima. The details of

calculation results are tabulated in Table 8.1. The non-uniform distortion was energetically
favorable over uniform distortion by 0.55 eV. Without initial perturbation, four Ce-atoms will
show a 3.5+ charge due to an equal share of electrons (GGA calculations symmetrically distribute
the electron cloud). Whereas in energetically favored case, two perturbed-Ce maintained 4+ charge
and the other two Ce became 3+ (Figure 8.1(a)). According to the bond distance calculation,
oxygen atoms (O-Ce-VO•• ) opposite to VO•• of Ce, moved closer to the Ce in <111>.
Table 8.1. Comparing calculated energy and bond lengths for different Vo-formation cases.
Ceria Structure
(2x2x2)
Perfect lattice
With
VO••

symmetric

With asymmetric
VO••

0K
Formation
Energy (eV)
-783.63
-775.20

-775.75

Ce – O bond
distance (Å)

Charge
on Ce
4+
3.5+

3+

4+

2.38
3x2.38,
3x2.37,
1x2.29
2x2.44,
3x2.42,
1x2.43,
1x2.33
2x2.31,
2x2.35,
1x2.33,
1x2.34,
1x2.27
95

Ce – VO••
bond distance
(Å)
2.53

VO•• /O – Ce –
O distance
(Å)
4.76
4.82

2.52

4.85

2.56

4.82

This observation is opposite to the hypothesis of Korobko et al.75 All Ce ions moved away from
the VO•• -site along <111>, Ce4+ moved further than Ce3+. Four O-atoms along the O-Ce-VO••
direction moved away from VO•• , O-Ce3+-VO•• moved further than O-Ce4+-VO•• (Table 8.1). The
complicated local distortion due to the VO•• formation is shown in Figure 8.1(b). It is the core of
lattice distortion from the 2x2x2 supercell after oxygen vacancy formation. The complicated
distortion was partially hypothesized before136,137 but here the complete description along with its’
effect on chemical strain is described. All the (seven) Ce4+-O2- bonds were shortened (blue). Six
Ce3+-O2- bonds were elongated (red) except the O-Ce3+ bond, with O-atom diagonally opposite to
the VO•• -site (blue). Green-O came closer to the VO•• -site and other red-O moved away from VO•• . OO distance in a perfect ceria lattice was 2.75 Å. Four O-atoms (green) connected to both Ce3+ and
Ce4+ came closer to VO•• and was 2.62 Å from the VO•• -site. O-atom (green) connected to two Ce4+
was 2.47 Å away from the VO•• -site. O-atom (red) connected to two Ce3+ moved away from VO•• ,
was 2.76 Å away from the VO•• -site.

Figure 8.1. Charge distribution on Ce and displacement of O-atoms around neutral oxygen
vacancy (𝑉𝑂•• ). Red-O moved away from oxygen vacancy site. Green-O came closer to oxygen
vacancy site. Red Ce-O bond signify elongation with respect to bond length in perfect structure
and blue Ce-O bond signify contraction.

96

8.2

PDOS of Perfect and Oxygen Vacancy Adjacent Ce

Partial density of states (PDOS) was calculated for Ce and O atoms in the perfect CeO2 lattice and
shown in Figure 8.2. Figure 8.3(a-c) show the Ce-5d and 4f orbitals in vacancy containing lattice
for two different scenarios. Analysis of the PDOS for Ce-5d and 4f orbitals before and after the
formation of oxygen vacancy has shown the charge localization in 4f orbital after VO•• formation.
Electron localization in Ce-4f orbital observed in Figure 8.3(b) is consistent with previous
computational studies.69,172 Whereas, Figure 8.3(a) shows the Ce4+ state which is comparable to
PDOS of Ce4+ present in the perfect lattice (Figure 8.2). PDOS plot for Ce3.5+ in Figure 8.3(c)
shows electron density in the 4f orbital of Ce is clearly in between of PDOS for Ce3+ and Ce4+.

Figure 8.2. Partial density of states plot for Ce-5s, 5p, 5d, 4f and O-2p orbitals in pure CeO2

97

(b)

(a)

(c)

Figure 8.3.(a) PDOS for Ce4+-5d and Ce4+-4f orbital for Ce4+ adjacent to 𝑉𝑂•• , after 𝑉𝑂•• formation
in asymmetric case. Orbitals show same shape as Ce4+ in perfect CeO2.(b) PDOS for Ce3+-5d and
Ce3+-4f orbital for Ce3+ adjacent to 𝑉𝑂•• , after 𝑉𝑂•• formation in asymmetric case. Figure shows
localization of excess electron in 4f orbital. .(c) PDOS for Ce3.5+-5d and Ce3.5+-4f orbital for Ce3.5+
adjacent to 𝑉𝑂•• , after 𝑉𝑂•• formation in symmetric case. Figure shows localization of excess
electron in 4f orbital. This scenario is computational artifact as the structure was energetically
unfavorable.

98

8.3

Calculating Anisotropic Chemical Strain
It was interesting to notice how the non-uniform local deformation is impacting the overall

chemical strain. Method to calculate the anisotropic chemical strain coefficient tensor is described
𝑓

in Section 2.6. 𝐺 tensor is calculated from the slope of 𝑑𝐸V•• ,𝜀 versus 𝜀 as shown in Figure 8.4.
O

Figure 8.4. The oxygen vacancy formation energy varies with an applied strain in different
directions. The slope of the strain versus energy curve gives the elastic dipole tensor mapped with
applied strain direction.
−9.633
𝐺 = [ 0.001
0.001

0.001
0.001
−9.633 1.734 ]
1.734 −9.633

Energy for the short and long rangee interaction was calculated using Equation 2.22 and 2.23
respectively. 𝑉𝑈 for this calculation is 41.50Å3 , volume per formula unit of the perfect CeO2. Due
to cubic symmetry of the perfect lattice components of the elastic stiffness tensor, ℂ can be defined
as 𝐶11 = 343 GPa, 𝐶22 = 103 GPa, and 𝐶44 = 54 GPa, was comparable to previous calculations.173
Calculated Young’s modulus of 198 GPa was also comparable to previous experimental study.174
99

Calculated 𝜀𝐶 at δ = 0.03 and 𝛼𝐶 is provided below. The mean chemical expansion coefficient
𝛼𝐶,𝑀 = 0.067 was comparable to experimentally obtained average 𝛼𝐶 .40,175

0.002
𝜀𝐶 = [0.000
0.000

0.000
0.002
−0.039

0.000
0.067
−0.039], 𝛼𝐶 = [0.000
0.002
0.000

0.000
0.000
0.067 −0.124]
−0.124 0.067

The 𝛼𝐶 calculated above was not able to show the anisotropic effect in principle directions.
On diagonalization of the 𝛼𝐶 , one can obtain the total effect of the chemical strain projected in
principle directions. Therefore, the Eigen values and Eigen vectors were calculated. Where, Eigen
vectors provides the rotation matrix, R and Eigen values show components of the chemical
expansion coefficient mapped to the principle directions, 𝛼𝐶,𝑃 .
0 0
𝑅 = [1̅ 1
1 1

1
0.191
0], 𝛼𝐶,𝑃 = [0.000
0
0.000

0.000
−0.057
0.000

0.000
0.000]
0.067

On close observation of the rotation matrix (or the Eigen vectors), one can notice that the x and yaxis of the diagonalized tensor was mapped to [01̅1] and [011] direction respectively. This is the
direction between Ce3+ and Ce4+ (Figure 8.1(a)). Therefore, orienting the lattice in <110> one can
maximize the anisotropic chemical strain effect. After diagonalizing the strain tensor, the
anisotropic effect is clearly observed. The maximum possible chemical expansion coefficient was
0.191, and in perpendicular direction, there was compressive strain effect with 𝛼𝐶,𝑦𝑦 = -0.057.
Principle positive chemical expansion coefficient, 𝛼𝐶,𝑃𝑃 for ceria was calculated as 0.129. The
average 𝜀𝐶 with δ, observed in the experiment is comparable to the calculated average 𝜀𝐶 with δ
in this study (Figure 8.5(a)). 𝛼𝐶,𝑃𝑃 provided the upper bound for all the experimental data, whereas
𝛼𝐶,𝑀 is comparable to previous experimental observations (Figure 8.5(a)).

100

Figure 8.5.(a) Average expansion of single crystal ceria with increasing oxygen nonstoichiometry.
(b) possible orientations of Ce3+ and Ce4+ atoms around oxygen vacancy.
There are six possible directions in nonstoichiometric ceria for local orientation in <110>
among Ce3+ and Ce4+ ions around the VO•• -site (Figure 8.5(b)) but without any applied electric field
they are randomly oriented giving the average effect of the chemical strain (𝛼𝐶,𝑀 = 0.067).
Assuming an electric field applied in <110> direction of the Figure 8.5(b), it can be noticed only
in Figure 8.5(b)(i-ii) Ce4+-Ce3+ are oriented in preferred direction of applied bias but other Ce
atoms in Figure 8.5(b)(iii-vi) are randomly oriented. On application of the electric field electrons
got transferred from Ce3+ to Ce4+ to align the polaron in the preferred orientation of the applied
bias so there are less degenerated states. The resultant effect of the preferred orientation produce
the amplified chemical strain in the material. If the crystal orientation is known then chemical
strain due to electric field, 𝐸𝛼𝐶 can be calculated as 𝐸𝛼𝐶 = 𝑂𝑟𝑖𝑒𝑛𝑡𝑒𝑑(𝛼𝐶 ) − 𝑡𝑟(𝛼𝐶 ).

101

9.

Conclusion and Future Work
9.1

Thesis Summary and Conclusion

This thesis performed DFT-based investigations on the oxygen vacancy formation, interaction,
transport, and its’ chemical strain effect in perovskite LSF, and fluorite CeO2-δ, as two important
members of oxide families for SOFC applications. Several computational challenges were solved
and many insights were developed from these calculations.
To elucidate the effect of oxygen vacancy polaron shape and size on the vacancy interaction
in LSF, three different LSF compositions were computationally studied with GGA+U. HubbardU (= 3) parameter was calibrated in order to predict the magnetic moment on Fe as a function of
La/Sr ratio and oxygen non-stoichiometry. The magnetic moment on Fe was utilized to interpret
the formal charge on Fe. The oxidation state change on Fe upon oxygen vacancy generation
outlined the shape and size of the oxygen polarons. When an oxygen vacancy is formed in a
perovskite lattice, the two neighboring Fe-O coordinated octahedra collapse into two square
pyramids. It was found that in LaFeO3 the excess electrons produced during oxygen vacancy
formation were localized at the square pyramidally (Fe-O5) coordinated Fe atoms adjacent to the
oxygen vacancy, resulting in a small polaron. In contrast, in SrFeO3 the excess electrons produced
during oxygen vacancy formation were transferred to the octahedrally (Fe-O6) coordinated secondnearest-neighbor Fe atoms, resulting in a pancake-shaped large polarons. This difference is well
explained by the electron occupancy and crystal-field-splitting-induced differences between the
Fe 3d-orbitals of the square pyramidally coordinated, oxygen-vacancy-adjacent Fe atoms and the
octahedrally-coordinated, oxygen-vacancy-second-nearest-neighbor Fe atoms. For La0.5Sr0.5FeO3δ,

the polaron size in the rhombohedral phase was smaller than that in cubic phase, due to the tilt

of Fe-O6 octahedra. In all cases, larger oxygen vacancy polarons resulted in higher oxygen vacancy
102

interactions. In fact, in multiple phases, the oxygen vacancy formation energy started to increase
with oxygen non-stoichiometry beyond a critical value when the polarons began to overlap. There
has been a long-standing debate in the literature on the mixed charge states of Fe in the SrFeO3-δ.
The magnetic moment interpreted charge state correlations presented here clearly show that square
pyramidal Fe always remains +4 in SrFeO3-δ. The charges of Fe in the distorted octahedra in
SrFeO3-δ vary from +4 to +3 depending on δ, and they are easier to reduce than those in square
pyramidal coordination, as explained with crystal field theory. As a result, in all strontium ferrite
structures 1) oxygen vacancies are always formed at the site shared by the two highest charged
octahedral Fe (except brownmillerite) and 2) the oxygen vacancy formation energy increases after
each phase transformation.
Overall, the calculated polaron shape changes and progressively increased polaron sizes
with increasing Sr-doping, resulted in increasing oxygen vacancy formation energies in La1xSrxFeO3-δ.

This was consistent with experimental results showing that Sr-rich LSF and highly

oxygen deficient compositions are more prone to oxygen-vacancy-ordering-induced phase
transformations, while Sr-poor and oxygen rich LSF compositions are not. Since oxygen vacancy
induced phase transformations cause a decrease in the mobile oxygen vacancy site fraction (X),
both δ and X were also predicted as a function of temperature and oxygen partial pressure, for
multiple LSF compositions and phases using a combined thermodynamics and DFT approach.
The oxygen vacancy migration barrier did not change with the oxygen nonstoichiometry
in rhombohedral LSF55 and orthorhombic LFO, due to smaller polaron volume. With stronger
oxygen vacancy interaction, SFO demonstrated increasing oxygen vacancy migration barrier with
phase change (cubic < tetragonal <orthorhombic). LSF55 also showed increased oxygen vacancy
migration barrier due to rhombohedral to cubic phase change. Finally, computationally predicted

103

ionic conductivity showed rhombohedral LSF has higher conductivity than cubic LSF at the same
SOFC operating conditions. Cubic SFO is promising at low temperature, as it possess the highest
ionic conductivity than other phases until its cubic to tetragonal phase transition at ~250oC.
Comparing all the phases and compositions of LSF, the overall effect of phase change on the
oxygen vacancy migration barrier is not as substantial (differ by ~0.7eV), as it was on oxygen
vacancy formation energy (differ by ~3eV). Therefore increasing oxygen vacancy site fraction can
be more beneficial in improving the conductivity.
Computationally predicted anisotropic chemical strain in cubic ceria explained the reason
behind observed anomalous strain effect in this material on application of electric field. The
presence of highly localized 4f valance orbital in Ce causes charge disproportionation on the
formation of neutral oxygen vacancy, producing anisotropic chemical strain in crystal with cubic
symmetry. Calculated average chemical expansion was comparable to experimental observation.
Without application of external bias, strain field polaron remain randomly oriented giving average
effect of anisotropic chemical strain. On the application of electric field, Ce3+-Ce4+ bonds orient
perpendicular to the thin film enhancing chemical strain effect.

9.2

Future Work

Cubic SFO has shown a good ionic conductivity at relatively low temperature due to very low
oxygen vacancy formation energy (0.4 eV) at the dilute condition.176 But due to vacancy ordering
induced phase change to tetragonal structure, it loses its ionic conductivity even at ~250oC. To
resist the vacancy ordering induced phase transition, a Ruddlesden-Popper (RP)177 type strontium
ferrite structure can be investigated in future.178 As shown in Figure 10.1 the RP type Sr3Fe2O7
(RP-SFO) has layered FeO6 octahedral structure, where Fe maintains its 4+ charge state.179 Further,
it was discussed in Chapter 5, that it was easier to form oxygen vacancy adjacent to Fe 4+ and
104

prevent vacancy ordered phase transition. However, its diffusion mechanism is largely unknown.
Maljuk et al.180 have grown high purity RP-SFO single crystal, experimentally. Its potential as a
SOFC cathode material is worth to investigate.

Figure 9.1. Ruddlesden-Popper type Sr3Fe2O7 lattice with single oxygen vacancy (black atom) at
the center. Blue, brown, and red atoms are Sr, Fe, and O respectively. Sr3Fe2O7 can be a promising
compound for low-temperature SOFC cathode.

105

Figure 9.2. Crystal structure of 12.5% Gd/Sm/Pr (purple atoms) doped ceria (GDC/SDC/PDC).
Here, Ce and O atoms are colored in yellow and red respectively. Gd/Sm/Pr atoms are randomly
distributed in the lattice.

Gd, Sm, or Pr doped ceria is commonly used in place of pure ceria due to increased oxygen
vacancy concentration in it via doping, under same operating environment. Chapter 8 elucidates
how the presence of dilute oxygen vacancy produces anisotropic chemical strain in the the cubic
lattice and described the complicated lattice distortion around oxygen vacancy. Further studies are
required to understand the effect of the highly doped ceria and how the presence of foreign element
around oxygen vacancy influences the shape and size of the polaron. In addition, studies are
required to understand the effect of oxygen vacancy interaction in doped ceria with higher oxygen
vacancy concentration.

106

APPENDIX

107

APPENDIX

Numerical Method to Find Temperature Corresponding to an Oxygen Nonstoichiometry for
Certain O2 Partial Pressure (MATLAB code) to Solve Equation 2.17.
% Temperature from nonstoichiometry
TC = 0:0.1:2000; % Temperature range in deg C
T = TC+273.15; % T in K
p = 0.21; % atm
R = 8.314e-3; % kJ/mol-K
cF = 96.485; % 1eV = 96.485 kJ/mol --> conversion Factor
Nv = 3; % # of available sites/ formula unit
x = 0.0625; % Nonstoichiometry
Evac_DFT = 0.56; % Vacancy formation energy
%
for i = 1:length(T)
mu_O2_corr = O2_correction(T(i),p)/cF;
%
G_vac(i) = Evac_DFT + 0.5 * mu_O2_corr;
%
A_fac = exp(-G_vac(i)*cF/(R*T(i)));
nonstoi(i) = Nv*A_fac/(1+A_fac);
%
if nonstoi(i) >= x
vac_form = E_vac(i);
TempK=T(i);
delta = nonstoi(i);
break;
end
end
%
TempC = TempK-273.15% Temperature in deg C
Function: O2_correction
function mu_O2_corr = O2_correction(T,p)
% This function returns connection energy in eV for oxygen MOLECULE
% between room temperature(298.15K) at 1 bar (0.987 atm) and
% any other temperature (K) and pressure (atm)
%
Tr = 298.15; % K
Sr = 205.152e-3; % kJ/mol-K
p0 = 0.986923267; % atm
R = 8.314e-3; % kJ/mol-K
cF = 96.485; % 1eV = 96.485 kJ/mol --> conversion Factor
108

%
delG_O2 = -0.28*2; % Ref # 123 Average for O2
%
if T < 700
% for Temp range 100 to 700 K
A = 31.32234; % J/(mol-K) ******
http://webbook.nist.gov/cgi/cbook.cgi?ID=C7782447&Units=SI&Mask=1#Thermo-Gas
*******
B = -20.23531; % J/(mol-K^2)
C = 57.86644; % J/(mol-K^3)
D = -36.50624; % J/(mol-K^4)
E = -0.007374; % J-K/mol
F = -8.903471; % kJ/mol
G = 246.7945; % J/(mol-K)
else
% for Temp range 700 to 2000 K
A = 30.03235; % J/(mol-K) ******
http://webbook.nist.gov/cgi/cbook.cgi?ID=C7782447&Units=SI&Mask=1#Thermo-Gas
*******
B = 8.772972; % J/(mol-K^2)
C = -3.988133; % J/(mol-K^3)
D = 0.788313; % J/(mol-K^4)
E = -0.741599; % J-K/mol
F = -11.32468; % kJ/mol
G = 236.1663; % J/(mol-K)
end
%
t = T/1000; % T in K
% del_H = H(pr,T) - H(pr,Tr)
del_H = A*t +(B/2)*t^2 +(C/3)*t^3 +(D/4)*t^4 -E/t+F; % kJ/mol
%
% del_S = S(pr,T) - S(pr,Tr)
S = A*log(t) +B*t +(C/2)*t^2 +(D/3)*t^3-E/(2*t^2)+G; % J/mol-K
S = S/1000; % kJ/mol-K
del_S = S - Sr;
%
mu_O2_corr = delG_O2*cF + del_H - T*del_S - (T-Tr)*Sr + R*T*log(p/p0); % as p0 = 1 bar
end

109

BIBLIOGRAPHY

110

BIBLIOGRAPHY

(1) Yadav, A. K.; Singh, R. K.; Singh, P. Fabrication of Lanthanum Ferrite Based Liquefied
Petroleum Gas Sensor. Sens. Actuators B Chem. 2016, 229, 25–30.
(2) Brisotto, M.; Cernuschi, F.; Drago, F.; Lenardi, C.; Rosa, P.; Meneghini, C.; Merlini, M.;
Rinaldi, C. High Temperature Stability of Ba0.5Sr0.5Co0.8Fe0.2O3−δ and
La0.6Sr0.4Co1−yFeyO3−δ Oxygen Separation Perovskite Membranes. J. Eur. Ceram. Soc.
2016, 36 (7), 1679–1690.
(3) Magyari-KĂśpe, B.; Park, S. G.; Lee, H.-D.; Nishi, Y. First Principles Calculations of Oxygen
Vacancy-Ordering Effects in Resistance Change Memory Materials Incorporating Binary
Transition Metal Oxides. J. Mater. Sci. 2012, 47 (21), 7498–7514.
(4) Baumann, F. S.; Fleig, J.; Cristiani, G.; Stuhlhofer, B.; Habermeier, H.-U.; Maier, J.
Quantitative Comparison of Mixed Conducting SOFC Cathode Materials by Means of Thin
Film Model Electrodes. J. Electrochem. Soc. 2007, 154 (9), B931–B941.
(5) Nicholas, J.; Highlights from the 2013 National Science Foundation Solid Oxide Fuel Cell
Promise, Progress, and Priorities (SOFC-PPP) Workshop. Electrochem. Soc. Interface 2013,
No. Winter 2013, 49.
(6) Yang, Q.; Burye, T. E.; Lunt, R. R.; Nicholas, J. D. In Situ Oxygen Surface Exchange
Coefficient Measurements on Lanthanum Strontium Ferrite Thin Films via the Curvature
Relaxation Method. Solid State Ionics. 2013, 249–250 (0), 123–128.
(7) Burye, T. E.; Tang, H.; Nicholas, J. D. The Effect of Precursor Solution Desiccation or NanoCeria Pre-Infiltration on Various La0.6Sr0.4FeyCo1-YO3-x Infiltrate Compositions. J.
Electrochem. Soc. 2016, 163 (9), F1017–F1022.
(8) Celorrio, V.; Bradley, K.; Weber, O. J.; Hall, S. R.; FermĂ­n, D. J. Photoelectrochemical
Properties of LaFeO3 Nanoparticles. ChemElectroChem 2014, 1 (10), 1667–1671.
(9) Wachsman, E. D.; Lee, K. T. Lowering the Temperature of Solid Oxide Fuel Cells. Science
2011, 334 (6058), 935–939.
(10) Mogensen, M.; Nigel, S.; Geoff, T. Physical, Chemical and Electrochemical Properties of
Pure and Doped Ceria. Solid State Ionics. 2000, 129 (1–4), 63–94.
(11) Yoo, J. Determination of the Equilibrium Oxygen Non-Stoichiometry and the Electrical
Conductivity of La0.5Sr0.5FeO3-X. Solid State Ionics. 2004, 175 (1–4), 55–58.
(12) Liu, Z.; Ding, D.; Liu, M.; Ding, X.; Chen, D.; Li, X.; Xia, C.; Liu, M. High-Performance,
Ceria-Based Solid Oxide Fuel Cells Fabricated at Low Temperatures. J. Power Sources 2013,
241, 454–459.

111

(13) Dalslet, B.; Blennow, P.; Hendriksen, P. V.; Bonanos, N.; Lybye, D.; Mogensen, M.
Assessment of Doped Ceria as Electrolyte. J. Solid State Electrochem. 2006, 10 (8), 547–561.
(14) Ruiz de Larramendi, I.; Ortiz-Vitoriano, N.; Acebedo, B.; Jimenez de Aberasturi, D.; Gil de
Muro, I.; Arango, A.; RodrĂ­guez-CastellĂłn, E.; Ruiz de Larramendi, J. I.; Rojo, T. Pr-Doped
Ceria Nanoparticles as Intermediate Temperature Ionic Conductors. Int. J. Hydrog. Energy
2011, 36 (17), 10981–10990.
(15) Spiridigliozzi, L.; Biesuz, M.; Dell’Agli, G.; Bartolomeo, E. D.; Zurlo, F.; Sglavo, V. M.
Microstructural and Electrical Investigation of Flash-Sintered Gd/Sm-Doped Ceria. J. Mater.
Sci. 2017, 52 (12), 7479–7488.
(16) Liu, T.; Zhang, X.; Wang, X.; Yu, J.; Li, L. A Review of Zirconia-Based Solid Electrolytes.
Ionics 2016, 22 (12), 2249–2262.
(17) Mizusaki, J.; Sasamoto, T.; Cannon, W. R.; Bowen, H. K. Electronic Conductivity, Seebeck
Coefficient, and Defect Structure of La1-XSrxFeO3 (X=0.l, 0.25). J. Am. Ceram. Soc. 1983,
66 (4), 247–252.
(18) Berenov, A. V.; Atkinson, A.; Kilner, J. A.; Bucher, E.; Sitte, W. Oxygen Tracer Diffusion
and Surface Exchange Kinetics in La0.6Sr0.4CoO3−δ. Solid State Ionics. 2010, 181 (17–18),
819–826.
(19) Grundy, A. N.; Hallstedt, B.; Gauckler, L. J. Assessment of the La–Sr–Mn–O System.
Calphad 2004, 28 (2), 191–201.
(20) De Souza, R. A.; Kilner, J. A. Oxygen Transport in La1−xSrxMn1−yCoyO3±δ Perovskites:
Part I. Oxygen Tracer Diffusion. Solid State Ionics. 1998, 106 (3–4), 175–187.
(21) Bansal, N. P.; Zhong, Z. Synthesis of Sm0.5Sr0.5CoO3-x and La0.6Sr0.4CoO3-x
Nanopowders by Solution Combustion Process. In Novel Processing of Ceramics and
Composites; Bansal, N. P., Singh, J. P., Smay, J. E., Ohji, T., Eds.; John Wiley & Sons, Inc.:
Hoboken, NJ, USA, 2011; pp 21–32.
(22) Ahmed, B.; Lee, S.-B.; Song, R.-H.; Lee, J.-W.; Lim, T.-H.; Park, S.-J. Development of
Novel LSM/GDC Composite and Electrochemical Characterization of LSM/GDC Based
Cathode-Supported Direct Carbon Fuel Cells. J. Solid State Electrochem. 2014, 18 (2), 435–
443.
(23) Liu, M.; Ding, D.; Blinn, K.; Li, X.; Nie, L.; Liu, M. Enhanced Performance of LSCF Cathode
Through Surface Modification. Int. J. Hydrog. Energy 2012, 37 (10), 8613–8620.
(24) Ritzmann, A. M.; MuĂąoz-GarcĂ­a, A. B.; Pavone, M.; Keith, J. A.; Carter, E. A. Ab Initio
DFT+U Analysis of Oxygen Vacancy Formation and Migration in La1-XSrxFeO3-δ = 0,
0.25, 0.50). Chem. Mater. 2013, 25 (15), 3011–3019.
(25) Cassir, M.; Gourba, E. Reduction in the Operating Temperature of Solid Oxide Fuel Cells—
Potential Use in Transport Applications. Ann. Chim. Sci. Matér. 2001, 26 (4), 49–58.
112

(26) Steele, B. C. H.; Heinzel, A. Materials for Fuel-Cell Technologies. Nature 2001, 414 (6861),
345–352.
(27) Das, T.; Nicholas, J. D.; Qi, Y. Long-Range Charge Transfer and Oxygen Vacancy
Interactions in Strontium Ferrite. J. Mater. Chem. A 2017, 5 (9), 4493–4506.
(28) Lee, Y.-L.; Morgan, D. Ab Initio and Empirical Defect Modeling of LaMnO3¹δ for Solid
Oxide Fuel Cell Cathodes. Phys Chem Chem Phys 2012, 14 (1), 290–302.
(29) Kuklja, M. M.; Mastrikov, Y. A.; Jansang, B.; Kotomin, E. A. The Intrinsic Defects,
Disordering, and Structural Stability of BaxSr1–xCoyFe1–yO3−δ Perovskite Solid Solutions.
J. Phys. Chem. C 2012, 116 (35), 18605–18611.
(30) Kuklja, M. M.; Kotomin, E. A.; Merkle, R.; Mastrikov, Y. A.; Maier, J. Combined
Theoretical and Experimental Analysis of Processes Determining Cathode Performance in
Solid Oxide Fuel Cells. Phys. Chem. Chem. Phys. 2013, 15 (15), 5443.
(31) Das, T.; Nicholas, J. D.; Qi, Y. First-Principles Studies of Oxygen Vacancy Interactions and
Their Impact on Oxygen Migration in Lanthanum Strontium Ferrite. ECS Trans. 2017, 78
(1), 2807–2814.
(32) Hodges, J. P.; Short, S.; Jorgensen, J. D.; Xiong, X.; Dabrowski, B.; Mini, S. M.; Kimball, C.
W. Evolution of Oxygen-Vacancy Ordered Crystal Structures in the Perovskite Series
SrnFenO3n−1 (N=2, 4, 8, and ∞), and the Relationship to Electronic and Magnetic Properties.
J. Solid State Chem. 2000, 151 (2), 190–209.
(33) Schmidt, M.; Hofmann, M.; Campbell, S. J. Magnetic Structure of Strontium Ferrite
Sr4Fe4O11. J. Phys. Condens. Matter 2003, 15 (50), 8691.
(34) Vidya, R.; Ravindran, P.; FjellvĂĽg, H.; Kjekshus, A. Spin- and Charge-Ordering in OxygenVacancy-Ordered Mixed-Valence Sr4Fe4O11. Phys. Rev. B 2006, 74 (5), 054422.
(35) Adler, P. Comment on ``Spin and Charge Ordering in Oxygen Vacancy Ordered MixedValence Sr4Fe4O11. Phys. Rev. B 2008, 77 (13), 136401.
(36) Ravindran, P.; Vidya, R.; FjellvĂĽg, H.; Kjekshus, A. Reply to ``Comment on `Spin- and
Charge-Ordering in Oxygen-Vacancy-Ordered Mixed-Valence Sr4Fe4O11. Phys. Rev. B
2008, 77 (13), 136402.
(37) Hombo, J.; Matsumoto, Y.; Kawano, T. Electrical Conductivities of SrFeO3−δ and
BaFeO3−δ Perovskites. J. Solid State Chem. 1990, 84 (1), 138–143.
(38) Kharton, V.; Kovalevsky, A.; Tsipis, E.; Viskup, A.; Naumovich, E.; Jurado, J.; Frade, J.
Mixed Conductivity and Stability of A-Site-Deficient Sr(Fe,Ti)O3-δ; Perovskites. J. Solid
State Electrochem. 2002, 7 (1), 30–36.

113

(39) Li, Y.; Kraynis, O.; Kas, J.; Weng, T.-C.; Sokaras, D.; Zacharowicz, R.; Lubomirsky, I.;
Frenkel, A. I. Geometry of Electromechanically Active Structures in Gadolinium - Doped
Cerium Oxides. AIP Adv. 2016, 6 (5), 055320.
(40) Marrocchelli, D.; Bishop, S. R.; Tuller, H. L.; Yildiz, B. Understanding Chemical Expansion
in Non-Stoichiometric Oxides: Ceria and Zirconia Case Studies. Adv. Funct. Mater. 2012, 22
(9), 1958–1965.
(41) Er, D.; Li, J.; Cargnello, M.; Fornasiero, P.; Gorte, R. J.; Shenoy, V. B. A Model to Determine
the Chemical Expansion in Non-Stoichiometric Oxides Based on the Elastic Force Dipole. J.
Electrochem. Soc. 2014, 161 (11), F3060–F3064.
(42) James, C.; Wu, Y.; Sheldon, B.; Qi, Y. Computational Analysis of Coupled Anisotropic
Chemical Expansion in Li2-XMnO3-δ. MRS Adv. 2016, 1 (15), 1037–1042.
(43) Ishigaki, T.; Yamauchi, S.; Kishio, K.; Mizusaki, J.; Fueki, K. Diffusion of Oxide Ion
Vacancies in Perovskite-Type Oxides. J. Solid State Chem. 1988, 73 (1), 179–187.
(44) Takeda, Y.; Kanno, K.; Takada, T.; Yamamoto, O.; Takano, M.; Nakayama, N.; Bando, Y.
Phase Relation in the Oxygen Nonstoichiometric System, SrFeOx (2.5 ≤ x ≤ 3.0). J. Solid
State Chem. 1986, 63 (2), 237–249.
(45) Choi, S.; Yoo, S.; Kim, J.; Park, S.; Jun, A.; Sengodan, S.; Kim, J.; Shin, J.; Jeong, H. Y.;
Choi, Y.; Kim, G.; Liu, M. Highly Efficient and Robust Cathode Materials for LowTemperature Solid Oxide Fuel Cells: PrBa0.5Sr0.5Co2−xFexO5+δ. Sci. Rep. 2013, 3.
(46) Kharton, V. V.; Kovalevsky, A. V.; Patrakeev, M. V.; Tsipis, E. V.; Viskup, A. P.; Kolotygin,
V. A.; Yaremchenko, A. A.; Shaula, A. L.; Kiselev, E. A.; Waerenborgh, J. C. Oxygen
Nonstoichiometry, Mixed Conductivity, and Mössbauer Spectra of Ln0.5A0.5FeO3−δ (Ln =
La−Sm, A = Sr, Ba): Effects of Cation Size. Chem. Mater. 2008, 20 (20), 6457–6467.
(47) Burye,
T.
E.;
Nicholas,
J.
D.
Nano-Ceria
Pre-Infiltration
Improves
La0.6Sr0.4Co0.8Fe0.2O3−x Infiltrated Solid Oxide Fuel Cell Cathode Performance. J.
Power Sources 2015, 300, 402–412.
(48) Pavone, M.; MuĂąoz-GarcĂ­a, A. B.; Ritzmann, A. M.; Carter, E. A. First-Principles Study of
Lanthanum Strontium Manganite: Insights into Electronic Structure and Oxygen Vacancy
Formation. J. Phys. Chem. C 2014, 118 (25), 13346–13356.
(49) MuĂąoz-GarcĂ­a, A. B.; Pavone, M.; Ritzmann, A. M.; Carter, E. A. Oxide Ion Transport in
Sr2Fe1.5Mo0.5O6−δ, a Mixed Ion-Electron Conductor: New Insights From First Principles
Modeling. Phys. Chem. Chem. Phys. 2013, 15 (17), 6250.
(50) MuĂąoz-GarcĂ­a, A. B.; Ritzmann, A. M.; Pavone, M.; Keith, J. A.; Carter, E. A. Oxygen
Transport in Perovskite-Type Solid Oxide Fuel Cell Materials: Insights from Quantum
Mechanics. Acc. Chem. Res. 2014, 47 (11), 3340–3348.

114

(51) Burbano, M.; Norberg, S. T.; Hull, S.; Eriksson, S. G.; Marrocchelli, D.; Madden, P. A.;
Watson, G. W. Oxygen Vacancy Ordering and the Conductivity Maximum in Y2O3-Doped
CeO2. Chem. Mater. 2012, 24 (1), 222–229.
(52) Kittel, C. Introduction to Solid State Physics, 8 edition.; Wiley: Hoboken, NJ, 2004.
(53) Ritzmann, A. M.; Pavone, M.; MuĂąoz-GarcĂ­a, A. B.; Keith, J. A.; Carter, E. A. Ab Initio
DFT+U Analysis of Oxygen Transport in LaCoO3: The Effect of Co3+ Magnetic States. J
Mater Chem A 2014, 2 (21), 8060–8074.
(54) Rothschild, A.; Menesklou, W.; Tuller, H. L.; Ivers-TiffĂŠe, E. Electronic Structure, Defect
Chemistry, and Transport Properties of SrTi1-XFexO3-y Solid Solutions. Chem. Mater. 2006,
18 (16), 3651–3659.
(55) Mastrikov, Y. A.; Merkle, R.; Kotomin, E. A.; Kuklja, M. M.; Maier, J. Formation and
Migration of Oxygen Vacancies in La1−xSrxCo1−yFeyO3−δ Perovskites: Insight from Ab
Initio Calculations and Comparison with Ba1−xSrxCo1−yFeyO3−δ. Phys Chem Chem Phys
2013, 15 (3), 911–918.
(56) Gryaznov, D.; Finnis, M. W.; Evarestov, R. A.; Maier, J. Oxygen Vacancy Formation
Energies in Sr-Doped Complex Perovskites: Ab Initio Thermodynamic Study. Solid State
Ionics. 2014, 254, 11–16.
(57) Vashuk, V. V.; Kokhanovskii, L. V.; Yushkevich, I. I. Electrical Conductivity and Oxygen
Stoichiometry of SrFeO3-δ. Neorganicheskie Mater. 2000, 36 (1), 90–96.
(58) Gao, M.; Li, C.-J.; Li, C.-X.; Yang, G.-J.; Fan, S.-Q. Microstructure, Oxygen Stoichiometry
and Electrical Conductivity of Flame-Sprayed Sm0.7Sr0.3CoO3−δ. J. Power Sources 2009,
191 (2), 275–279.
(59) Lankhorst, M. H. R.; Bouwmeester, H. J. M. Determination of Oxygen Nonstoichiometry
and Diffusivity in Mixed Conducting Oxides by Oxygen Coulometric Titration II. Oxygen
Nonstoichiometry and Defect Model for La0.8Sr0.2CoO3-Delta. J. Electrochem. Soc. 1997,
144 (4), 1268–1273.
(60) Ravindran, P.; Vidya, R.; FjellvĂĽg, H.; Kjekshus, A. Validity of Bond-Length and MĂśssbauer
Parameters to Assign Oxidation States in Multicomponent Oxides: Case Study of Sr4Fe4O11.
Phys. Rev. B 2008, 77 (13), 134448.
(61) Mizusaki, J.; Yoshihiro, M.; Yamauchi, S.; Fueki, K. Nonstoichiometry and Defect Structure
of the Perovskite-Type Oxides La1−xSrxFeO3−δ. J. Solid State Chem. 1985, 58 (2), 257–
266.
(62) Patrakeev, M. V.; Bahteeva, J. A.; Mitberg, E. B.; Leonidov, I. A.; Kozhevnikov, V. L.;
Poeppelmeier, K. R. Electron/Hole and Ion Transport in La1−xSrxFeO3−δ. J. Solid State
Chem. 2003, 172 (1), 219–231.

115

(63) Bahteeva, J. A.; Leonidov, I. A.; Patrakeev, M. V.; Mitberg, E. B.; Kozhevnikov, V. L.;
Poeppelmeier, K. R. High-Temperature Ion Transport in La1−xSrxFeO3-δ. J. Solid State
Electrochem. 2004, 8 (9), 578–584.
(64) Dann, S. E.; Currie, D. B.; Weller, M. T.; Thomas, M. F.; Al-Rawwas, A. D. The Effect of
Oxygen Stoichiometry on Phase Relations and Structure in the System La1-XSrxFeO3-δ (0
≤ x ≤ 1, 0 ≤ δ ≤ 0.5). J. Solid State Chem. 1994, 109 (1), 134–144.
(65) Fossdal, A.; Menon, M.; WĂŚrnhus, I.; Wiik, K.; Einarsrud, M.-A.; Grande, T. Crystal
Structure and Thermal Expansion of La1−xSrxFeO3−δ Materials. J. Am. Ceram. Soc. 2004,
87 (10), 1952–1958.
(66) Taylor, F. H.; Buckeridge, J.; Catlow, C. R. A. Defects and Oxide Ion Migration in the Solid
Oxide Fuel Cell Cathode Material LaFeO3. Chem. Mater. 2016, 28 (22), 8210–8220.
(67) Landau, L. D. Über die Bewegung der Elektronen in Kristalgitter. Phys Z Sowjetunion 1933,
3, 644–645.
(68) Jung, W.-H. Transport Mechanisms in La0.7Sr0.3FeO3: Evidence for Small Polaron
Formation. Phys. B Condens. Matter 2001, 299 (1), 120–123.
(69) Castleton, C. W. M.; Lee, A. L.; Kullgren, J.; Hermansson, K. Description of Polarons in
Ceria Using Density Functional Theory. J. Phys. Conf. Ser. 2014, 526 (1), 012002.
(70) Das, T.; Nicholas, J. D.; Qi, Y. In Preparation: Anisotropic Chemical Expansion in
Nonstoichiometric Ceria. 2017.
(71) Taylor, D. D.; Schreiber, N. J.; Levitas, B. D.; Xu, W.; Whitfield, P. S.; Rodriguez, E. E.
Oxygen Storage Properties of La1–xSrxFeO3−δ for Chemical-Looping Reactions—An In
Situ Neutron and Synchrotron X-Ray Study. Chem. Mater. 2016, 28 (11), 3951–3960.
(72) Ren, Y.; KĂźngas, R.; Gorte, R. J.; Deng, C. The Effect of A-Site Cation (Ln=La, Pr, Sm) on
the Crystal Structure, Conductivity and Oxygen Reduction Properties of Sr-Doped Ferrite
Perovskites. Solid State Ionics. 2012, 212, 47–54.
(73) Bishop, S. R. Chemical Expansion of Solid Oxide Fuel Cell Materials: A Brief Overview.
Acta Mech. Sin. 2013, 29 (3), 312–317.
(74) Bishop, S. R.; Stefanik, T. S.; Tuller, H. L. Electrical Conductivity and Defect Equilibria of
Pr0.1Ce0.9O2−δ. Phys. Chem. Chem. Phys. 2011, 13 (21), 10165.
(75) Korobko, R.; Patlolla, A.; Kossoy, A.; Wachtel, E.; Tuller, H. L.; Frenkel, A. I.; Lubomirsky,
I. Giant Electrostriction in Gd-Doped Ceria. Adv. Mater. 2012, 24 (43), 5857–5861.
(76) Nicholas, J. D.; Qi, Y.; Bishop, S. R.; Mukherjee, P. P. Introduction to Mechano-ElectroChemical Coupling in Energy Related Materials and Devices. J. Electrochem. Soc. 2014, 161
(11), Y11–Y12.

116

(77) Kresse, G.; FurthmĂźller, J. Efficient Iterative Schemes for Ab Initio Total-Energy
Calculations Using a Plane-Wave Basis Set. Phys. Rev. B 1996, 54 (16), 11169–11186.
(78) Kresse, G.; Joubert, D. From Ultrasoft Pseudopotentials to the Projector Augmented-Wave
Method. Phys. Rev. B 1999, 59 (3), 1758–1775.
(79) Blöchl, P. E. Projector Augmented-Wave Method. Phys. Rev. B 1994, 50 (24), 17953–17979.
(80) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple.
Phys. Rev. Lett. 1996, 77 (18), 3865–3868.
(81) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple
[Phys. Rev. Lett. 77, 3865 (1996)]. Phys. Rev. Lett. 1997, 78 (7), 1396–1396.
(82) Terakura, K.; Oguchi, T.; Williams, A. R.; Kübler, J. Band Theory of Insulating TransitionMetal Monoxides: Band-Structure Calculations. Phys. Rev. B 1984, 30 (8), 4734–4747.
(83) Anisimov, V. I.; Zaanen, J.; Andersen, O. K. Band Theory and Mott Insulators: Hubbard U
Instead of Stoner I. Phys. Rev. B 1991, 44 (3), 943–954.
(84) BlĂśchl, P. E.; Jepsen, O.; Andersen, O. K. Improved Tetrahedron Method for Brillouin-Zone
Integrations. Phys. Rev. B 1994, 49 (23), 16223–16233.
(85) Lee, C.-W.; Behera, R. K.; Wachsman, E. D.; Phillpot, S. R.; Sinnott, S. B. Stoichiometry of
the LaFeO3 (010) Surface Determined from First-Principles and Thermodynamic
Calculations. Phys. Rev. B 2011, 83 (11), 115418.
(86) Marezio, M.; Dernier, P. D. The Bond Lengths in LaFeO3. Mater. Res. Bull. 1971, 6 (1), 23–
29.
(87) Shimony, U.; Knudsen, J. M. MĂśssbauer Studies on Iron in the Perovskites La1-XSrxFeO3
(0<~x<~1). Phys. Rev. 1966, 144 (1), 361–366.
(88) Tang, W.; Sanville, E.; Henkelman, G. A Grid-Based Bader Analysis Algorithm without
Lattice Bias. J. Phys. Condens. Matter 2009, 21 (8), 084204.
(89) Bouibes, A.; Zaoui, A. High-Pressure Polymorphs of ZnCO3: Evolutionary Crystal Structure
Prediction. Sci. Rep. 2014, 4.
(90) Dudarev, S. L.; Botton, G. A.; Savrasov, S. Y.; Humphreys, C. J.; Sutton, A. P. ElectronEnergy-Loss Spectra and the Structural Stability of Nickel Oxide: An LSDA+U Study. Phys.
Rev. B 1998, 57 (3), 1505–1509.
(91) 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
(4), 045115.

117

(92) Wang, L.; Maxisch, T.; Ceder, G. Oxidation Energies of Transition Metal Oxides within the
GGA+U Framework. Phys. Rev. B 2006, 73 (19), 195107.
(93) Ceder, G.; Hautier, G.; Jain, A.; Ong, S. P. Recharging Lithium Battery Research with FirstPrinciples Methods. MRS Bull. 2011, 36 (3), 185–191.
(94) Anisimov, V. I.; Aryasetiawan, F.; Lichtenstein, A. I. First-Principles Calculations of the
Electronic Structure and Spectra of Strongly Correlated Systems: The LDA+ U Method. J.
Phys. Condens. Matter 1997, 9 (4), 767.
(95) 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 (1), 014103.
(96) Akriche, A.; Bouafia, H.; Hiadsi, S.; Abidri, B.; Sahli, B.; Elchikh, M.; Timaoui, M. A.;
Djebour, B. First-Principles Study of Mechanical, Exchange Interactions and the Robustness
in Co2MnSi Full Heusler Compounds. J. Magn. Magn. Mater. 2017, 422, 13–19.
(97) Bentouaf, A.; Mebsout, R.; Rached, H.; Amari, S.; Reshak, A. H.; AĂŻssa, B. Theoretical
Investigation of the Structural, Electronic, Magnetic and Elastic Properties of Binary Cubic
C15-Laves Phases TbX2 (X = Co and Fe). J. Alloys Compd. 2016, 689, 885–893.
(98) Gallagher, P. K.; MacChesney, J. B. Mössbauer Effect in the System Sr1–xLaxFeO3. Symp
Faraday Soc 1967, 1 (0), 40–47.
(99) Gryaznov, D.; Heifets, E.; Kotomin, E. The First-Principles Treatment of the ElectronCorrelation and Spin–orbital Effects in Uranium Mononitride Nuclear Fuels. Phys. Chem.
Chem. Phys. 2012, 14 (13), 4482–4490.
(100) Koehler, W. C.; Wollan, E. O. Neutron-Diffraction Study of the Magnetic Properties of
Perovskite-like Compounds LaBO3. J. Phys. Chem. Solids 1957, 2 (2), 100–106.
(101) Shein, I. R.; Shein, K. I.; Kozhevnikov, V. L.; Ivanovskii, A. L. Band Structure and the
Magnetic and Elastic Properties of SrFeO3 and LaFeO3 Perovskites. Phys. Solid State 2005,
47 (11), 2082–2088.
(102) Zhang, Q.; Yunoki, S. A First-Principles Study for Electronic and Magnetic Properties of
LaFeO3 /LaCrO3 Superlattices. J. Phys. Conf. Ser. 2012, 400 (3), 032126.
(103) Hoffmann, A.; Seo, J. W.; Fitzsimmons, M. R.; Siegwart, H.; Fompeyrine, J.; Locquet, J.P.; Dura, J. A.; Majkrzak, C. F. Induced Magnetic Moments at a FerromagnetAntiferromagnet Interface. Phys. Rev. B 2002, 66 (22), 220406.
(104) Takano, M.; Takeda, Y. Electronic State of Fe4+ Ions in Perovskite-Type Oxides. Bull. Inst.
Chem. Res. Kyoto Univ. 1983, 61 (5–6), 406–425.
(105) Reehuis, M.; Ulrich, C.; Maljuk, A.; Niedermayer, C.; Ouladdiaf, B.; Hoser, A.; Hofmann,
T.; Keimer, B. Neutron Diffraction Study of Spin and Charge Ordering in SrFeO3−δ. Phys.
Rev. B 2012, 85 (18), 184109.
118

(106) Xiao, G.; Liu, Q.; Wang, S.; Komvokis, V. G.; Amiridis, M. D.; Heyden, A.; Ma, S.; Chen,
F. Synthesis and Characterization of Mo-Doped SrFeO3−δ as Cathode Materials for Solid
Oxide Fuel Cells. J. Power Sources 2012, 202, 63–69.
(107) Takeda, T; Komura, S; Watanabe, N. In FERRITES: Proc. Int. Conf. Sept.-Oct. 1980;
Center for Academic Publication, Japan: Japan; p 385.
(108) Watanabe, H.; Oda, H.; Nakamura, E.; Yamaguchi, Y.; Takei, H. In FERRITES: Proc. Int.
Conf. Sept.-Oct. 1980; Center for Academic Publication, Japan: Japan; p 381.
(109) Takeda, T.; Watanabe, T.; Komura, S.; Fujii, H. Magnetic Properties of SrFe1-XCoxO3. J.
Phys. Soc. Jpn. 1987, 56 (2), 731–735.
(110) MuĂąoz-GarcĂ­a, A. B.; Pavone, M.; Carter, E. A. Effect of Antisite Defects on the Formation
of Oxygen Vacancies in Sr2FeMoO6 : Implications for Ion and Electron Transport. Chem.
Mater. 2011, 23 (20), 4525–4536.
(111) Yang, Z.; Huang, Z.; Ye, L.; Xie, X. Influence of Parameters U and J in the LSDA+U
Method on Electronic Structure of the Perovskites LaMO3 (M=Cr,Mn,Fe,Co,Ni). Phys. Rev.
B 1999, 60 (23), 15674–15682.
(112) Tsujimoto, Y.; Tassel, C.; Hayashi, N.; Watanabe, T.; Kageyama, H.; Yoshimura, K.;
Takano, M.; Ceretti, M.; Ritter, C.; Paulus, W. Infinite-Layer Iron Oxide with a SquarePlanar Coordination. Nature 2007, 450 (7172), 1062–1065.
(113) Wiβmann, S.; Becker, K. D. Localization of Electrons in Nonstoichiometric SrFeO3 − Δ.
Solid State Ionics. 1996, 85 (1–4), 279–283.
(114) David Sholl; Janice A Steckel. Density Functional Theory: A Practical Introduction, April
2009.; WILEY, 2009.
(115) Giannozzi, P.; Car, R.; Scoles, G. Oxygen Adsorption on Graphite and Nanotubes. J. Chem.
Phys. 2003, 118 (3), 1003–1006.
(116) Herzberg, G.; Huber, K. P. Molecular Spectra and Molecular Structure; Van Nostrand
Reinhold: New York ;London, 1979.
(117) Greg Mills, K. W. J. Nudged Elastic Band Method for Finding Minimum Energy Paths of
Transitions. 1998.
(118) Sheppard, D.; Terrell, R.; Henkelman, G. Optimization Methods for Finding Minimum
Energy Paths. J. Chem. Phys. 2008, 128 (13), 134106.
(119) Tucker, M. C. Progress in Metal-Supported Solid Oxide Fuel Cells: A Review. J. Power
Sources 2010, 195 (15), 4570–4582.

119

(120) Henkelman, G.; Uberuaga, B. P.; JĂłnsson, H. A Climbing Image Nudged Elastic Band
Method for Finding Saddle Points and Minimum Energy Paths. J. Chem. Phys. 2000, 113
(22), 9901–9904.
(121) Zhang, W.; Smith, J. R.; Wang, X.-G. Thermodynamics from Ab Initio Computations.
Phys. Rev. B 2004, 70 (2), 024103.
(122) 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
(2), 024112.
(123) Lee, Y.-L.; Kleis, J.; Rossmeisl, J.; Morgan, D. Ab Initio Energetics of LaBO3(001)
(B=Mn, Fe, Co, and Ni) for Solid Oxide Fuel Cell Cathodes. Phys. Rev. B 2009, 80 (22),
224101.
(124) Reuter, K.; Scheffler, M. Composition, Structure, and Stability of RuO2(110) as a Function
of Oxygen Pressure. Phys. Rev. B 2001, 65 (3), 035406.
(125) Lide, D. R. CRC Handbook of Chemistry and Physics; CRC Pr: Boca Raton, 1994.
(126) Chase, M. W. NIST-JANAF Themochemical Tables, Fourth Edition. J. Phys. Chem. Ref.
Data, Monograph 9 1998, 1–1951.
(127) JANAF Thermochemical Tables,2nd Edition, 2nd edition.; Stull, D. R., Prophet, H., Eds.;
National Bureau of Standards U.S, 1971.
(128) Xi, L.; Qiu, Y.; Zheng, S.; Shi, X.; Yang, J.; Chen, L.; Singh, D. J.; Yang, J.; Zhang, W.
Complex Doping of Group 13 Elements In and Ga in Caged Skutterudite CoSb3. Acta Mater.
2015, 85, 112–121.
(129) Loftager, S.; GarcĂ­a-Lastra, J. M.; Vegge, T. A Density Functional Theory Study of the
Ionic and Electronic Transport Mechanisms in LiFeBO3 Battery Electrodes. J. Phys. Chem.
C 2016, 120 (33), 18355–18364.
(130) Mantina, M.; Chen, L. Q.; Liu, Z. K. Predicting Diffusion Coefficients from First Principles
via Eyring’s Reaction Rate Theory. Defect Diffus. Forum 2009, 294, 1–13.
(131) Mantina, M.; Wang, Y.; Chen, L. Q.; Liu, Z. K.; Wolverton, C. First Principles Impurity
Diffusion Coefficients. Acta Mater. 2009, 57 (14), 4102–4108.
(132) Compaan, K.; Haven, Y. Correlation Factors for Diffusion in Solids. Trans. Faraday Soc.
1956, 52 (0), 786–801.
(133) Fabris, S.; de Gironcoli, S.; Baroni, S.; Vicario, G.; Balducci, G. Taming Multiple Valency
with Density Functionals: A Case Study of Defective Ceria. Phys. Rev. B 2005, 71 (4),
041102.

120

(134) Chen, L.-J.; Tang, Y.; Cui, L.; Ouyang, C.; Shi, S. Charge Transfer and Formation of Ce3+
upon Adsorption of Metal Atom M (M = Cu, Ag, Au) on CeO2 (100) Surface. J. Power
Sources 2013, 234, 69–81.
(135) Balaji Gopal, C.; GarcĂ­a-Melchor, M.; Lee, S. C.; Shi, Y.; Shavorskiy, A.; Monti, M.; Guan,
Z.; Sinclair, R.; Bluhm, H.; Vojvodic, A.; Chueh, W. C. Equilibrium Oxygen Storage
Capacity of Ultrathin CeO2-δ Depends Non-Monotonically on Large Biaxial Strain. Nat.
Commun. 2017, 8, 15360.
(136) Scavini, M.; Coduri, M.; Allieta, M.; Brunelli, M.; Ferrero, C. Probing Complex Disorder
in Ce 1- x Gd x O 2- x /2 Using the Pair Distribution Function Analysis. Chem. Mater. 2012,
24 (7), 1338–1345.
(137) Tang, Y.; Zhang, H.; Cui, L.; Ouyang, C.; Shi, S.; Tang, W.; Li, H.; Lee, J.-S.; Chen, L.
First-Principles Investigation on Redox Properties of M-Doped CeO2 ( M = Mn , Pr , Sn ,
Zr ). Phys. Rev. B 2010, 82 (12).
(138) Esch, F.; Fabris, S.; Zhou, L.; Montini, T.; Africh, C.; Fornasiero, P.; Comelli, G.; Rosei,
R. Electron Localization Determines Defect Formation on Ceria Substrates. Science 2005,
309 (5735), 752–755.
(139) Nolan, M.; Fearon, J.; Watson, G. Oxygen Vacancy Formation and Migration in Ceria.
Solid State Ionics. 2006, 177 (35–36), 3069–3074.
(140) Tang, Y.; Zhang, H.; Cui, L.; Ouyang, C.; Shi, S.; Tang, W.; Li, H.; Chen, L. Electronic
States of Metal (Cu, Ag, Au) Atom on CeO2(111) Surface: The Role of Local Structural
Distortion. J. Power Sources 2012, 197, 28–37.
(141) Castleton, C. W. M.; Kullgren, J.; Hermansson, K. Tuning LDA+U for Electron
Localization and Structure at Oxygen Vacancies in Ceria. J. Chem. Phys. 2007, 127 (24),
244704.
(142) Takeda, T.; Yamaguchi, Y.; Watanabe, H. Magnetic Structure of SrFeO3. J. Phys. Soc.
Jpn. 1972, 33 (4), 967–969.
(143) Yang, J. B.; Yelon, W. B.; James, W. J.; Chu, Z.; Kornecki, M.; Xie, Y. X.; Zhou, X. D.;
Anderson, H. U.; Joshi, A. G.; Malik, S. K. Crystal Structure, Magnetic Properties, and
Mössbauer Studies of La0.6Sr0.4FeO3−δ Prepared by Quenching in Different Atmospheres.
Phys. Rev. B 2002, 66 (18).
(144) Seo, J. W.; Fullerton, E. E.; Nolting, F.; Scholl, A.; Fompeyrine, J.; Locquet, J.-P.
Antiferromagnetic LaFeO3 Thin Films and Their Effect on Exchange Bias. J. Phys. Condens.
Matter 2008, 20 (26), 264014.
(145) Haas, O.; Vogt, U. F.; Soltmann, C.; Braun, A.; Yoon, W.-S.; Yang, X. Q.; Graule, T. The
Fe K-Edge X-Ray Absorption Characteristics of La1−xSrxFeO3−δ Prepared by Solid State
Reaction. Mater. Res. Bull. 2009, 44 (6), 1397–1404.

121

(146) Arima, T.; Tokura, Y.; Torrance, J. B. Variation of Optical Gaps in Perovskite-Type 3d
Transition-Metal Oxides. Phys. Rev. B 1993, 48 (23), 17006–17009.
(147) Galakhov, V. R.; Kurmaev, E. Z.; Kuepper, K.; Neumann, M.; McLeod, J. A.; Moewes,
A.; Leonidov, I. A.; Kozhevnikov, V. L. Valence Band Structure and X-Ray Spectra of
Oxygen-Deficient Ferrites SrFeOx. J. Phys. Chem. C 2010, 114 (11), 5154–5159.
(148) Waerenborgh, J. C.; Tsipis, E. V.; Auckett, J. E.; Ling, C. D.; Kharton, V. V. Magnetic
Structure of Sr2Fe2O5 Brownmillerite by Single-Crystal MĂśssbauer Spectroscopy. J. Solid
State Chem. 2013, 205, 5–9.
(149) Haavik, C.; Atake, T.; Stølen, S. On the Enthalpic Contribution to the Redox Energetics of
SrFeO3–δ. Phys. Chem. Chem. Phys. 2002, 4 (6), 1082–1087.
(150) Greaves, C.; Jacobson, A. J.; Tofield, B. C.; Fender, B. E. F. A Powder Neutron Diffraction
Investigation of the Nuclear and Magnetic Structure of Sr2Fe2O5. Acta Crystallogr. B 1975,
31 (3), 641–646.
(151) Nakamura, S.; Iida, S. New Data on Electrical Properties and Antiferromagnetism of
Highly Oxidized Perovskite SrFeOx ( 2.5<x<3.0). Jpn. J. Appl. Phys. 1995, 34 (Part 2, No.
3A), L291–L293.
(152) Okube, M.; Furukawa, Y.; Yoshiasa, A.; Hashimoto, T.; Sugahara, M.; Nakatsuka, A.
Oxidation State and Effective Pair Potential of Fe 4+ Ions in Perovskite-Type SrFeO3
Annealed under High Oxygen Pressure. J. Phys. Conf. Ser. 2008, 121 (9), 092004.
(153) Williams, G. V. M.; Hemery, E. K.; McCann, D. Magnetic and Transport Properties of
SrFeOx. Phys. Rev. B 2009, 79 (2), 024412.
(154) Wang, Q.-B.; Zhou, C.; Wu, J.; LĂź, T.; He, K.-H. GGA+U Study of the Electronic and
Optical Properties of Hexagonal BN Phase ZnO Under Pressure. Comput. Mater. Sci. 2015,
102, 196–201.
(155) Fuks, D.; Mastrikov, Y.; Kotomin, E.; Maier, J. Ab Initio Thermodynamic Study of
(Ba,Sr)(Co,Fe)O3 Perovskite Solid Solutions for Fuel Cell Applications. J. Mater. Chem. A
2013, 1 (45), 14320.
(156) Cheng, J.; Navrotsky, A.; Zhou, X.-D.; Anderson, H. U. Thermochemistry of La1XSrxFeO3-δ Solid Solutions (0.0 ≤ x ≤ 1.0, 0.0 ≤ δ ≤ 0.5). Chem. Mater. 2005, 17 (8), 2197–
2207.
(157) MacChesney, J. B.; Sherwood, R. C.; Potter, J. F. Electric and Magnetic Properties of the
Strontium Ferrates. J. Chem. Phys. 1965, 43 (6), 1907–1913.
(158) Mizusaki, J.; Okayasu, M.; Yamauchi, S.; Fueki, K. Nonstoichiometry and Phase
Relationship of the SrFeO2.5-SrFeO3 System at High Temperature. J. Solid State Chem.
1992, 99 (1), 166–172.

122

(159) Navrotsky, A. Energetics and Crystal Chemical Systematics among Ilmenite, Lithium
Niobate, and Perovskite Structures. Chem. Mater. 1998, 10 (10), 2787–2793.
(160) Das, T.; Nicholas, J. D.; Qi, Y. Polaron Size and Shape Effects on Oxygen Vacancy
Interactions in Lanthanum Strontium Ferrite. JMCA submitted 2017.
(161) Kawakami, T.; Nasu, S.; Sasaki, T.; Morimoto, S.; Endo, S.; Kawasaki, S.; Takano, M.
Charge Disproportionation and Magnetic Order of CaFeO3 under High Pressure up to 65
GPa. J. Phys. Soc. Jpn. 2001, 70 (6), 1491–1494.
(162) Mueller, D. N.; Machala, M. L.; Bluhm, H.; Chueh, W. C. Redox Activity of Surface
Oxygen Anions in Oxygen-Deficient Perovskite Oxides During Electrochemical Reactions.
Nat. Commun. 2015, 6, 6097.
(163) Hirai, K.; Kan, D.; Ichikawa, N.; Mibu, K.; Yoda, Y.; Andreeva, M.; Shimakawa, Y.
Strain-Induced Significant Increase in Metal-Insulator Transition Temperature in OxygenDeficient Fe Oxide Epitaxial Thin Films. Sci. Rep. 2015, 5, 7894.
(164) Mastrikov, Y. A.; Merkle, R.; Heifets, E.; Kotomin, E. A.; Maier, J. Pathways for Oxygen
Incorporation in Mixed Conducting Perovskites: A DFT-Based Mechanistic Analysis for (La,
Sr)MnO3−δ. J. Phys. Chem. C 2010, 114 (7), 3017–3027.
(165) Das, T.; Qi, Y.; Nicholas, J. D. In Preparation.
(166) ten Elshof, J. E.; Lankhorst, M. H. R.; Bouwmeester, H. J. M. Chemical Diffusion and
Oxygen Exchange of La0.6Sr0.4Co0.6Fe0.4O3−δ. Solid State Ionics. 1997, 99 (1–2), 15–22.
(167) Mosleh, M.; Søgaard, M.; Hendriksen, P. V. Kinetics and Mechanisms of Oxygen Surface
Exchange on La0.6Sr0.4FeO3 − δ Thin Films. J. Electrochem. Soc. 2009, 156 (4), B441–
B457.
(168) Søgaard, M.; Vang Hendriksen, P.; Mogensen, M. Oxygen Nonstoichiometry and
Transport Properties of Strontium Substituted Lanthanum Ferrite. J. Solid State Chem. 2007,
180 (4), 1489–1503.
(169) Søgaard, M.; Bieberle-Hßtter, A.; Hendriksen, P. V.; Mogensen, M.; Tuller, H. L. Oxygen
Incorporation in Porous Thin Films of Strontium Doped Lanthanum Ferrite. J.
Electroceramics 2011, 27 (3–4), 134–142.
(170) Chroneos, A.; Yildiz, B.; TarancĂłn, A.; Parfitt, D.; Kilner, J. A. Oxygen Diffusion in Solid
Oxide Fuel Cell Cathode and Electrolyte Materials: Mechanistic Insights from Atomistic
Simulations. Energy Environ. Sci. 2011, 4 (8), 2774.
(171) Wærnhus, I.; Grande, T.; Wiik, K. Surface Exchange of Oxygen in La1−x Sr x FeO3−δ (x
= 0, 0.1). Top. Catal. 2011, 54 (13–15), 1009–1015.

123

(172) Loschen, C.; Carrasco, J.; Neyman, K. M.; Illas, F. First-Principles LDA + U and GGA +
U Study of Cerium Oxides: Dependence on the Effective U Parameter. Phys. Rev. B 2007,
75 (3).
(173) Shi, S.; Tang, Y.; Ouyang, C.; Cui, L.; Xin, X.; Li, P.; Zhou, W.; Zhang, H.; Lei, M.; Chen,
L. O-Vacancy and Surface on CeO2: A First-Principles Study. J. Phys. Chem. Solids 2010,
71 (5), 788–796.
(174) Yavo, N.; Noiman, D.; Wachtel, E.; Kim, S.; Feldman, Y.; Lubomirsky, I.; Yeheskel, O.
Elastic Moduli of Pure and Gadolinium Doped Ceria Revisited: Sound Velocity
Measurements. Scr. Mater. 2016, 123, 86–89.
(175) Bishop, S. R.; Duncan, K. L.; Wachsman, E. D. Defect Equilibria and Chemical Expansion
in Non-Stoichiometric Undoped and Gadolinium-Doped Cerium Oxide. Electrochimica Acta
2009, 54 (5), 1436–1443.
(176) Khare, A.; Shin, D.; Yoo, T. S.; Kim, M.; Kang, T. D.; Lee, J.; Roh, S.; Jung, I.-H.; Hwang,
J.; Kim, S. W.; Noh, T. W.; Ohta, H.; Choi, W. S. Topotactic Metal-Insulator Transition in
Epitaxial SrFeOx Thin Films. Adv. Mater. 2017, 1606566.
(177) Amow, G.; Skinner, S. J. Recent Developments in Ruddlesden–Popper Nickelate Systems
for Solid Oxide Fuel Cell Cathodes. J. Solid State Electrochem. 2006, 10 (8), 538–546.
(178) Mori, K.; Kamiyama, T.; Kobayashi, H.; Torii, S.; Izumi, F.; Asano, H. Crystal Structure
of Sr3Fe2O7−δ. J. Phys. Chem. Solids 1999, 60 (8–9), 1443–1446.
(179) Adler, P. Electronic State, Magnetism, and Electrical Transport Behavior of
Sr3−xAxFe2O7(X≤0.4,A=Ba, La). J. Solid State Chem. 1997, 130 (1), 129–139.
(180) Maljuk, A.; Strempfer, J.; Ulrich, C.; Sofin, M.; Capogna, L.; Lin, C. T.; Keimer, B.
Growth of Sr3Fe2O7−x Single Crystals by the Floating Zone Method. J. Cryst. Growth 2004,
273 (1), 207–212.

124