disturbed
Disturbing Member
Full Member
- Messages
- 870
- Reaction score
- 56
This ones going to be difficult to explain. I have some more information to regurgitate(love that)..
I am now in contact with chemists at the local college and talking to Ivoclar to get the specifics on the molecular exchange between Zirliner and Zirconium, as I have proof that IF there is one it is VERY weak. more likely there is not.
FROM THE JOURNAL OF MATERIALS CHEMISTRY accepted march 2010 sooo... recent
this is complicated info..start with Section IV (fig 4) Migration barrier as a function of biaxial lattice strain. than read conclusion.
this is above a lot of peoples understanding so I will be producing studies from colleges across the globe and doing my best to explain the chemistry over the next few months...
and yes.. its another attempt to show you all that a bond other than mechanical retention is IMPOSSIBLE to zirconium and that their will be more chipping and failures to come. if you only care about your restorations lasting a max of 5 years than go ahead and ignore this.. it will probably only confuse and uspet you..
If I find that materials are added to zirconium oxide that allow for a chemical bond that gives an adequate bond strength I will post.. but as the vitrification of the material happens around 2700 degrees C I doubt it is possible.
View attachment KushimaYildiz_JMaterChem_2010.pdf
Oxygen ion diffusivity in strained yttria stabilized zirconia: where is the
fastest strain?
Akihiro Kushima and Bilge Yildiz*
Received 8th January 2010, Accepted 9th March 2010
First published as an Advance Article on the web 29th April 2010
DOI: 10.1039/c000259c
We present the mechanism and the extent of increase in the oxygen anion diffusivity in Y2O3 stabilized
ZrO2 (YSZ) under biaxial lattice strain. The oxygen vacancy migration paths and barriers in YSZ as
a function of lattice strain was assessed computationally using density functional theory (DFT) and
nudged elastic band (NEB) method. Two competing and non-linear processes acting in parallel were
identified to alter the migration barrier upon applied strain: (1) the change in the space, or electronic
density, along the migration path, and (2) the change in the strength of the interatomic bonds between the
migrating oxygen and the nearest neighbor cations that keep the oxygen from migrating. The increase of
the migration space and the weakening of the local oxygen–cation bonds correspond to a decrease of the
migration barrier, and vice versa. The contribution of the bond strength to the changes in the migration
barrier is more significant than that of the opening of migration space in strained YSZ. A database of
migration barrier energies as a function of lattice strain for a set of representative defect distributions in
the vicinity of the migration path inYSZwas constructed. This database was used in kinetic Monte Carlo
(KMC) simulations to estimate the effective oxygen diffusivity in strained YSZ. The oxygen diffusivity
exhibits an exponential increase up to a critical value of tensile strain, or the fastest strain. This increase is
more significant at the lower temperatures. At the strain states higher than the critical strain, the
diffusivity decreases. This is attributed to the local relaxations at large strain states beyond a limit of
elastic bond strain, resulting in the strengthening of the local oxygen–cation bonds that increases the
migration barrier. The highest enhancement of diffusivity in 9%-YSZ compared to its unstrained state is
6.8 103 times at 4% strain and at 400 K. The results indicate that inducing an optimal strain state by
direct mechanical load or by creating a coherent hetero-interface with lattice mismatch can enable
desirably high ionic conductivity in YSZ at reduced temperatures. The insights gained here particularly
on the nonlinear and competing consequences of lattice strain on the local bonding structure and charge
transport process are of importance for tuning the ionic transport properties in a variety of solid-state
conducting material applications, including but not limited to fuel cells.
Introduction
There is much interest in decreasing the Solid Oxide Fuel Cell
(SOFC) operating temperature to an intermediate-to-low
temperature range to improve chemical and mechanical stability
and reduce the cost. Enhancing the ionic transport in the electrolyte
and cathode materials is a key for this purpose. Over the
past several decades, developing new compositions with high
ionic mobility at low temperatures has been the primary activity.1
In recent years, however, nanoscale structures were shown as
means for potentially increasing the transport properties in
oxides of relevance to SOFCs.2–9 The importance of utilizing the
favorable interface properties of nanoscale structures shifted the
paradigm from developing new material compositions to
tailoring nano-structures using existing materials for the above
purpose.5,10 Among these, nanoscale hetero-layered oxides were
suggested to enhance the ionic conductivity along the interfaces
by orders of magnitude.5–9 The presence of an interface between
dissimilar oxides is commonly the source of charged mobile
defects and local anisotropic distortions and strained interatomic
bonds. These contribute to changes in the local electronic
structure; formation, coordination and redistribution of charged
defects; and mobility of defects in the vicinity of interfaces. Most
recently, an eight orders of magnitude increase in ionic conductance
was reported by Barriocanal et al.8 for 1–30 nm-thick yttria
stabilized zirconia (YSZ) layers coherently ‘‘strained’’ between
dielectric SrTiO3 (STO) layers. The exact nature of the ionic vs.
electronic conductance induced in the vicinity of the YSZ/STO
interface remains debatable.11 While these observations imply
that the hetero-interfaces with controlled lattice-strain and
defect-chemistry can play an important role in improving the
ionic mobility, results to date remain phenomenological without
systematic theoretical efforts. Beyond the SOFC materials
domain, similar questions for the role of hetero-interfaces on the
local bonding structure and the charge transport process are
open for tuning the transport properties in a variety of solid-state
conducting material applications.12–15 Systematic experimental
and theoretical investigation of the role of local strain states on
tuning the mechanism and magnitude of charge transfer properties
in these materials is prone to pursuit.
Department of Nuclear Science and Engineering, Massachusetts Institute
of Technology, 77 Massachusetts Avenue, Cambridge, MA, 02139, USA.
E-mail: [email protected]; Tel: +617-324-4009
This journal is ª The Royal Society of Chemistry 2010 J. Mater. Chem., 2010, 20, 4809–4819 | 4809
PAPER Journal of Materials Chemistry Home-High impact applications, properties and synthesis of exciting new materials | Journal of Materials Chemistry
Two main mechanisms contribute to increase the ionic transport
at the hetero-interface: (1) a favorable strain state at the
interface to shift and/or change the symmetry of electron energy
levels to provide for improved charge transfer16 and ion
mobility,17 and (2) the alteration of the defect chemistry7,18 near
the interface to enhance the density and distribution of desired
charge carriers. To understand the underlying mechanism of
oxygen anion transport in anion conducting ceramics, a significant
number of theoretical studies using atomistic simulations
based on both density functional theory (DFT) and empirical
potentials have been conducted (examples related to this work in
ref. 19–25). However, most of them are focused on the ‘‘bulk’’ or
the ‘‘surface’’ ionic transport. Only few theoretical studies exist
for the anion transport at hetero-interfaces, on specific cases,
i.e. CaF2/BaF2, CeO2/YSZ,20 with no generalization of results.
Of particular interest for this paper is the effect of lattice strain
on the ionic conductivity in the YSZ, as may be encountered due
to a lattice mismatch near the YSZ/STO interface. In probing the
effect of strain on conductivity in YSZ, a number of experimental
and simulation studies have been conducted in the past.21–24 Even
though the past results imply a potential increase in the
conductivity, the magnitude of increase was far lower than what
was reported recently in ref. 8 for the YSZ/STO interface.
Furthermore, the reported increases in magnitude of the
enhancement are inconsistent among these past studies because
of the differences in conditions, structures and methods that were
considered. For example, results were based on examining either
incoherent hetero-interfaces which cannot sustain large lattice
strain,21 or low strain conditions,24 or molecular dynamics
simulation with empirical potentials23,24 which are not appropriate
for interfaces, or phenomenological modeling using
macroscopic elastic properties of the material.22 In essence, the
mechanism by which a biaxial lattice strain alters the anion
transport in the YSZ fluorite structure was not uncovered. Given
these past inconsistencies and open questions, we adopt a firstprinciples-
based approach in probing the mechanistic picture of
the ionic conduction in strained YSZ. Our specific objectives are
to elucidate the mechanism and to identify the possible extent of
increase in the oxygen ion conductivity in YSZ due to a biaxial
lattice strain.
In order to attain reliable results in probing the effect of strain
on ionic transport in YSZ, it is important to consider a realistic
representation of the defect structure, in particular the distribution
of vacancies and dopant cations in YSZ. This is because the
oxygen vacancy migration barrier, which can be altered by strain,
depends strongly on the local defect distributions.25 On the other
hand, it is difficult to characterize experimentally the vacancy
and cation distribution in YSZ at atomic resolution, therefore
this type of information is not available from literature
straightforwardly. In spite of the difficulties, past work on this
topic provided useful insights on the vacancy–vacancy, vacancy–
cation and cation–cation interactions, enabling a reasonable
description of the defect structure in YSZ.
Here we provide a summary of defect–defect interactions that
we considered for constructing the YSZ model. The final YSZ
model used in our work should reasonably capture these reported
defect–defect interactions. First of these interactions
considered, namely the vacancy–cation interaction, determines
the relative position of the vacancy with respect to the Zr and Y
in YSZ. In earlier studies for this interaction, there were contradicting
experimental results; whether the vacancy is first
nearest neighbor26–28 (1NN) or second nearest neighbor29–31
(2NN) to the Y cations. However, later experimental30–33 and
theoretical34–40 studies indicated that the vacancy favors to bind
to the host Zr cations, resulting in a configuration such that the
vacancy is 1NN to the Zr and 2NN to the Y cations. Using DFT
calculations, Bogicevic and Wolverton37 showed that the local
relaxation of atoms was responsible for the vacancy to favor the
2NN position to the Y cation; opposite to what is expected from
the electrostatic interactions alone. The same trend was reported
in doped CeO2 when the dopant radius is larger than the host
cation,41 as in YSZ. For the vacancy–vacancy and the cation–
cation interactions, they concluded that the vacancies locate
themselves as far away as possible from each other, and that the
1NN Y–Y pair configurations are slightly more favorable than
the 2NN Y–Y pairs.36,38 More recent experimental work using
selected area diffraction also indicated a random distribution of
oxygen vacancies, maximizing the distance between them, in up
to 10% yttria doped YSZ.42 Consistent with the findings of
Bogicevic and Wolverton, theoretical results by Predith et al.43
suggest even the existence of the ordering of the dopant cations
(Y) in YSZ at high doping concentrations, and this cation
ordering could exist locally in YSZ at the lower Y concentration
range. Among these three types of defect–defect interactions, the
vacancy–vacancy interaction is the strongest followed by the
vacancy–cation interaction, and the cation–cation interaction is
the weakest of the three.36 Supporting the importance of the
vacancy–vacancy interactions, it was shown recently that the
vacancy distribution around the migration path could significantly
impact the value of the migration barriers in YSZ.44
In this study, we considered only the strain state in the bulk
to provide for improved ionic mobility in YSZ and did not
include an interface explicitly as in Barriocanal’s work.8
Therefore, the model in this paper excludes the effect of altered
defect chemistry near the interfaces of YSZ as a source of
increase in oxygen diffusivity. In the following sections, we
describe the simulation approach to calculate the oxygen
vacancy migration path and energy barriers. We report the
energetically favorable distributions of oxygen vacancies and Y
cations in the YSZ model structure to enable a realistic representation
of YSZ in this study. The mechanism by which the
lattice strain increases or decreases the anion transport in YSZ
are discussed. Because the vacancy migration barrier depends
on the local defect configurations, a database of migration
barrier energies as a function of lattice strain for a set of
representative defect distributions in the vicinity of the migration
path are constructed and presented. Finally, the effective
tracer diffusivity of oxygen in YSZ at different strain states is
presented and interpreted in comparison to the recent results
reported for the YSZ/STO interface.8
Simulation approach
The YSZ model was established to capture the charged defect
interactions and the defect distribution in the bulk YSZ, as
summarized in the Introduction. Upon identifying the stable
vacancy–cation distribution, the effect of the biaxial lattice
strain in the xy-plane on the migration barrier and on the
4810 | J. Mater. Chem., 2010, 20, 4809–4819 This journal is ª The Royal Society of Chemistry 2010
effective vacancy diffusivity were calculated and interpreted.
We used first-principles-based simulations to directly and
accurately probe the atomistic-scale nature of the vacancy
migration in the presence of lattice strain. The kinetic Monte
Carlo (KMC) method was used for calculating the effective
diffusivities.
The first principles calculations for identifying the energies
were performed by the Vienna ab initio simulation package
(VASP),45,46 employing density functional theory (DFT) using
a plane-wave basis set. Projector-augmented wave (PAW)
method47,48 with plane waves up to the energy cutoff at 400 eV
was used and the exchange–correlation energy was evaluated
by the generalized gradient approximation (GGA) using the
Perdew–Burke–Ernzerhof (PBE) function.49 Y, Zr and O atoms
were described by 11 (4s24p65s24d1),12 (4s24p65s24d2),and 6
(2s22p4) valence electrons, respectively. The G point in the Brillouin
zone was selected in the calculation. We confirmed the
energy convergence by comparing the energy obtained from each
of the k-point meshes of 3 3 3, 5 5 5 and 7 7 7. The
calculated energies are within 0.6% error compared to the energy
obtained using only the G point. The same condition was applied
to all of the simulations in this study unless stated otherwise.
YSZ has the fluorite crystal structure with the cations
(Y3+, Zr4+) occupying the FCC lattice sites and the O2 ions
occupying the tetragonal sites. For every two Y atoms, one
oxygen vacancy is created to satisfy charge neutrality. The
oxygen vacancy migrates by hopping to the adjacent tetragonal
site as schematically illustrated in Fig. 1. The migration barrier
depends on the neighboring atoms 1–6 in the figure.25 The
surrounding cations within a larger cutoff radius to the migration
path must also be accounted for in order to accurately quantify
the effect of the strain on the migration barrier.50 However, this
leads to an impractical number of patterns (well above 10 000
configurations just by including the next nearest neighbor cations
to the migration path) to be handled by the DFT simulation.
Therefore, we selected a subset of defect structure patterns near
the migration path, forming a database of migration paths and
barriers. This database serves as input for the KMC calculations
in estimating the effective diffusivity as a function of lattice strain
and temperature. In constructing this database, first, we
considered only two combinations of atoms at the sites 3 and 4,
as Zr–Zr and Zr–Y, in the migration path shown in Fig. 1. The
migration barrier through the Y–Y pair was found to be significantly
higher than these two combinations, and this makes the
Y–Y pair unlikely on a vacancy migration path. This selection
was also validated by the work of Car et al.,50 which implied the
Zr–Zr and Zr–Y subset at sites 3 and 4 to be sufficient to capture
the quantitative trends in ionic conductivity in YSZ as a function
of Y2O3 doping concentration. Second, we considered the defect
configurations outside of the first neighbor cations to the
migration path to affect the barrier. The particularly important
impact of the vacancy interactions on the migration barrier was
recently reported by Pietrucci et al.44 Therefore, in creating the
migration path and barrier database, we varied the relative
positions of the Y cations and the vacancies with respect to the
initial and final positions of the vacancy on the migration path.
The details of the migration path configurations in this database
are provided in the following section. For each migration path
configuration, nudged elastic band (NEB) method51 was used to
find the vacancy migration barrier as a function of biaxial strain
on the xy-plane. Given that the vacancy migration path in YSZ is
rather simple with single hops between the adjacent tetragonal
sites, the saddle point is expected near the midst of the migration
path.25,50 Three images between the initial and the final states
were taken in the simulations. Here we used the climbing image
NEB method52 to calculate not only the minimum energy path
but also the transition state configuration at the saddle point.
The method drives the image with the highest energy up to the
saddle point, by maximizing its energy along the band connected
to the image and minimizing in all other directions. This way, the
method allows us to accurately calculate the saddle point energy
and configuration with a small number of images needed in the
calculation. The NEB simulations were conducted with fixed
lattice vectors.
The database created for the migration barrier values at
different lattice strains and the defect configurations was
provided as input for kinetic Monte Carlo (KMC) simulations.
KMC calculations were performed to quantify the effective
diffusivity of oxygen in YSZ as a function of strain and
temperature. The diffusivity calculated in this work is the oxygen
tracer diffusion constant. We expect that the lattice strain has
a similar order of magnitude impact on the mobility of an oxygen
ion in both the chemical and the tracer diffusion constants,
especially in the range of low vacancy concentrations.53 In the
simulation model, oxygen atoms were assumed to hop to the
nearest neighbor site on the oxygen sublattice through the path
shown in Fig. 1. The rate of hops in the migration paths were
obtained by the Boltzman relationship as,
vAB ¼ n0 expEAB
kBT (1)
where EAB, found by the DFT-NEB simulations, denotes the
migration barrier across the A–B cation pair at sites 3 and 4 in
Fig. 1. The attempt frequency was fixed at n0 ¼ 1013 s1, which is
appropriate for most metal oxide systems.50 The periodic
boundary conditions were applied to all three dimensions of the
simulation cell. The vacancy diffusivity, Dv was extracted from
mean square displacement hR2i calculated in the KMC simulation
using the Einstein relation,
hR2i ¼ 6Dvt. (2)
The oxygen diffusivity DO was then obtained using Dv
considering the balance between the fluxes of oxygen atoms and
vacancies as,
Fig. 1 Migration path of the oxygen vacancy. Positions 1–6 are occupied
by the cations Zr or Y.
This journal is ª The Royal Society of Chemistry 2010 J. Mater. Chem., 2010, 20, 4809–4819 | 4811
DO ¼
cv
1 cv
Dv (3)
where cv is the vacancy concentration fraction. Dv exponentially
depends on the effective barrier for oxygen vacancy migration,
Eeff, as,
Dv ¼ dv
0 exp
Eeff
kBT (4)
where, d0v
is a constant. At small values of cv, D0 scales almost
linearly with cv. On the other hand, the change in the migration
barrier, Eeff, contributes exponentially to Dv and D0.
Results and discussion
This section consists of five sub-sections, I–V. Sub-sections I–II
demonstrate that the YSZ model in this work captures the
reported vacancy and cation interactions in stabilizing the cubic
structure in YSZ. The results in I and II are important for
obtaining reliable values of the energy barriers for vacancy
migration in YSZ. Sub-section III reports the effect of defect
configurations in the vicinity of the migrating vacancy on the
migration energy barrier. This proves that, in order to assess the
effect of strain on migration barriers, one must consider a realistic
subset of defect configurations in YSZ, and we report these
in sub-section IV. Sub-section IV in particular illustrates the
mechanism by which strain alters the diffusion barriers for
oxygen vacancies in YSZ. Finally sub-section V quantifies the
effective increase in oxygen tracer diffusivity as a function of
biaxial lattice strain, based on the subset of defect configurations
considered in this work for YSZ.
I. Defect structure in YSZ: oxygen vacancy and yttrium
distribution
To determine the favorable location of vacancies with respect to
the Y cations in the YSZ model, we tested a number of defect
configurations. These models constitute 30 Zr, 2 Y and 63 O
atoms, corresponding to a low doping ratio, 3% Y2O3, allowing
only one vacancy in the simulation cell, in order to eliminate the
vacancy–vacancy interaction in this part of the simulations. At
this low Y2O3 doping, the YSZ oxygen sublattice cannot keep its
cubic structure.34 Therefore, only in this part of the simulations,
the atoms beyond the second neighbor cation shell to Y cations
and beyond the third neighbor oxygen shell to the vacancy were
fixed to retain the cubic phase.34 The comparison of the energies
for the different YSZ models and defect configurations is presented
in Fig. 2. The energies are shifted such that the lowest
energy found is zero. The results show that the total energy can
vary by up to 0.5 eV when the Y cation is located at the third or
first neighbor site to the vacancy, instead of the second. This
result is consistent with the prior reports on cation–vacancy
interactions in YSZ, as discussed in the Introduction. Since this
large energy difference arises from local configuration differences,
the result suggests that the cation configuration
surrounding the diffusion path of the oxygen–vacancy pair
should be taken into account for accurate evaluation of the
migration barriers.
II. Distortion from the cubic lattice sites in YSZ
The atoms are not stabilized at the exact fluorite structure lattice
centers in YSZ, and large displacements of the cations and the
anions from the fluorite lattice sites are present.55–57 This is in
part due to the relaxation of the nearest neighbor oxygen anions
and cations towards and away from the vacancy, respectively.32
For an acceptable YSZ structure in this study, it is important to
capture this deviation quantitatively consistently with past
experimental and simulation reports. For this purpose, we
evaluated the dependence of the atomic displacements from the
ideal fluorite lattice sites on the local defect structures in a set of
candidate YSZ models. The fully relaxed structure of the models
constituted 9% Y2O3 doped YSZ, made of 26 Zr, 6 Y and 61 O
atoms in the simulation cell. The relaxation of the atoms were not
constrained in this and the following simulations because the
cubic structure is stable in the 9% Y2O3 doped YSZ, as long as
Fig. 2 Energy as a function of the neighboring location of the oxygen
vacancy with respect to the Y cations in YSZ. The energies were calculated
for different cation/vacancy distributions, and plotted such that the
lowest energy found is zero.
Fig. 3 Fully relaxed atomic configurations of the 9% yttria doped
zirconia with (a) nearest neighbor Y–Y pairs, and (b) uniformly
distributed Y. Yellow, gray and red spheres indicate Y, Zr and O ions.
Large distortion of the oxygen anion sublattice is evident when Y cations
are uniformly spread in the YSZ without any short-range ordering of the
Y cations.
4812 | J. Mater. Chem., 2010, 20, 4809–4819 This journal is ª The Royal Society of Chemistry 2010
the vacancies and dopant cations are distributed to enable the
stability. Fig. 3 shows the relaxed atomic configurations for two
models of YSZ. The first one (Fig. 3(a)) has first nearest neighbor
(1NN) Y–Y pairs and the second one (Fig. 3(b)) has a uniform
distribution of Y cations in the cell. In both models, the vacancies
were placed as far away as possible from each other and located
second nearest neighbor (2NN) to Y cations if allowed by the
available space. The first model retains an overall cubic structure
with oxygen anions displaced by 0.1–0.27 A from the ideal
fluorite lattice sites. The extent of the displacements calculated
here is well within the experimentally measured value of
0.36 A28 and prior simulation results of 0.20–0.31 A.34 On the
other hand, the second model shows large distortions, exemplified
at the circled region in Fig. 3(b),and has 1.2 eV higher energy
compared to the first one. This difference arises due to local
configuration differences, and thus, can affect the vacancy
migration barrier in YSZ which ranges from 0.2 to 1.4 eV.25
The lack of stability of the cubic oxygen sublattice in the second
model is due to the limited available space to accommodate the
relative positions of the vacancy and Y cations at the favorable
distances when the Y atoms are distributed uniformly into the
model. In the circled region in Fig. 3(b),a vacancy is located at
the 1NN position to the Y cations. This is because the 2NN
position for the vacancy is not available to all the Y atoms when
the Y atoms are distributed uniformly in the model. The
restriction of some vacancies to the unfavorable 1NN positions
with respect to the Y in YSZ destabilizes the cubic phase. On the
contrary, the first model configuration with the Y–Y pairs
introduces more space for the vacancies to be placed at the 2NN
positions with respect to the Y cations. This enables more
stability on the cubic oxygen sublattice. This result is consistent
with the locally ordered nature of Y cations in YSZ at high
concentrations of Y2O3, as found by Predith et al.,43 even if the
YSZ has only 9% Y2O3, here. While the effect of Y–Y pairing on
the migration barrier in a cubic YSZ is small compared to that of
the vacancy–vacancy and vacancy–cation interactions,36 it is
important to consider it to enable a stabilized cubic fluorite
structure of YSZ in the model.
III. Vacancy migration barriers in YSZ: effect of the local
defect structure
As implied in the previous section, the local defect structure
around the oxygen vacancy migration path in YSZ affects the
stability of the structure and the migration barrier,54 and this in
turn can determine the favorable paths for vacancy to migrate
through. To illustrate this, we present the differences in the
migration barrier in a model with 1NN Y–Y pairs, which lower
the energy and enable the stabilization of cubic phase. The
consequent migration paths considered were A / B / C in
Fig. 4(a). All six nearest neighbor cations (1–6 as shown in Fig. 1)
are the same, Zr, for each migration step in A/B/C. On the
other hand, the distribution of vacancies and Y cations near A
/B and B/C differ beyond the nearest neighbor cations. This
enables to investigate the variation of the migration barrier as
a function of the defect structure near the migration path beyond
the nearest neighbor cations. Fig. 4(b) shows the relative energy
in the NEB trajectories for the migration through A/B/C.
The resulting migration barriers range from 0.55 to 0.95 eV, even
though all the cations in the nearest neighbor position to the
migration path are the same. This significant difference exemplified
in the migration barriers in A/B/C indicates that it is
important to systematically take into account the longer-range
interactions between the migrating vacancy and its surrounding
cations and vacancies in YSZ.
We evaluated a subset of defect configurations near the
migration path to form a representative database of migration
barriers in YSZ. The five different configurations, models A–E,
that yield ten forward and backward jump barriers considered
here are shown in Fig. 5. This subset aims to capture the isolated
Fig. 4 (a) The vacancy positions and local structure considered for the
vacancy migration path, A / B / C. Large black spheres indicate
positions of the vacancies. Yellow, gray and red spheres indicate Y, Zr
and O ions. Numbers of Y cations at the second nearest neighbor position
to the vacancy is two at sites A and B, and one at site C. (b) NEB barriers
calculated along the A / B and B / C paths. The variation in the
migration barrier is clear in (b),even though all the six nearest neighbor
cations in the migration path are Zr. The difference in barriers is due to
the dissimilar vacancy and Y distribution near the A / B and B / C
paths beyond the nearest neighbor cations.
Fig. 5 Models A–E with varying defect configurations surrounding the
migration path. ‘‘# of 2NNY’’ denotes the number of second nearest
neighbor Y cations to the vacancy positions. ‘‘NN vacancy position’’
denotes that the nearest neighbor vacancy to the migrating vacancy
is located at n-th nearest neighbor position on the anion sublattice
(n-th ¼ 3rd or 4th in these models). Lower panel indicates the 2NN Y
position, and the 3rd and 4th nearest neighbor vacancies on the anion
sublattice with respect to the migrating vacancy V.
This journal is ª The Royal Society of Chemistry 2010 J. Mater. Chem., 2010, 20, 4809–4819 | 4813
and coupled effects of vacancy–vacancy and vacancy–cation
interactions in the local environment around the migration path.
The vacancy–vacancy interaction is considered up to the fourth
nearest neighbor on the anion sublattice, which was the
maximum that could be accommodated in our 2 2 2 unit cell
model when the nearest neighbor vacancy pairs are avoided. The
cation–vacancy interaction was taken up to the second nearest
neighbor. First, for considering the effect of the nearest neighbor
cations in the migration path (sites 1–6 in Fig. 1),we only
included the variation of the bridging cation pairs (sites 3 and 4
in Fig. 1). This was shown to be the most important contribution
to the vacancy migration in YSZ50 when the long-range interactions
were not considered. In models A, B, and E, the vacancy
migrates through the Zr–Zr cation pair, and in models C and D,
through the Zr–Y cation pair in the midst of the path. All the
other cation sites closest to the migration path (sites 1, 2, 5 and
6 in Fig. 1) are Zr. The path through Y–Y pair was excluded from
the simulations since it has a much higher migration barrier
compared to the Zr–Zr and Zr–Y pairs,25,50,54 making it significantly
less likely to migrate through the Y–Y cation pair. Second,
for vacancy–vacancy interaction, models A and C consider the
position of nearest neighbor vacancies, denoted as ‘‘NN vacancy
position’’ in Fig. 5, to the initial and the final position of the
migrating vacancy. The distance between the migrating vacancy
and the nearest neighbor vacancy changes from the fourth to the
third on the anion sublattice upon migrating from the initial to
the final position (and vice versa in the backward migration). # of
2NN Y are the same in models A and C. Third, for vacancy–
cation interaction beyond the nearest neighbor cations, models B
and D consider the number of Y cations that are at the second
nearest neighbor (2NN) position with respect to the initial and
the final state of the vacancy. This is denoted as ‘‘# of 2NN Y’’ in
Fig. 5. The number of 2NN Y cations changes from two to three
upon migrating from the initial to the final state (and vice versa in
the backward migration). Lastly, in model E, both the number of
2NN Y cations and the distance to the nearest vacancies at the
initial and the final state were varied. The former changes from
three to two and the latter from the fourth to the third nearest
neighbor (and vice versa in the backward migration).
In creating the database of migration paths and barriers,
a model to represent a symmetric defect configuration around the
initial and final states of the migrating vacancy was not possible.
This is because of the computational limitation in the simulation
cell size used in our DFT/NEB calculations. In computing the
effective diffusivity in the following sections, we will show that
we represented the barrier for a symmetric migration path to be
the average of the barriers from models A–E.
A simulation unit cell of 2 2 2 dimensions, containing
26 Zr, 6 Y and 61 O atoms corresponding to a 9% yttria doped
zirconia was used in the calculation of the migration energy
barriers for each of the models A–E described in Fig. 5. The
oxygen vacancies and the Y cations were distributed in the cell
satisfying the defect interactions discussed in the previous
sections. Atomic configurations and the cell vectors were fully
relaxed such that the total energy was minimized and the net
stress acting on the cell was zero. The migration barriers for each
model path were calculated using the climbing NEB technique.
The calculated barriers for each path (A–E) are shown in Fig. 6.
The energy values are relative to the initial state configuration.
The reaction coordinate of each NEB image was normalized with
respect to the distance between the initial and the final state of the
vacancy. The oxygen vacancy migration process requires
climbing over the highest energy (the saddle point) between the
initial and the final state. E!B and E B in Fig. 6 are the migration
barriers for the forward and the backward jumps, respectively,
corresponding to the migration from the initial to the final
positions of the vacancy shown in Fig. 5, and vice versa.
As a result, the migration configurations in models C and D
have higher barrier than those in models A, B and E. This is
primarily because of the difficulty in passing through the Zr–Y
pair compared to that through Zr–Zr pair. The forward jump
barrier is higher than the backward jump barrier for models
A and C, and lower for models B and D. For models A and C,
this is because the distance from the migrating vacancy to the
nearest neighbor vacancy is shorter at the final state, making it
more difficult for the migrating vacancy to approach to the final
state, primarily due the electrostatic interaction of vacancies. For
models B and D, the number of Y cations at the second nearest
neighbor position to the migrating vacancy is larger at the final
state, making this final state energetically more favorable, and
thus the reduced migration barrier in the forward migration.
In model E, the forward jump barrier is higher than the backward
one. Considering the vacancy–cation interactions, the final
state in model E is energetically more favorable than the initial
state because of the presence of a larger number of 2NN Y
cations to the migrating vacancy at the final state. This suggests
the migration barrier for the forward jump to be lower than the
backward one. However, the result shows the opposite behavior,
because the migrating vacancy is actually getting closer to its
nearest neighbor vacancy at its final state. This indicates that the
vacancy–vacancy interaction dominates the vacancy–cation
interaction in determining the energy barriers in model E. All
implications related to the defect–defect interactions leading to
the resulting migration barriers in models A–E in the unstrained
state of YSZ are consistent with those suggested by previous
Fig. 6 NEB paths and the corresponding energies in the migration path
calculated for models A–E. The relative energy and the coordinates are
normalized with respect to the initial energy, and the distance between the
initial and the final state, respectively. The saddle point energy is the
highest energy state between the initial and the final states, and E!B
and E B stand for barriers in the forward and the backward jumps,
respectively.
4814 | J. Mater. Chem., 2010, 20, 4809–4819 This journal is ª The Royal Society of Chemistry 2010
reports by Bogicevic et al.36 Therefore, the results presented here
enable a reasonable model of YSZ for assessing the effect of
biaxial lattice strain on the ionic mobility.
IV. Migration barrier as a function of biaxial lattice strain
In assessing the role of the biaxial xy-plane lattice strain on the
oxygen vacancy transport, we investigated the corresponding
changes in the migration barriers, and the underlying reasons for
these changes on models A–E in Fig. 5. The xy-plane strains
considered here were 0.00, 0.02, 0.04, 0.06 and 0.08. The range of
strain considered in this study is far above the experimentally
measured fracture strain of a single-crystal YSZ that is less than
0.5%.58 However, large strain in brittle ceramics could be achieved
by creating coherent hetero-interfaces at the nanoscale.
A maximum 7% tensile lattice strain was suggested to exist in the
YSZ thin film interfaces coherently strained between STO layers
due to the mismatch between STO and YSZ.8 On the other hand,
ability of thin-film (10–30 nm) YSZ to sustain this large tensile
strain is unclear. We considered only tensile strain to increase the
interatomic distances in YSZ, allowing more space for the
oxygen migration. Results for model A are shown in Fig. 7.
The migration barrier monotonically reduces with the increase in
the biaxial lattice strain. The barrier in the forward jump
decreases from 0.47 eV at 3 ¼ 0.00 to 0.07 eV at 3 ¼ 0.04, suggesting
a significant enhancement of the vacancy diffusivity
through this path under strain. Beyond this strain level, the
vacancy was not stable at the initial site and favored to form at
the defined final site in model A, with no barrier in the suggested
transition path. This result implies that the defect interactions
influencing the stability and energetics can also change with the
strain state of YSZ. While not explicitly included in our DFT
calculations, we implicitly took into account the effect of the
strain on the stable defect configurations in the KMC simulations,
as explained in the following section. We considered both
the forward and the backward migration barriers, and the system
has a higher probability of taking the lower energy migration
paths favoring the more stable structures.
The picture is not as straight-forward in models B–E as in
model A that shows a monotonic reduction in the migration
barrier with strain. Depending on the extent of the lattice strain,
the consequent barrier either decreases or increases in the
configurations represented by models B–E, as shown in Fig. 8.
The lack of a straight-forward trend in the magnitude of
migration barrier as a function of lattice strain suggests that there
is more than one mechanism acting simultaneously upon strain
in YSZ. Here, we explain the competing effects of strain on the
vacancy migration barrier in YSZ from an atomistic and electronic
perspective. The oxygen vacancy migration process
involves bond-breaking and bond-making between the migrating
oxygen and the surrounding cations. Therefore, we investigated
the changes in the electronic density distribution of the valence
electrons describing the Zr, Y and O atoms near the migration
path to explain the local effects of strain in YSZ. We focused on
the two planes in the migration path as shown in Fig. 9. The
charge distribution on the ‘‘oxygen plane’’ in Fig. 9 correlates to
the ‘‘space available’’ for the oxygen, denoted as O, to migrate
along the prescribed path. The migration space is also related to
the cation-cation distance midst the migration path. On the other
hand, the ‘‘bonding plane’’ in Fig. 9 shows a measure of the
‘‘bonding strength’’ between the oxygen and the nearest cations
at the initial site, quantified in terms of the oxygen–cation ‘‘bond
thickness’’. It is this set of bonds, shown between the oxygen
denoted as O and the cation denoted as C, which are broken in
the initiation of oxygen diffusion. Both the migration space and
the bond strength change when YSZ is subjected to biaxial strain.
Fig. 10 illustrates an example in which the increase in strain
increases the open space, and also increases the O–C bonding
strength to the extent of local inelastic relaxation near the initial
state of oxygen. In this case, the lattice strain is beyond just
Fig. 7 (a) NEB paths and the corresponding energies in the migration
path calculated for model A. The relative energy and the migration
coordinates are normalized with respect to the initial energy, and the
distance between the initial and the final state, respectively. The saddle
point energy is the highest energy state between the initial and the final
states. (b) The vacancy migration barrier, EB, as a function of biaxial
lattice strain. The filled circles and the open circles indicate forward and
backward jump, respectively. The barrier decrease as the strain increases
in the migration path in model A.
Fig. 8 The vacancy migration barrier, EB, as a function of biaxial lattice
strain for models B–E (in a–d). The filled circles and the open circles
indicate the forward and the backward jump barriers, respectively.
This journal is ª The Royal Society of Chemistry 2010 J. Mater. Chem., 2010, 20, 4809–4819 | 4815
elastically extending the oxygen–cation bond length; the oxygen,
O, was actually attracted towards the cation C and broke its
bond from the opposite cation upon local relaxation of atoms. As
a result, a stronger binding of O to C arises, shown with the
thicker O–C bond in Fig. 10. These two consequences of strain
act oppositely on the migration barrier – the former (opening of
migration space) favors the migration and the latter (local
inelastic relaxations near the oxygen) inhibits the migration.
A specific illustration of these competing consequences of
biaxial strain in YSZ represented by model E is shown in Fig. 11.
The resulting change in the migration space was quantified as the
distance between the cations midst the migration path, D(Zr–Zr)
in this configuration. The bond O–C in Fig. 11(b) is one of the
oxygen–zirconium bonds that break upon the migration of
oxygen. The large reduction in the forward migration barrier
from 0.51 eV at 3 ¼ 0.00 to 0.10 eV at 3 ¼ 0.06 as shown in
Fig. 8(d) corresponds to the increase in the D(Zr–Zr) distance
(shown in Fig. 11(a)) and the weakening of the O–C bond (shown
in Fig. 11(b)). On the other hand, this trend is not monotonic. At
higher strain, 3 ¼ 0.08, the cation–cation distance decreases and
the oxygen–zirconium bond strengthens again, resulting in
a significant increase of the migration barrier from 0.10 eV at
3 ¼ 0.06 to 0.47 eV at 3 ¼ 0.08. At this strain, the local relaxation
of the atoms particularly re-strengthens the oxygen-zirconium
bond at the initial site of the oxygen, also affecting the migration
space. If assumed independent, the oxygen–cation bond strength
is more strongly correlated to this large increase in the barrier at
large strain seen in Fig. 8(d). The increase in the cation–oxygen
bond strength is evident at 3 ¼ 0.08 in Fig. 11(b),but a very small
change occurs in the D(Zr–Zr) in Fig. 11(a). While we present this
set of results explicitly for the model E, the same qualitative
correlation of the competing effects of the migration space and
the cation–oxygen bond strength governed the migration barrier
changes with strain consistently also in the other models A–D.
V. Effective oxygen diffusivity as a function of biaxial lattice
strain
Models A–E constituted the database of the migration paths and
barriers as a function of the biaxial strain in YSZ. Using this
database to govern the overall oxygen diffusion, the macroscopic
oxygen diffusivity in the bulk of YSZ was estimated by KMC
simulations as described in the Simulation Approach section.
A 15 15 15 simulation cell of 8% yttria doped zirconia was
used, consistent with the doping ratio in YSZ reported in the
Ref. 8 for the YSZ/STO heterolayers. Y cations and vacancies
were distributed randomly in the KMC simulation cell. At a first
sight, this configuration contradicts with the implications for the
presence of 1NN Y–Y pairs in YSZ as discussed in Introduction.
On the other hand, the rules in selecting the diffusion barriers in
the KMC simulations were taken consistently with the local
defect configurations, including the effect of Y–Y pairs as evaluated
in the DFT-NEB calculations in this work. Therefore, we
considered this model in KMC to be reasonable for assessing the
overall oxygen diffusivity in YSZ. Supporting this approach,
past DFT and KMC studies of oxygen diffusivity in YSZ with
randomly distributed Y atoms and vacancies were shown to give
Fig. 9 The electronic density distribution on the ‘‘oxygen plane’’ and the
‘‘bonding plane’’ at the initial and the saddle states. The migrating
oxygen, O, moves towards the open migration space, and the bond
between oxygen O and cation C breaks in the migration.
Fig. 10 Change in the electronic density distribution on the oxygen
plane and the bonding plane under biaxial lattice strain, 3, from 0.00 to
0.07. Gray and red spheres represent Zr and O atoms, respectively.
4816 | J. Mater. Chem., 2010, 20, 4809–4819 This journal is ª The Royal Society of Chemistry 2010
the correct yttria doping ratio for the maximum oxygen
diffusivity.25,50
The effect of the local defect environment around the migration
path on the barriers input to the KMC simulations was
taken into account as follows:
Definitions. E!a and E a are the forward and backward jump
energy barriers in each model a (a ¼ A, B, C, D, E),respectively.
ri
v and rf
v stand for the distance of the first nearest neighbor (1NN)
vacancy to the migrating vacancy at the initial and final state,
respectively. ni
Y, and nf
Y are the number of second nearest
neighbor (2NN) Y cations to the migrating vacancy at the initial
and final state, respectively.
Rules. The rules in determining the energy barrier for each hop
in the KMC simulations treat the local defect configurations in
a coarse grained fashion, based on the models A–E (Fig. 5)
assumed to govern the oxygen vacancy diffusion in YSZ.
In models A and C, the vacancy forward hopping (from initial to
final state) results in a decrease of the 1NN vacancy distance to
the migrating vacancy, that is ri
v > rf
v. Then, in the KMC simulation,
E!a and E a (a ¼ A, C) are chosen for the cases of ri
v > rf
v,
and ri
v<rf
v, respectively. In models B and D, the vacancy forward
hopping (from initial to final state) results in an increase in the
number of 2NN Y to the migrating vacancy, that is ni
Y < nf
Y.
Then, in the KMC simulation, E!a and E a (a ¼ B, D) are chosen
for the cases of ni
Y <nf
Y and ni
Y > nf
Y, respectively. The explicit
definition of the rules implemented in the KMC simulations is as
follows:
(1) For the vacancy migration through a Zr–Zr pair midst the
migration path, we select the barrier among E!A, E A E!B and E B
according to the following criteria:
If ri
v > rf
v, select E!A; and if ri
v < rf
v, select E A.
If ri
v ¼ rf
v and ni
Y > nf
Y, select E!B; or if ri
v ¼ rf
v and ni
Y < nf
Y,
select E B.
If ri
v ¼ rf
v and ni
Y ¼ nf
Y, randomly choose barriers from E!A,
E A, E!B, and E B.
(2) For the vacancy migration through a Zr–Y pair midst the
migration path, we select the barrier among E!C, E C, E!D and E D
according to the following criteria:
If ri
v > rf
v, select E!C; if ri
v < rf
v, select E C.
If ri
v ¼ rf
v and ni
Y > nf
Y, select E!D; or if ri
v ¼ rf
v and ni
Y < nf
Y,
select E D.
If ri
v ¼ rf
v and ni
Y ¼ nf
Y, randomly choose barriers from E!C,
E C, E!D and E D.
(3) Vacancy migration through Y–Y pair is prohibited.
When ri
v ¼ rf
v and ni
Y ¼ nf
Y, the migration path is surrounded by
a symmetric distribution of defects. This case was not possible to
account for in the DFT-NEB calculations due to the computational
limitations on the simulation size. We assumed such
a symmetric path to be governed by a migration barrier that is
the average of the barriers from the non-symmetric models represented
in A–E here. Therefore, the random selections among
the entire set of E!a and E a (for a ¼ A, B, C, D, E) during the 108
KMC steps should provide sufficient statistics to represent the
average migration barrier for the symmetric path, ri
v ¼ rf
v and
ni
Y ¼ nf
Y. Moreover, the probability of having symmetric path in
the KMC runs was less than 3%. Therefore, the error that may
arise due to selecting an average migration barrier for the
symmetric configuration is expected to have little significance on
the effective diffusivity.
As seen from the rules above, the energy barrier for each hop
in the KMC simulations was chosen based on a coarse grained
treatment of the local defect configurations around the migrating
vacancy. We believe this approach gives a reasonable result,
particularly based on the validity of the models A–E (Fig. 5) to
govern the oxygen vacancy diffusion in YSZ. On the other hand,
the reader can refer to a more accurate treatment of the diffusion
barriers in KMC simulations using the cluster expansion
method,43,59 which can account for the energy deviations due to
the local atomic structures around each hop. In that case, the
KMC model itself has to take into account the stable distributions
of vacancies and cations (as in subsections I and II),rather
than a random distribution, making the task exhaustive.
The KMC simulations were performed as a function of
temperature ranging as 400, 600, 800 and 1000 K, and of biaxial
strain ranging as 0.00, 0.02, 0.04, 0.06 and 0.08. Ten simulations
with different distributions of dopant cations and vacancies were
conducted for each temperature and strain state. The resulting
strain dependence of the effective oxygen tracer diffusivity DO is
shown in Fig. 12. y-axis in the figure is on a relative scale with
Fig. 11 (a) Change in the cation–cation distance, D(Zr–Zr),and (b) the
bond thickness (strength) for O–C in the migration under biaxial lattice
strain from 0.00 to 0.08 for the forward jump in model E (see also
Fig. 8(d)). Increase in the migration space D(Zr–Zr),and the weakening
of the O–C bond correspond to the decrease of the migration barrier in
Fig. 8(d),and vice versa. Large gray spheres and small red spheres are Zr
and O atoms, respectively. The isoline interval is 0.5 A3.
This journal is ª The Royal Society of Chemistry 2010 J. Mater. Chem., 2010, 20, 4809–4819 | 4817
respect to the oxygen diffusivity, D0
O at 3 ¼ 0.00 for the lowest
temperature simulated, 400 K. The reported results show the
average values and the error bars over the ten simulations. The
result is not a simple monotonic increase of diffusivity with
increasing strain in the entire range considered. The oxygen
diffusivity exhibits an exponential increase up to a critical value,
4% of tensile strain. This increase is more significant at the lower
temperatures. The maximum enhancement in diffusivity
compared to the unstrained state in YSZ is 6.8 103 times higher
at 4% strain and at 400 K. At the higher strain states, the
diffusivity decreases. This is attributed to the local relaxations
beyond a limit of elastic bond strain, resulting in the strengthening
of the local oxygen–cation bonds that increases the
migration barrier.
Our first-principles-based simulation results indicate that
a biaxial strain in YSZ induced by a lattice-mismatch coherent
interface can significantly increase the diffusivity up to a critical
strain state, and worsen it beyond that strain state. One could
view this critical strain state as the optimal strain to induce in the
material to attain the maximum benefit of the biaxial lattice
strain on the oxygen diffusivity. We note that a recent prediction
by Schichtel et al.22 using a model based on the elastic properties
of YSZ estimated 2.5 orders of magnitude increase in the ionic
conductivity of YSZ at 7% strain at 573 K. The model of
Schichtel et al. relies on the isotropic pressure induced by the
strain to enhance the diffusivity, and does not take into account
the competing relaxations at these large tensile strain states in
YSZ. Therefore, it cannot capture the local bond-strength
changes that may mitigate the diffusivity at large strain values.
In spite of the differences in our and Schichtel’s model, the
magnitude of the maximum increase in the oxygen diffusivity in
strained YSZ predicted by these models range from 2.5 to 3.8
orders of magnitude. Therefore, the lattice strain (elastic and
inelastic) alone cannot explain the eight orders of magnitude
increase of ionic conductivity in YSZ reported by Barriocanal
et al.8 These results do not take into account the space-charge
effects in altering the anion concentrations near the YSZ interface.
On the other hand, the contribution to the increase in the
anion diffusivity due to a possible space-charge-induced increase
in the vacancy concentration is almost linear at low concentration
of vacancies (see eqn (3)) as in the 9% Y2O3 doped YSZ.
Therefore, this contribution is expected to be small compared to
the exponential increase found here due to the applied biaxial
lattice strain.60
Conclusion
This work provides a mechanistic and systematic assessment of
the effect of biaxial lattice strain on the oxygen diffusivity in
YSZ, independently from the defect chemistry near a heterointerface
with a dissimilar oxide, e.g. YSZ/STO as reported by
Barriocanal et al.8 Using DFT/NEB simulations, we uncovered
the underlying microscopic mechanism by which the biaxial
tensile strain acts on the oxygen vacancy transport in YSZ. The
migration barrier correlates fundamentally with two competing
processes acting in parallel in the presence of lattice strain: (1) the
migration space, quantified as the cation–cation distance midst
the migration path for the oxygen, and (2) the strength of the
bond between the oxygen and the nearest cation preventing the
oxygen from migrating. Increasing of the migration space and
weakening of the oxygen–cation bond contribute to the decrease
of the vacancy migration barrier, and vice versa. In assessing the
effective diffusivity as a function of strain in YSZ, we constructed
a data base of migration paths and barriers to account for the
important effects of local distributions of vacancies and dopant
cations near the migration path. KMC simulations using this
database showed that the macroscopic oxygen diffusivity
increases exponentially in YSZ up to a critical value of biaxial
tensile strain, with a more significant impact at the lower
temperatures. Beyond this critical value, the strain state in YSZ
worsens the diffusivity due to the local relaxations that
strengthen the oxygen–cation bond at the initiation of migration.
One could view this critical strain state as the optimal strain, or
the fastest strain, to attain the highest benefit of the biaxial lattice
strain on the oxygen diffusivity in the material system. The
maximum enhancement in diffusivity in YSZ was 6.8 103 times
higher at 4% strain and temperature of 400 K. This is a significant
and desirable increase. Nevertheless, it is far below the eight
orders of magnitude increase reported by Barriocanal et al. for
the thin-film YSZ layers coherently strained between STO
layers.8 Therefore, we believe that the lattice strain alone cannot
be responsible for the colossally high ionic conductivity suggested
to occur in YSZ in ref. 8. Nevertheless, given our results,
controlling the lattice strain by direct mechanical load or by
creating a coherent hetero-interface to induce an optimal strain
state can enable desirably high ionic conductivity at reduced
temperatures in the fluorite material, the YSZ in this case.
The insights gained here for the role of lattice strain on the
local bonding structure and ionic charge transport process are of
importance for modulating the transport properties in a variety
of solid state conducting material applications, beyond the
specific SOFC electrolyte material studied in this work. Interesting
examples may extend to Li-ion batteries, oxygen gas
sensors, semi-conductors and photovoltaics. We believe that the
mechanistic results obtained in this study, particularly on the
Fig. 12 The relative enhancement of the oxygen diffusivity calculated by
the kinetic Monte Carlo simulations at different temperatures and strain
states. D0
O is the oxygen diffusivity at 3 ¼ 0.00 at 400 K. The diffusivity
increases exponentially with strain, and shows a maximum at (or near)
3 ¼ 0.04 followed by a decrease at higher strains. The increase in diffusivity
is more significant at lower temperatures.
4818 | J. Mater. Chem., 2010, 20, 4809–4819 This journal is ª The Royal Society of Chemistry 2010
non-linear and competing consequences of strain on the bonding
structure in the vicinity of the ionic carrier species, provides
insights for creating functional strained nano-structures for ionic
and electronic conducting devices, including but not limited to
SOFCs.
Acknowledgements
We thank the Nuclear Regulatory Commission for financial
support, and the National Science Foundation for computational
support through the TeraGrid Advanced Support
Program under Grant No. TG-ASC090058.
References
1 J. C. Boivin and G. Mairesse, Chem. Mater., 1998, 10, 2870.
2 T. Suzuki, I. Kosacki and H. U. Anderson, Solid State Ionics, 2002,
151, 111.
3 I. Kosacki, C. M. Rouleau, P. J. Becher, J. Bentley and
D. H. Lowndes, Solid State Ionics, 2005, 176, 1319.
4 J. H. Shim, C.-C. Chao, H. Huang and F. B. Printz, Chem. Mater.,
2007, 19, 3850.
5 J. Maier, Nat. Mater., 2005, 4, 805.
6 N. Sata, K. Eberman, K. Eberl and J. Maier, Nature, 2000, 408, 946.
7 X. X. Guo, I. Matei, J.-S. Lee and J. Maier, Appl. Phys. Lett., 2007,
91, 103102.
8 J. G. Barriocanal, A. R. Calzada, M. Varela, Z. Sefrioui, E. Iborra,
C. Leon, S. J. Pennycook and J. Santamaria, Science, 2008, 321, 676.
9 S. Ramanathan, J. Vac. Sci. Technol., A, 2009, 27, 1126.
10 A. J. Jacobson, Chem. Mater., 2010, 22, 660.
11 X. Guo, Science, 2009, 324, 465; J. G. Barriocanal, A. R. Calzada,
M. Varela, Z. Sefrioui, E. Iborra, C. Leon, S. J. Pennycook and
J. Santamaria, Science, 2009, 324, 465.
12 M. P. O’Callaghan, A. S. Powell, J. J. Titman, G. Z. Chen and
E. J. Cussen, Chem. Mater., 2008, 20, 2360.
13 Y. Hoshina, K. Iwasaki, A. Yamada and M. Konagai, Jpn. J. Appl.
Phys., 2009, 48, 04C125.
14 J. Munguıa, G. Bremond, J. M. Bluet, J. M. Hartmann and
M. Mermoux, Appl. Phys. Lett., 2008, 93, 102101.
15 V. R. D’Costa, Y.-Y. Fang, J. Tolle, J. Kouvetakis and J. Menendez,
Phys. Rev. Lett., 2009, 102, 107403.
16 M. Mavrikakis, P. Stoltze and J. K. Nørskov, Catal. Lett., 2000, 64,
101.
17 J. Kilner, Nat. Mater., 2008, 7, 838.
18 J. Maier, Prog. Solid State Chem., 1995, 23, 171.
19 A. K. Ivanov-****z, Crystallogr. Rep., 2007, 52, 129.
20 D. C. Sayle, J. A. Doig, S. C. Parker and G. W. Watson, Chem.
Commun., 2003, 1804.
21 C. Korte, A. Peters, J. Janek, D. Hesse and N. Zakharov, Phys. Chem.
Chem. Phys., 2008, 10, 4623.
22 N. Schichtel, C. Korte, D. Hesse and J. Janek, Phys. Chem. Chem.
Phys., 2009, 11, 3043.
23 K. Suzuki, M. Kubo, Y. Oumi, R. Miura, H. Takaba, A. Fahmi,
A. Chatterjee, K. Teraishi and A. Miyamoto, Appl. Phys. Lett.,
1998, 73, 1502.
24 W. Araki, Y. Imai and T. Adachi, J. Eur. Ceram. Soc., 2009, 29, 2275.
25 R. Pornprasertsuk, P. Ramanarayanan, C. B. Musgrave and
F. B. J. Prinz, J. Appl. Phys., 2005, 98, 103513.
26 W. L. Roth, R. Wong, A. I. Goldman, E. Canova, Y. H. Kao and
B. Dunn, Solid State Ionics, 1986, 18–19, 1115; H. Morikawa,
Y. Shimizugawa, F. Marumo, T. Harasawa, H. Ikawa, K. Tohji
and Y. Udagawa, J. Ceram. Soc. Jpn, 1988, 96, 253; M. H. Tuilier,
J. Dexpert-Ghys, H. Dexpert and P. Lagarde, J. Solid State Chem.,
1987, 69, 153.
27 M. Weller, Z. Metallkd., 1993, 84, 6.
28 D. Steele and B. E. F. Fender, J. Phys. C: Solid State Phys., 1974, 7, 1.
29 B. W. Veal, A. G. McKale, A. P. Paulikas, S. J. Rothman and
L. J. Nowicki, Physica B+C, 1988, 150, 234; D. Komyoji,
A. Yoshiasa, T. Moriga, S. Emura, F. Kanamaru and K. Koto,
Solid State Ionics, 1992, 50, 291.
30 C. R. A. Catlow, A. V. Chadwick, G. N. Greaves and L. M. Moroney,
J. Am. Ceram. Soc., 1986, 69, 272.
31 P. Li and J. E. Penner-Hahn, Phys. Rev. B: Condens. Matter, 1993, 48,
10063; P. Li and J. E. Penner-Hahn, Phys. Rev. B: Condens. Matter,
1993, 48, 10074; P. Li and J. E. Penner-Hahn, Phys. Rev. B:
Condens. Matter, 1993, 48, 10082.
32 J. P. Goff, W. Hayes, S. Hull, M. T. Hutchings and K. N. Clausen,
Phys. Rev. B, 1999, 59, 14202.
33 K. Kawata, H. Maekawa, T. Nemoto and T. Yamamura, Solid State
Ionics, 2006, 177, 1687.
34 G. Stapper, M. Bernasconi, N. Nicoloso and M. Parrinello, Phys.
Rev. B: Condens. Matter Mater. Phys., 1999, 59, 797.
35 A. Bogicevic, C. Wolverton, G. M. Crosbie and E. B. Stechel, Phys.
Rev. B: Condens. Matter Mater. Phys., 2001, 64, 014106.
36 A. Bogicevic and C. Wolverton, Phys. Rev. B: Condens. Matter
Mater. Phys., 2003, 67, 024106.
37 A. Bogicevic and C. Wolverton, Europhys. Lett., 2001, 56, 393.
38 X. Xia, R. Oldman and R. Catlow, Chem. Mater., 2009, 21, 3576.
39 M. Sakib Khan, M. Saiful Islam and D. R. Bates, J. Mater. Chem.,
1998, 8, 2299.
40 M. O. Zacate, L. Minervini, D. J. Bradfield, R. W. Grimes and
K. E. Sickafus, Solid State Ionics, 2000, 128, 243.
41 D. A. Andersson, S. I. Simak, N. V. Skorodumova, I. A. Abrikosov
and B. Johansson, Proc. Natl. Acad. Sci. U. S. A., 2006, 103, 3518.
42 S. Garcıa-Martın, D. P. Fagg and J. T. S. Irvine, Chem. Mater., 2008,
20, 5933.
43 A. Predith, G. Ceder, C. Wolverton, K. Persson and T. Mueller, Phys.
Rev. B: Condens. Matter Mater. Phys., 2008, 77, 144104.
44 F. Pietrucci, M. Bernasconi, A. Laio and M. Parrinello, Phys. Rev. B:
Condens. Matter Mater. Phys., 2008, 78, 094301.
45 G. Kresse and J. Hafner, Phys. Rev.B: Condens. Matter, 1993, 47, 558.
46 G. Kresse and J. Furthm€uller, Phys. Rev. B: Condens. Matter, 1996,
54, 11169.
47 P. E. Bl€ochl, Phys. Rev. B: Condens. Matter, 1994, 50, 17953.
48 G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater.
Phys., 1999, 59, 1758.
49 J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77,
3865.
50 R. Krishnamurthy, Y.-G. Yoon, D. J. Srolovitz and R. Car, J. Am.
Ceram. Soc., 2004, 87, 1821.
51 G. Mills and H. Jonsson, Phys. Rev. Lett., 1994, 72, 1124.
52 G. Henkelman, B. P. Uberuaga and H. Jonsson, J. Chem. Phys., 2000,
113, 9901.
53 A. Van der Ven, G. Ceder, M. Asta and P. D. Tepesch, Phys. Rev. B:
Condens. Matter Mater. Phys., 2001, 64, 184307.
54 M. Martin, J. Electroceram., 2006, 17, 765.
55 F. Frey, H. Boysen and I. K- Bischoff, Z. Kristallogr., 2005, 220, 1017.
56 D. N. Argyriou, M. M. Elcombe and A. C. Larson, J. Phys. Chem.
Solids, 1996, 57, 183.
57 N. Ishizawa, Y. Matsushima, M. Hayashi and M. Ueki, Acta
Crystallogr., Sect. B: Struct. Sci., 1999, 55, 726.
58 D. Baither, M. Bartsch, B. Baufeld, A. Tikhonovsky, M. R€uhle and
U. Messerschmidta, J. Am. Ceram. Soc., 2001, 84, 1755.
59 A. Van der Ven, J. C. Thomas, Q. C. Xu, B. Swoboda and
D. Morgan, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 78,
104306.
60 A. Kushima and B. Yildiz, ECS Trans., 2009, 25, 1599.
This journal is ª The Royal Society of Chemistry 2010 J. Mater. Chem., 2010, 20, 4809–4819 | 4819
heres the first one..more to come..
I am now in contact with chemists at the local college and talking to Ivoclar to get the specifics on the molecular exchange between Zirliner and Zirconium, as I have proof that IF there is one it is VERY weak. more likely there is not.
FROM THE JOURNAL OF MATERIALS CHEMISTRY accepted march 2010 sooo... recent
this is complicated info..start with Section IV (fig 4) Migration barrier as a function of biaxial lattice strain. than read conclusion.
this is above a lot of peoples understanding so I will be producing studies from colleges across the globe and doing my best to explain the chemistry over the next few months...
and yes.. its another attempt to show you all that a bond other than mechanical retention is IMPOSSIBLE to zirconium and that their will be more chipping and failures to come. if you only care about your restorations lasting a max of 5 years than go ahead and ignore this.. it will probably only confuse and uspet you..
If I find that materials are added to zirconium oxide that allow for a chemical bond that gives an adequate bond strength I will post.. but as the vitrification of the material happens around 2700 degrees C I doubt it is possible.
View attachment KushimaYildiz_JMaterChem_2010.pdf
Oxygen ion diffusivity in strained yttria stabilized zirconia: where is the
fastest strain?
Akihiro Kushima and Bilge Yildiz*
Received 8th January 2010, Accepted 9th March 2010
First published as an Advance Article on the web 29th April 2010
DOI: 10.1039/c000259c
We present the mechanism and the extent of increase in the oxygen anion diffusivity in Y2O3 stabilized
ZrO2 (YSZ) under biaxial lattice strain. The oxygen vacancy migration paths and barriers in YSZ as
a function of lattice strain was assessed computationally using density functional theory (DFT) and
nudged elastic band (NEB) method. Two competing and non-linear processes acting in parallel were
identified to alter the migration barrier upon applied strain: (1) the change in the space, or electronic
density, along the migration path, and (2) the change in the strength of the interatomic bonds between the
migrating oxygen and the nearest neighbor cations that keep the oxygen from migrating. The increase of
the migration space and the weakening of the local oxygen–cation bonds correspond to a decrease of the
migration barrier, and vice versa. The contribution of the bond strength to the changes in the migration
barrier is more significant than that of the opening of migration space in strained YSZ. A database of
migration barrier energies as a function of lattice strain for a set of representative defect distributions in
the vicinity of the migration path inYSZwas constructed. This database was used in kinetic Monte Carlo
(KMC) simulations to estimate the effective oxygen diffusivity in strained YSZ. The oxygen diffusivity
exhibits an exponential increase up to a critical value of tensile strain, or the fastest strain. This increase is
more significant at the lower temperatures. At the strain states higher than the critical strain, the
diffusivity decreases. This is attributed to the local relaxations at large strain states beyond a limit of
elastic bond strain, resulting in the strengthening of the local oxygen–cation bonds that increases the
migration barrier. The highest enhancement of diffusivity in 9%-YSZ compared to its unstrained state is
6.8 103 times at 4% strain and at 400 K. The results indicate that inducing an optimal strain state by
direct mechanical load or by creating a coherent hetero-interface with lattice mismatch can enable
desirably high ionic conductivity in YSZ at reduced temperatures. The insights gained here particularly
on the nonlinear and competing consequences of lattice strain on the local bonding structure and charge
transport process are of importance for tuning the ionic transport properties in a variety of solid-state
conducting material applications, including but not limited to fuel cells.
Introduction
There is much interest in decreasing the Solid Oxide Fuel Cell
(SOFC) operating temperature to an intermediate-to-low
temperature range to improve chemical and mechanical stability
and reduce the cost. Enhancing the ionic transport in the electrolyte
and cathode materials is a key for this purpose. Over the
past several decades, developing new compositions with high
ionic mobility at low temperatures has been the primary activity.1
In recent years, however, nanoscale structures were shown as
means for potentially increasing the transport properties in
oxides of relevance to SOFCs.2–9 The importance of utilizing the
favorable interface properties of nanoscale structures shifted the
paradigm from developing new material compositions to
tailoring nano-structures using existing materials for the above
purpose.5,10 Among these, nanoscale hetero-layered oxides were
suggested to enhance the ionic conductivity along the interfaces
by orders of magnitude.5–9 The presence of an interface between
dissimilar oxides is commonly the source of charged mobile
defects and local anisotropic distortions and strained interatomic
bonds. These contribute to changes in the local electronic
structure; formation, coordination and redistribution of charged
defects; and mobility of defects in the vicinity of interfaces. Most
recently, an eight orders of magnitude increase in ionic conductance
was reported by Barriocanal et al.8 for 1–30 nm-thick yttria
stabilized zirconia (YSZ) layers coherently ‘‘strained’’ between
dielectric SrTiO3 (STO) layers. The exact nature of the ionic vs.
electronic conductance induced in the vicinity of the YSZ/STO
interface remains debatable.11 While these observations imply
that the hetero-interfaces with controlled lattice-strain and
defect-chemistry can play an important role in improving the
ionic mobility, results to date remain phenomenological without
systematic theoretical efforts. Beyond the SOFC materials
domain, similar questions for the role of hetero-interfaces on the
local bonding structure and the charge transport process are
open for tuning the transport properties in a variety of solid-state
conducting material applications.12–15 Systematic experimental
and theoretical investigation of the role of local strain states on
tuning the mechanism and magnitude of charge transfer properties
in these materials is prone to pursuit.
Department of Nuclear Science and Engineering, Massachusetts Institute
of Technology, 77 Massachusetts Avenue, Cambridge, MA, 02139, USA.
E-mail: [email protected]; Tel: +617-324-4009
This journal is ª The Royal Society of Chemistry 2010 J. Mater. Chem., 2010, 20, 4809–4819 | 4809
PAPER Journal of Materials Chemistry Home-High impact applications, properties and synthesis of exciting new materials | Journal of Materials Chemistry
Two main mechanisms contribute to increase the ionic transport
at the hetero-interface: (1) a favorable strain state at the
interface to shift and/or change the symmetry of electron energy
levels to provide for improved charge transfer16 and ion
mobility,17 and (2) the alteration of the defect chemistry7,18 near
the interface to enhance the density and distribution of desired
charge carriers. To understand the underlying mechanism of
oxygen anion transport in anion conducting ceramics, a significant
number of theoretical studies using atomistic simulations
based on both density functional theory (DFT) and empirical
potentials have been conducted (examples related to this work in
ref. 19–25). However, most of them are focused on the ‘‘bulk’’ or
the ‘‘surface’’ ionic transport. Only few theoretical studies exist
for the anion transport at hetero-interfaces, on specific cases,
i.e. CaF2/BaF2, CeO2/YSZ,20 with no generalization of results.
Of particular interest for this paper is the effect of lattice strain
on the ionic conductivity in the YSZ, as may be encountered due
to a lattice mismatch near the YSZ/STO interface. In probing the
effect of strain on conductivity in YSZ, a number of experimental
and simulation studies have been conducted in the past.21–24 Even
though the past results imply a potential increase in the
conductivity, the magnitude of increase was far lower than what
was reported recently in ref. 8 for the YSZ/STO interface.
Furthermore, the reported increases in magnitude of the
enhancement are inconsistent among these past studies because
of the differences in conditions, structures and methods that were
considered. For example, results were based on examining either
incoherent hetero-interfaces which cannot sustain large lattice
strain,21 or low strain conditions,24 or molecular dynamics
simulation with empirical potentials23,24 which are not appropriate
for interfaces, or phenomenological modeling using
macroscopic elastic properties of the material.22 In essence, the
mechanism by which a biaxial lattice strain alters the anion
transport in the YSZ fluorite structure was not uncovered. Given
these past inconsistencies and open questions, we adopt a firstprinciples-
based approach in probing the mechanistic picture of
the ionic conduction in strained YSZ. Our specific objectives are
to elucidate the mechanism and to identify the possible extent of
increase in the oxygen ion conductivity in YSZ due to a biaxial
lattice strain.
In order to attain reliable results in probing the effect of strain
on ionic transport in YSZ, it is important to consider a realistic
representation of the defect structure, in particular the distribution
of vacancies and dopant cations in YSZ. This is because the
oxygen vacancy migration barrier, which can be altered by strain,
depends strongly on the local defect distributions.25 On the other
hand, it is difficult to characterize experimentally the vacancy
and cation distribution in YSZ at atomic resolution, therefore
this type of information is not available from literature
straightforwardly. In spite of the difficulties, past work on this
topic provided useful insights on the vacancy–vacancy, vacancy–
cation and cation–cation interactions, enabling a reasonable
description of the defect structure in YSZ.
Here we provide a summary of defect–defect interactions that
we considered for constructing the YSZ model. The final YSZ
model used in our work should reasonably capture these reported
defect–defect interactions. First of these interactions
considered, namely the vacancy–cation interaction, determines
the relative position of the vacancy with respect to the Zr and Y
in YSZ. In earlier studies for this interaction, there were contradicting
experimental results; whether the vacancy is first
nearest neighbor26–28 (1NN) or second nearest neighbor29–31
(2NN) to the Y cations. However, later experimental30–33 and
theoretical34–40 studies indicated that the vacancy favors to bind
to the host Zr cations, resulting in a configuration such that the
vacancy is 1NN to the Zr and 2NN to the Y cations. Using DFT
calculations, Bogicevic and Wolverton37 showed that the local
relaxation of atoms was responsible for the vacancy to favor the
2NN position to the Y cation; opposite to what is expected from
the electrostatic interactions alone. The same trend was reported
in doped CeO2 when the dopant radius is larger than the host
cation,41 as in YSZ. For the vacancy–vacancy and the cation–
cation interactions, they concluded that the vacancies locate
themselves as far away as possible from each other, and that the
1NN Y–Y pair configurations are slightly more favorable than
the 2NN Y–Y pairs.36,38 More recent experimental work using
selected area diffraction also indicated a random distribution of
oxygen vacancies, maximizing the distance between them, in up
to 10% yttria doped YSZ.42 Consistent with the findings of
Bogicevic and Wolverton, theoretical results by Predith et al.43
suggest even the existence of the ordering of the dopant cations
(Y) in YSZ at high doping concentrations, and this cation
ordering could exist locally in YSZ at the lower Y concentration
range. Among these three types of defect–defect interactions, the
vacancy–vacancy interaction is the strongest followed by the
vacancy–cation interaction, and the cation–cation interaction is
the weakest of the three.36 Supporting the importance of the
vacancy–vacancy interactions, it was shown recently that the
vacancy distribution around the migration path could significantly
impact the value of the migration barriers in YSZ.44
In this study, we considered only the strain state in the bulk
to provide for improved ionic mobility in YSZ and did not
include an interface explicitly as in Barriocanal’s work.8
Therefore, the model in this paper excludes the effect of altered
defect chemistry near the interfaces of YSZ as a source of
increase in oxygen diffusivity. In the following sections, we
describe the simulation approach to calculate the oxygen
vacancy migration path and energy barriers. We report the
energetically favorable distributions of oxygen vacancies and Y
cations in the YSZ model structure to enable a realistic representation
of YSZ in this study. The mechanism by which the
lattice strain increases or decreases the anion transport in YSZ
are discussed. Because the vacancy migration barrier depends
on the local defect configurations, a database of migration
barrier energies as a function of lattice strain for a set of
representative defect distributions in the vicinity of the migration
path are constructed and presented. Finally, the effective
tracer diffusivity of oxygen in YSZ at different strain states is
presented and interpreted in comparison to the recent results
reported for the YSZ/STO interface.8
Simulation approach
The YSZ model was established to capture the charged defect
interactions and the defect distribution in the bulk YSZ, as
summarized in the Introduction. Upon identifying the stable
vacancy–cation distribution, the effect of the biaxial lattice
strain in the xy-plane on the migration barrier and on the
4810 | J. Mater. Chem., 2010, 20, 4809–4819 This journal is ª The Royal Society of Chemistry 2010
effective vacancy diffusivity were calculated and interpreted.
We used first-principles-based simulations to directly and
accurately probe the atomistic-scale nature of the vacancy
migration in the presence of lattice strain. The kinetic Monte
Carlo (KMC) method was used for calculating the effective
diffusivities.
The first principles calculations for identifying the energies
were performed by the Vienna ab initio simulation package
(VASP),45,46 employing density functional theory (DFT) using
a plane-wave basis set. Projector-augmented wave (PAW)
method47,48 with plane waves up to the energy cutoff at 400 eV
was used and the exchange–correlation energy was evaluated
by the generalized gradient approximation (GGA) using the
Perdew–Burke–Ernzerhof (PBE) function.49 Y, Zr and O atoms
were described by 11 (4s24p65s24d1),12 (4s24p65s24d2),and 6
(2s22p4) valence electrons, respectively. The G point in the Brillouin
zone was selected in the calculation. We confirmed the
energy convergence by comparing the energy obtained from each
of the k-point meshes of 3 3 3, 5 5 5 and 7 7 7. The
calculated energies are within 0.6% error compared to the energy
obtained using only the G point. The same condition was applied
to all of the simulations in this study unless stated otherwise.
YSZ has the fluorite crystal structure with the cations
(Y3+, Zr4+) occupying the FCC lattice sites and the O2 ions
occupying the tetragonal sites. For every two Y atoms, one
oxygen vacancy is created to satisfy charge neutrality. The
oxygen vacancy migrates by hopping to the adjacent tetragonal
site as schematically illustrated in Fig. 1. The migration barrier
depends on the neighboring atoms 1–6 in the figure.25 The
surrounding cations within a larger cutoff radius to the migration
path must also be accounted for in order to accurately quantify
the effect of the strain on the migration barrier.50 However, this
leads to an impractical number of patterns (well above 10 000
configurations just by including the next nearest neighbor cations
to the migration path) to be handled by the DFT simulation.
Therefore, we selected a subset of defect structure patterns near
the migration path, forming a database of migration paths and
barriers. This database serves as input for the KMC calculations
in estimating the effective diffusivity as a function of lattice strain
and temperature. In constructing this database, first, we
considered only two combinations of atoms at the sites 3 and 4,
as Zr–Zr and Zr–Y, in the migration path shown in Fig. 1. The
migration barrier through the Y–Y pair was found to be significantly
higher than these two combinations, and this makes the
Y–Y pair unlikely on a vacancy migration path. This selection
was also validated by the work of Car et al.,50 which implied the
Zr–Zr and Zr–Y subset at sites 3 and 4 to be sufficient to capture
the quantitative trends in ionic conductivity in YSZ as a function
of Y2O3 doping concentration. Second, we considered the defect
configurations outside of the first neighbor cations to the
migration path to affect the barrier. The particularly important
impact of the vacancy interactions on the migration barrier was
recently reported by Pietrucci et al.44 Therefore, in creating the
migration path and barrier database, we varied the relative
positions of the Y cations and the vacancies with respect to the
initial and final positions of the vacancy on the migration path.
The details of the migration path configurations in this database
are provided in the following section. For each migration path
configuration, nudged elastic band (NEB) method51 was used to
find the vacancy migration barrier as a function of biaxial strain
on the xy-plane. Given that the vacancy migration path in YSZ is
rather simple with single hops between the adjacent tetragonal
sites, the saddle point is expected near the midst of the migration
path.25,50 Three images between the initial and the final states
were taken in the simulations. Here we used the climbing image
NEB method52 to calculate not only the minimum energy path
but also the transition state configuration at the saddle point.
The method drives the image with the highest energy up to the
saddle point, by maximizing its energy along the band connected
to the image and minimizing in all other directions. This way, the
method allows us to accurately calculate the saddle point energy
and configuration with a small number of images needed in the
calculation. The NEB simulations were conducted with fixed
lattice vectors.
The database created for the migration barrier values at
different lattice strains and the defect configurations was
provided as input for kinetic Monte Carlo (KMC) simulations.
KMC calculations were performed to quantify the effective
diffusivity of oxygen in YSZ as a function of strain and
temperature. The diffusivity calculated in this work is the oxygen
tracer diffusion constant. We expect that the lattice strain has
a similar order of magnitude impact on the mobility of an oxygen
ion in both the chemical and the tracer diffusion constants,
especially in the range of low vacancy concentrations.53 In the
simulation model, oxygen atoms were assumed to hop to the
nearest neighbor site on the oxygen sublattice through the path
shown in Fig. 1. The rate of hops in the migration paths were
obtained by the Boltzman relationship as,
vAB ¼ n0 expEAB
kBT (1)
where EAB, found by the DFT-NEB simulations, denotes the
migration barrier across the A–B cation pair at sites 3 and 4 in
Fig. 1. The attempt frequency was fixed at n0 ¼ 1013 s1, which is
appropriate for most metal oxide systems.50 The periodic
boundary conditions were applied to all three dimensions of the
simulation cell. The vacancy diffusivity, Dv was extracted from
mean square displacement hR2i calculated in the KMC simulation
using the Einstein relation,
hR2i ¼ 6Dvt. (2)
The oxygen diffusivity DO was then obtained using Dv
considering the balance between the fluxes of oxygen atoms and
vacancies as,
Fig. 1 Migration path of the oxygen vacancy. Positions 1–6 are occupied
by the cations Zr or Y.
This journal is ª The Royal Society of Chemistry 2010 J. Mater. Chem., 2010, 20, 4809–4819 | 4811
DO ¼
cv
1 cv
Dv (3)
where cv is the vacancy concentration fraction. Dv exponentially
depends on the effective barrier for oxygen vacancy migration,
Eeff, as,
Dv ¼ dv
0 exp
Eeff
kBT (4)
where, d0v
is a constant. At small values of cv, D0 scales almost
linearly with cv. On the other hand, the change in the migration
barrier, Eeff, contributes exponentially to Dv and D0.
Results and discussion
This section consists of five sub-sections, I–V. Sub-sections I–II
demonstrate that the YSZ model in this work captures the
reported vacancy and cation interactions in stabilizing the cubic
structure in YSZ. The results in I and II are important for
obtaining reliable values of the energy barriers for vacancy
migration in YSZ. Sub-section III reports the effect of defect
configurations in the vicinity of the migrating vacancy on the
migration energy barrier. This proves that, in order to assess the
effect of strain on migration barriers, one must consider a realistic
subset of defect configurations in YSZ, and we report these
in sub-section IV. Sub-section IV in particular illustrates the
mechanism by which strain alters the diffusion barriers for
oxygen vacancies in YSZ. Finally sub-section V quantifies the
effective increase in oxygen tracer diffusivity as a function of
biaxial lattice strain, based on the subset of defect configurations
considered in this work for YSZ.
I. Defect structure in YSZ: oxygen vacancy and yttrium
distribution
To determine the favorable location of vacancies with respect to
the Y cations in the YSZ model, we tested a number of defect
configurations. These models constitute 30 Zr, 2 Y and 63 O
atoms, corresponding to a low doping ratio, 3% Y2O3, allowing
only one vacancy in the simulation cell, in order to eliminate the
vacancy–vacancy interaction in this part of the simulations. At
this low Y2O3 doping, the YSZ oxygen sublattice cannot keep its
cubic structure.34 Therefore, only in this part of the simulations,
the atoms beyond the second neighbor cation shell to Y cations
and beyond the third neighbor oxygen shell to the vacancy were
fixed to retain the cubic phase.34 The comparison of the energies
for the different YSZ models and defect configurations is presented
in Fig. 2. The energies are shifted such that the lowest
energy found is zero. The results show that the total energy can
vary by up to 0.5 eV when the Y cation is located at the third or
first neighbor site to the vacancy, instead of the second. This
result is consistent with the prior reports on cation–vacancy
interactions in YSZ, as discussed in the Introduction. Since this
large energy difference arises from local configuration differences,
the result suggests that the cation configuration
surrounding the diffusion path of the oxygen–vacancy pair
should be taken into account for accurate evaluation of the
migration barriers.
II. Distortion from the cubic lattice sites in YSZ
The atoms are not stabilized at the exact fluorite structure lattice
centers in YSZ, and large displacements of the cations and the
anions from the fluorite lattice sites are present.55–57 This is in
part due to the relaxation of the nearest neighbor oxygen anions
and cations towards and away from the vacancy, respectively.32
For an acceptable YSZ structure in this study, it is important to
capture this deviation quantitatively consistently with past
experimental and simulation reports. For this purpose, we
evaluated the dependence of the atomic displacements from the
ideal fluorite lattice sites on the local defect structures in a set of
candidate YSZ models. The fully relaxed structure of the models
constituted 9% Y2O3 doped YSZ, made of 26 Zr, 6 Y and 61 O
atoms in the simulation cell. The relaxation of the atoms were not
constrained in this and the following simulations because the
cubic structure is stable in the 9% Y2O3 doped YSZ, as long as
Fig. 2 Energy as a function of the neighboring location of the oxygen
vacancy with respect to the Y cations in YSZ. The energies were calculated
for different cation/vacancy distributions, and plotted such that the
lowest energy found is zero.
Fig. 3 Fully relaxed atomic configurations of the 9% yttria doped
zirconia with (a) nearest neighbor Y–Y pairs, and (b) uniformly
distributed Y. Yellow, gray and red spheres indicate Y, Zr and O ions.
Large distortion of the oxygen anion sublattice is evident when Y cations
are uniformly spread in the YSZ without any short-range ordering of the
Y cations.
4812 | J. Mater. Chem., 2010, 20, 4809–4819 This journal is ª The Royal Society of Chemistry 2010
the vacancies and dopant cations are distributed to enable the
stability. Fig. 3 shows the relaxed atomic configurations for two
models of YSZ. The first one (Fig. 3(a)) has first nearest neighbor
(1NN) Y–Y pairs and the second one (Fig. 3(b)) has a uniform
distribution of Y cations in the cell. In both models, the vacancies
were placed as far away as possible from each other and located
second nearest neighbor (2NN) to Y cations if allowed by the
available space. The first model retains an overall cubic structure
with oxygen anions displaced by 0.1–0.27 A from the ideal
fluorite lattice sites. The extent of the displacements calculated
here is well within the experimentally measured value of
0.36 A28 and prior simulation results of 0.20–0.31 A.34 On the
other hand, the second model shows large distortions, exemplified
at the circled region in Fig. 3(b),and has 1.2 eV higher energy
compared to the first one. This difference arises due to local
configuration differences, and thus, can affect the vacancy
migration barrier in YSZ which ranges from 0.2 to 1.4 eV.25
The lack of stability of the cubic oxygen sublattice in the second
model is due to the limited available space to accommodate the
relative positions of the vacancy and Y cations at the favorable
distances when the Y atoms are distributed uniformly into the
model. In the circled region in Fig. 3(b),a vacancy is located at
the 1NN position to the Y cations. This is because the 2NN
position for the vacancy is not available to all the Y atoms when
the Y atoms are distributed uniformly in the model. The
restriction of some vacancies to the unfavorable 1NN positions
with respect to the Y in YSZ destabilizes the cubic phase. On the
contrary, the first model configuration with the Y–Y pairs
introduces more space for the vacancies to be placed at the 2NN
positions with respect to the Y cations. This enables more
stability on the cubic oxygen sublattice. This result is consistent
with the locally ordered nature of Y cations in YSZ at high
concentrations of Y2O3, as found by Predith et al.,43 even if the
YSZ has only 9% Y2O3, here. While the effect of Y–Y pairing on
the migration barrier in a cubic YSZ is small compared to that of
the vacancy–vacancy and vacancy–cation interactions,36 it is
important to consider it to enable a stabilized cubic fluorite
structure of YSZ in the model.
III. Vacancy migration barriers in YSZ: effect of the local
defect structure
As implied in the previous section, the local defect structure
around the oxygen vacancy migration path in YSZ affects the
stability of the structure and the migration barrier,54 and this in
turn can determine the favorable paths for vacancy to migrate
through. To illustrate this, we present the differences in the
migration barrier in a model with 1NN Y–Y pairs, which lower
the energy and enable the stabilization of cubic phase. The
consequent migration paths considered were A / B / C in
Fig. 4(a). All six nearest neighbor cations (1–6 as shown in Fig. 1)
are the same, Zr, for each migration step in A/B/C. On the
other hand, the distribution of vacancies and Y cations near A
/B and B/C differ beyond the nearest neighbor cations. This
enables to investigate the variation of the migration barrier as
a function of the defect structure near the migration path beyond
the nearest neighbor cations. Fig. 4(b) shows the relative energy
in the NEB trajectories for the migration through A/B/C.
The resulting migration barriers range from 0.55 to 0.95 eV, even
though all the cations in the nearest neighbor position to the
migration path are the same. This significant difference exemplified
in the migration barriers in A/B/C indicates that it is
important to systematically take into account the longer-range
interactions between the migrating vacancy and its surrounding
cations and vacancies in YSZ.
We evaluated a subset of defect configurations near the
migration path to form a representative database of migration
barriers in YSZ. The five different configurations, models A–E,
that yield ten forward and backward jump barriers considered
here are shown in Fig. 5. This subset aims to capture the isolated
Fig. 4 (a) The vacancy positions and local structure considered for the
vacancy migration path, A / B / C. Large black spheres indicate
positions of the vacancies. Yellow, gray and red spheres indicate Y, Zr
and O ions. Numbers of Y cations at the second nearest neighbor position
to the vacancy is two at sites A and B, and one at site C. (b) NEB barriers
calculated along the A / B and B / C paths. The variation in the
migration barrier is clear in (b),even though all the six nearest neighbor
cations in the migration path are Zr. The difference in barriers is due to
the dissimilar vacancy and Y distribution near the A / B and B / C
paths beyond the nearest neighbor cations.
Fig. 5 Models A–E with varying defect configurations surrounding the
migration path. ‘‘# of 2NNY’’ denotes the number of second nearest
neighbor Y cations to the vacancy positions. ‘‘NN vacancy position’’
denotes that the nearest neighbor vacancy to the migrating vacancy
is located at n-th nearest neighbor position on the anion sublattice
(n-th ¼ 3rd or 4th in these models). Lower panel indicates the 2NN Y
position, and the 3rd and 4th nearest neighbor vacancies on the anion
sublattice with respect to the migrating vacancy V.
This journal is ª The Royal Society of Chemistry 2010 J. Mater. Chem., 2010, 20, 4809–4819 | 4813
and coupled effects of vacancy–vacancy and vacancy–cation
interactions in the local environment around the migration path.
The vacancy–vacancy interaction is considered up to the fourth
nearest neighbor on the anion sublattice, which was the
maximum that could be accommodated in our 2 2 2 unit cell
model when the nearest neighbor vacancy pairs are avoided. The
cation–vacancy interaction was taken up to the second nearest
neighbor. First, for considering the effect of the nearest neighbor
cations in the migration path (sites 1–6 in Fig. 1),we only
included the variation of the bridging cation pairs (sites 3 and 4
in Fig. 1). This was shown to be the most important contribution
to the vacancy migration in YSZ50 when the long-range interactions
were not considered. In models A, B, and E, the vacancy
migrates through the Zr–Zr cation pair, and in models C and D,
through the Zr–Y cation pair in the midst of the path. All the
other cation sites closest to the migration path (sites 1, 2, 5 and
6 in Fig. 1) are Zr. The path through Y–Y pair was excluded from
the simulations since it has a much higher migration barrier
compared to the Zr–Zr and Zr–Y pairs,25,50,54 making it significantly
less likely to migrate through the Y–Y cation pair. Second,
for vacancy–vacancy interaction, models A and C consider the
position of nearest neighbor vacancies, denoted as ‘‘NN vacancy
position’’ in Fig. 5, to the initial and the final position of the
migrating vacancy. The distance between the migrating vacancy
and the nearest neighbor vacancy changes from the fourth to the
third on the anion sublattice upon migrating from the initial to
the final position (and vice versa in the backward migration). # of
2NN Y are the same in models A and C. Third, for vacancy–
cation interaction beyond the nearest neighbor cations, models B
and D consider the number of Y cations that are at the second
nearest neighbor (2NN) position with respect to the initial and
the final state of the vacancy. This is denoted as ‘‘# of 2NN Y’’ in
Fig. 5. The number of 2NN Y cations changes from two to three
upon migrating from the initial to the final state (and vice versa in
the backward migration). Lastly, in model E, both the number of
2NN Y cations and the distance to the nearest vacancies at the
initial and the final state were varied. The former changes from
three to two and the latter from the fourth to the third nearest
neighbor (and vice versa in the backward migration).
In creating the database of migration paths and barriers,
a model to represent a symmetric defect configuration around the
initial and final states of the migrating vacancy was not possible.
This is because of the computational limitation in the simulation
cell size used in our DFT/NEB calculations. In computing the
effective diffusivity in the following sections, we will show that
we represented the barrier for a symmetric migration path to be
the average of the barriers from models A–E.
A simulation unit cell of 2 2 2 dimensions, containing
26 Zr, 6 Y and 61 O atoms corresponding to a 9% yttria doped
zirconia was used in the calculation of the migration energy
barriers for each of the models A–E described in Fig. 5. The
oxygen vacancies and the Y cations were distributed in the cell
satisfying the defect interactions discussed in the previous
sections. Atomic configurations and the cell vectors were fully
relaxed such that the total energy was minimized and the net
stress acting on the cell was zero. The migration barriers for each
model path were calculated using the climbing NEB technique.
The calculated barriers for each path (A–E) are shown in Fig. 6.
The energy values are relative to the initial state configuration.
The reaction coordinate of each NEB image was normalized with
respect to the distance between the initial and the final state of the
vacancy. The oxygen vacancy migration process requires
climbing over the highest energy (the saddle point) between the
initial and the final state. E!B and E B in Fig. 6 are the migration
barriers for the forward and the backward jumps, respectively,
corresponding to the migration from the initial to the final
positions of the vacancy shown in Fig. 5, and vice versa.
As a result, the migration configurations in models C and D
have higher barrier than those in models A, B and E. This is
primarily because of the difficulty in passing through the Zr–Y
pair compared to that through Zr–Zr pair. The forward jump
barrier is higher than the backward jump barrier for models
A and C, and lower for models B and D. For models A and C,
this is because the distance from the migrating vacancy to the
nearest neighbor vacancy is shorter at the final state, making it
more difficult for the migrating vacancy to approach to the final
state, primarily due the electrostatic interaction of vacancies. For
models B and D, the number of Y cations at the second nearest
neighbor position to the migrating vacancy is larger at the final
state, making this final state energetically more favorable, and
thus the reduced migration barrier in the forward migration.
In model E, the forward jump barrier is higher than the backward
one. Considering the vacancy–cation interactions, the final
state in model E is energetically more favorable than the initial
state because of the presence of a larger number of 2NN Y
cations to the migrating vacancy at the final state. This suggests
the migration barrier for the forward jump to be lower than the
backward one. However, the result shows the opposite behavior,
because the migrating vacancy is actually getting closer to its
nearest neighbor vacancy at its final state. This indicates that the
vacancy–vacancy interaction dominates the vacancy–cation
interaction in determining the energy barriers in model E. All
implications related to the defect–defect interactions leading to
the resulting migration barriers in models A–E in the unstrained
state of YSZ are consistent with those suggested by previous
Fig. 6 NEB paths and the corresponding energies in the migration path
calculated for models A–E. The relative energy and the coordinates are
normalized with respect to the initial energy, and the distance between the
initial and the final state, respectively. The saddle point energy is the
highest energy state between the initial and the final states, and E!B
and E B stand for barriers in the forward and the backward jumps,
respectively.
4814 | J. Mater. Chem., 2010, 20, 4809–4819 This journal is ª The Royal Society of Chemistry 2010
reports by Bogicevic et al.36 Therefore, the results presented here
enable a reasonable model of YSZ for assessing the effect of
biaxial lattice strain on the ionic mobility.
IV. Migration barrier as a function of biaxial lattice strain
In assessing the role of the biaxial xy-plane lattice strain on the
oxygen vacancy transport, we investigated the corresponding
changes in the migration barriers, and the underlying reasons for
these changes on models A–E in Fig. 5. The xy-plane strains
considered here were 0.00, 0.02, 0.04, 0.06 and 0.08. The range of
strain considered in this study is far above the experimentally
measured fracture strain of a single-crystal YSZ that is less than
0.5%.58 However, large strain in brittle ceramics could be achieved
by creating coherent hetero-interfaces at the nanoscale.
A maximum 7% tensile lattice strain was suggested to exist in the
YSZ thin film interfaces coherently strained between STO layers
due to the mismatch between STO and YSZ.8 On the other hand,
ability of thin-film (10–30 nm) YSZ to sustain this large tensile
strain is unclear. We considered only tensile strain to increase the
interatomic distances in YSZ, allowing more space for the
oxygen migration. Results for model A are shown in Fig. 7.
The migration barrier monotonically reduces with the increase in
the biaxial lattice strain. The barrier in the forward jump
decreases from 0.47 eV at 3 ¼ 0.00 to 0.07 eV at 3 ¼ 0.04, suggesting
a significant enhancement of the vacancy diffusivity
through this path under strain. Beyond this strain level, the
vacancy was not stable at the initial site and favored to form at
the defined final site in model A, with no barrier in the suggested
transition path. This result implies that the defect interactions
influencing the stability and energetics can also change with the
strain state of YSZ. While not explicitly included in our DFT
calculations, we implicitly took into account the effect of the
strain on the stable defect configurations in the KMC simulations,
as explained in the following section. We considered both
the forward and the backward migration barriers, and the system
has a higher probability of taking the lower energy migration
paths favoring the more stable structures.
The picture is not as straight-forward in models B–E as in
model A that shows a monotonic reduction in the migration
barrier with strain. Depending on the extent of the lattice strain,
the consequent barrier either decreases or increases in the
configurations represented by models B–E, as shown in Fig. 8.
The lack of a straight-forward trend in the magnitude of
migration barrier as a function of lattice strain suggests that there
is more than one mechanism acting simultaneously upon strain
in YSZ. Here, we explain the competing effects of strain on the
vacancy migration barrier in YSZ from an atomistic and electronic
perspective. The oxygen vacancy migration process
involves bond-breaking and bond-making between the migrating
oxygen and the surrounding cations. Therefore, we investigated
the changes in the electronic density distribution of the valence
electrons describing the Zr, Y and O atoms near the migration
path to explain the local effects of strain in YSZ. We focused on
the two planes in the migration path as shown in Fig. 9. The
charge distribution on the ‘‘oxygen plane’’ in Fig. 9 correlates to
the ‘‘space available’’ for the oxygen, denoted as O, to migrate
along the prescribed path. The migration space is also related to
the cation-cation distance midst the migration path. On the other
hand, the ‘‘bonding plane’’ in Fig. 9 shows a measure of the
‘‘bonding strength’’ between the oxygen and the nearest cations
at the initial site, quantified in terms of the oxygen–cation ‘‘bond
thickness’’. It is this set of bonds, shown between the oxygen
denoted as O and the cation denoted as C, which are broken in
the initiation of oxygen diffusion. Both the migration space and
the bond strength change when YSZ is subjected to biaxial strain.
Fig. 10 illustrates an example in which the increase in strain
increases the open space, and also increases the O–C bonding
strength to the extent of local inelastic relaxation near the initial
state of oxygen. In this case, the lattice strain is beyond just
Fig. 7 (a) NEB paths and the corresponding energies in the migration
path calculated for model A. The relative energy and the migration
coordinates are normalized with respect to the initial energy, and the
distance between the initial and the final state, respectively. The saddle
point energy is the highest energy state between the initial and the final
states. (b) The vacancy migration barrier, EB, as a function of biaxial
lattice strain. The filled circles and the open circles indicate forward and
backward jump, respectively. The barrier decrease as the strain increases
in the migration path in model A.
Fig. 8 The vacancy migration barrier, EB, as a function of biaxial lattice
strain for models B–E (in a–d). The filled circles and the open circles
indicate the forward and the backward jump barriers, respectively.
This journal is ª The Royal Society of Chemistry 2010 J. Mater. Chem., 2010, 20, 4809–4819 | 4815
elastically extending the oxygen–cation bond length; the oxygen,
O, was actually attracted towards the cation C and broke its
bond from the opposite cation upon local relaxation of atoms. As
a result, a stronger binding of O to C arises, shown with the
thicker O–C bond in Fig. 10. These two consequences of strain
act oppositely on the migration barrier – the former (opening of
migration space) favors the migration and the latter (local
inelastic relaxations near the oxygen) inhibits the migration.
A specific illustration of these competing consequences of
biaxial strain in YSZ represented by model E is shown in Fig. 11.
The resulting change in the migration space was quantified as the
distance between the cations midst the migration path, D(Zr–Zr)
in this configuration. The bond O–C in Fig. 11(b) is one of the
oxygen–zirconium bonds that break upon the migration of
oxygen. The large reduction in the forward migration barrier
from 0.51 eV at 3 ¼ 0.00 to 0.10 eV at 3 ¼ 0.06 as shown in
Fig. 8(d) corresponds to the increase in the D(Zr–Zr) distance
(shown in Fig. 11(a)) and the weakening of the O–C bond (shown
in Fig. 11(b)). On the other hand, this trend is not monotonic. At
higher strain, 3 ¼ 0.08, the cation–cation distance decreases and
the oxygen–zirconium bond strengthens again, resulting in
a significant increase of the migration barrier from 0.10 eV at
3 ¼ 0.06 to 0.47 eV at 3 ¼ 0.08. At this strain, the local relaxation
of the atoms particularly re-strengthens the oxygen-zirconium
bond at the initial site of the oxygen, also affecting the migration
space. If assumed independent, the oxygen–cation bond strength
is more strongly correlated to this large increase in the barrier at
large strain seen in Fig. 8(d). The increase in the cation–oxygen
bond strength is evident at 3 ¼ 0.08 in Fig. 11(b),but a very small
change occurs in the D(Zr–Zr) in Fig. 11(a). While we present this
set of results explicitly for the model E, the same qualitative
correlation of the competing effects of the migration space and
the cation–oxygen bond strength governed the migration barrier
changes with strain consistently also in the other models A–D.
V. Effective oxygen diffusivity as a function of biaxial lattice
strain
Models A–E constituted the database of the migration paths and
barriers as a function of the biaxial strain in YSZ. Using this
database to govern the overall oxygen diffusion, the macroscopic
oxygen diffusivity in the bulk of YSZ was estimated by KMC
simulations as described in the Simulation Approach section.
A 15 15 15 simulation cell of 8% yttria doped zirconia was
used, consistent with the doping ratio in YSZ reported in the
Ref. 8 for the YSZ/STO heterolayers. Y cations and vacancies
were distributed randomly in the KMC simulation cell. At a first
sight, this configuration contradicts with the implications for the
presence of 1NN Y–Y pairs in YSZ as discussed in Introduction.
On the other hand, the rules in selecting the diffusion barriers in
the KMC simulations were taken consistently with the local
defect configurations, including the effect of Y–Y pairs as evaluated
in the DFT-NEB calculations in this work. Therefore, we
considered this model in KMC to be reasonable for assessing the
overall oxygen diffusivity in YSZ. Supporting this approach,
past DFT and KMC studies of oxygen diffusivity in YSZ with
randomly distributed Y atoms and vacancies were shown to give
Fig. 9 The electronic density distribution on the ‘‘oxygen plane’’ and the
‘‘bonding plane’’ at the initial and the saddle states. The migrating
oxygen, O, moves towards the open migration space, and the bond
between oxygen O and cation C breaks in the migration.
Fig. 10 Change in the electronic density distribution on the oxygen
plane and the bonding plane under biaxial lattice strain, 3, from 0.00 to
0.07. Gray and red spheres represent Zr and O atoms, respectively.
4816 | J. Mater. Chem., 2010, 20, 4809–4819 This journal is ª The Royal Society of Chemistry 2010
the correct yttria doping ratio for the maximum oxygen
diffusivity.25,50
The effect of the local defect environment around the migration
path on the barriers input to the KMC simulations was
taken into account as follows:
Definitions. E!a and E a are the forward and backward jump
energy barriers in each model a (a ¼ A, B, C, D, E),respectively.
ri
v and rf
v stand for the distance of the first nearest neighbor (1NN)
vacancy to the migrating vacancy at the initial and final state,
respectively. ni
Y, and nf
Y are the number of second nearest
neighbor (2NN) Y cations to the migrating vacancy at the initial
and final state, respectively.
Rules. The rules in determining the energy barrier for each hop
in the KMC simulations treat the local defect configurations in
a coarse grained fashion, based on the models A–E (Fig. 5)
assumed to govern the oxygen vacancy diffusion in YSZ.
In models A and C, the vacancy forward hopping (from initial to
final state) results in a decrease of the 1NN vacancy distance to
the migrating vacancy, that is ri
v > rf
v. Then, in the KMC simulation,
E!a and E a (a ¼ A, C) are chosen for the cases of ri
v > rf
v,
and ri
v<rf
v, respectively. In models B and D, the vacancy forward
hopping (from initial to final state) results in an increase in the
number of 2NN Y to the migrating vacancy, that is ni
Y < nf
Y.
Then, in the KMC simulation, E!a and E a (a ¼ B, D) are chosen
for the cases of ni
Y <nf
Y and ni
Y > nf
Y, respectively. The explicit
definition of the rules implemented in the KMC simulations is as
follows:
(1) For the vacancy migration through a Zr–Zr pair midst the
migration path, we select the barrier among E!A, E A E!B and E B
according to the following criteria:
If ri
v > rf
v, select E!A; and if ri
v < rf
v, select E A.
If ri
v ¼ rf
v and ni
Y > nf
Y, select E!B; or if ri
v ¼ rf
v and ni
Y < nf
Y,
select E B.
If ri
v ¼ rf
v and ni
Y ¼ nf
Y, randomly choose barriers from E!A,
E A, E!B, and E B.
(2) For the vacancy migration through a Zr–Y pair midst the
migration path, we select the barrier among E!C, E C, E!D and E D
according to the following criteria:
If ri
v > rf
v, select E!C; if ri
v < rf
v, select E C.
If ri
v ¼ rf
v and ni
Y > nf
Y, select E!D; or if ri
v ¼ rf
v and ni
Y < nf
Y,
select E D.
If ri
v ¼ rf
v and ni
Y ¼ nf
Y, randomly choose barriers from E!C,
E C, E!D and E D.
(3) Vacancy migration through Y–Y pair is prohibited.
When ri
v ¼ rf
v and ni
Y ¼ nf
Y, the migration path is surrounded by
a symmetric distribution of defects. This case was not possible to
account for in the DFT-NEB calculations due to the computational
limitations on the simulation size. We assumed such
a symmetric path to be governed by a migration barrier that is
the average of the barriers from the non-symmetric models represented
in A–E here. Therefore, the random selections among
the entire set of E!a and E a (for a ¼ A, B, C, D, E) during the 108
KMC steps should provide sufficient statistics to represent the
average migration barrier for the symmetric path, ri
v ¼ rf
v and
ni
Y ¼ nf
Y. Moreover, the probability of having symmetric path in
the KMC runs was less than 3%. Therefore, the error that may
arise due to selecting an average migration barrier for the
symmetric configuration is expected to have little significance on
the effective diffusivity.
As seen from the rules above, the energy barrier for each hop
in the KMC simulations was chosen based on a coarse grained
treatment of the local defect configurations around the migrating
vacancy. We believe this approach gives a reasonable result,
particularly based on the validity of the models A–E (Fig. 5) to
govern the oxygen vacancy diffusion in YSZ. On the other hand,
the reader can refer to a more accurate treatment of the diffusion
barriers in KMC simulations using the cluster expansion
method,43,59 which can account for the energy deviations due to
the local atomic structures around each hop. In that case, the
KMC model itself has to take into account the stable distributions
of vacancies and cations (as in subsections I and II),rather
than a random distribution, making the task exhaustive.
The KMC simulations were performed as a function of
temperature ranging as 400, 600, 800 and 1000 K, and of biaxial
strain ranging as 0.00, 0.02, 0.04, 0.06 and 0.08. Ten simulations
with different distributions of dopant cations and vacancies were
conducted for each temperature and strain state. The resulting
strain dependence of the effective oxygen tracer diffusivity DO is
shown in Fig. 12. y-axis in the figure is on a relative scale with
Fig. 11 (a) Change in the cation–cation distance, D(Zr–Zr),and (b) the
bond thickness (strength) for O–C in the migration under biaxial lattice
strain from 0.00 to 0.08 for the forward jump in model E (see also
Fig. 8(d)). Increase in the migration space D(Zr–Zr),and the weakening
of the O–C bond correspond to the decrease of the migration barrier in
Fig. 8(d),and vice versa. Large gray spheres and small red spheres are Zr
and O atoms, respectively. The isoline interval is 0.5 A3.
This journal is ª The Royal Society of Chemistry 2010 J. Mater. Chem., 2010, 20, 4809–4819 | 4817
respect to the oxygen diffusivity, D0
O at 3 ¼ 0.00 for the lowest
temperature simulated, 400 K. The reported results show the
average values and the error bars over the ten simulations. The
result is not a simple monotonic increase of diffusivity with
increasing strain in the entire range considered. The oxygen
diffusivity exhibits an exponential increase up to a critical value,
4% of tensile strain. This increase is more significant at the lower
temperatures. The maximum enhancement in diffusivity
compared to the unstrained state in YSZ is 6.8 103 times higher
at 4% strain and at 400 K. At the higher strain states, the
diffusivity decreases. This is attributed to the local relaxations
beyond a limit of elastic bond strain, resulting in the strengthening
of the local oxygen–cation bonds that increases the
migration barrier.
Our first-principles-based simulation results indicate that
a biaxial strain in YSZ induced by a lattice-mismatch coherent
interface can significantly increase the diffusivity up to a critical
strain state, and worsen it beyond that strain state. One could
view this critical strain state as the optimal strain to induce in the
material to attain the maximum benefit of the biaxial lattice
strain on the oxygen diffusivity. We note that a recent prediction
by Schichtel et al.22 using a model based on the elastic properties
of YSZ estimated 2.5 orders of magnitude increase in the ionic
conductivity of YSZ at 7% strain at 573 K. The model of
Schichtel et al. relies on the isotropic pressure induced by the
strain to enhance the diffusivity, and does not take into account
the competing relaxations at these large tensile strain states in
YSZ. Therefore, it cannot capture the local bond-strength
changes that may mitigate the diffusivity at large strain values.
In spite of the differences in our and Schichtel’s model, the
magnitude of the maximum increase in the oxygen diffusivity in
strained YSZ predicted by these models range from 2.5 to 3.8
orders of magnitude. Therefore, the lattice strain (elastic and
inelastic) alone cannot explain the eight orders of magnitude
increase of ionic conductivity in YSZ reported by Barriocanal
et al.8 These results do not take into account the space-charge
effects in altering the anion concentrations near the YSZ interface.
On the other hand, the contribution to the increase in the
anion diffusivity due to a possible space-charge-induced increase
in the vacancy concentration is almost linear at low concentration
of vacancies (see eqn (3)) as in the 9% Y2O3 doped YSZ.
Therefore, this contribution is expected to be small compared to
the exponential increase found here due to the applied biaxial
lattice strain.60
Conclusion
This work provides a mechanistic and systematic assessment of
the effect of biaxial lattice strain on the oxygen diffusivity in
YSZ, independently from the defect chemistry near a heterointerface
with a dissimilar oxide, e.g. YSZ/STO as reported by
Barriocanal et al.8 Using DFT/NEB simulations, we uncovered
the underlying microscopic mechanism by which the biaxial
tensile strain acts on the oxygen vacancy transport in YSZ. The
migration barrier correlates fundamentally with two competing
processes acting in parallel in the presence of lattice strain: (1) the
migration space, quantified as the cation–cation distance midst
the migration path for the oxygen, and (2) the strength of the
bond between the oxygen and the nearest cation preventing the
oxygen from migrating. Increasing of the migration space and
weakening of the oxygen–cation bond contribute to the decrease
of the vacancy migration barrier, and vice versa. In assessing the
effective diffusivity as a function of strain in YSZ, we constructed
a data base of migration paths and barriers to account for the
important effects of local distributions of vacancies and dopant
cations near the migration path. KMC simulations using this
database showed that the macroscopic oxygen diffusivity
increases exponentially in YSZ up to a critical value of biaxial
tensile strain, with a more significant impact at the lower
temperatures. Beyond this critical value, the strain state in YSZ
worsens the diffusivity due to the local relaxations that
strengthen the oxygen–cation bond at the initiation of migration.
One could view this critical strain state as the optimal strain, or
the fastest strain, to attain the highest benefit of the biaxial lattice
strain on the oxygen diffusivity in the material system. The
maximum enhancement in diffusivity in YSZ was 6.8 103 times
higher at 4% strain and temperature of 400 K. This is a significant
and desirable increase. Nevertheless, it is far below the eight
orders of magnitude increase reported by Barriocanal et al. for
the thin-film YSZ layers coherently strained between STO
layers.8 Therefore, we believe that the lattice strain alone cannot
be responsible for the colossally high ionic conductivity suggested
to occur in YSZ in ref. 8. Nevertheless, given our results,
controlling the lattice strain by direct mechanical load or by
creating a coherent hetero-interface to induce an optimal strain
state can enable desirably high ionic conductivity at reduced
temperatures in the fluorite material, the YSZ in this case.
The insights gained here for the role of lattice strain on the
local bonding structure and ionic charge transport process are of
importance for modulating the transport properties in a variety
of solid state conducting material applications, beyond the
specific SOFC electrolyte material studied in this work. Interesting
examples may extend to Li-ion batteries, oxygen gas
sensors, semi-conductors and photovoltaics. We believe that the
mechanistic results obtained in this study, particularly on the
Fig. 12 The relative enhancement of the oxygen diffusivity calculated by
the kinetic Monte Carlo simulations at different temperatures and strain
states. D0
O is the oxygen diffusivity at 3 ¼ 0.00 at 400 K. The diffusivity
increases exponentially with strain, and shows a maximum at (or near)
3 ¼ 0.04 followed by a decrease at higher strains. The increase in diffusivity
is more significant at lower temperatures.
4818 | J. Mater. Chem., 2010, 20, 4809–4819 This journal is ª The Royal Society of Chemistry 2010
non-linear and competing consequences of strain on the bonding
structure in the vicinity of the ionic carrier species, provides
insights for creating functional strained nano-structures for ionic
and electronic conducting devices, including but not limited to
SOFCs.
Acknowledgements
We thank the Nuclear Regulatory Commission for financial
support, and the National Science Foundation for computational
support through the TeraGrid Advanced Support
Program under Grant No. TG-ASC090058.
References
1 J. C. Boivin and G. Mairesse, Chem. Mater., 1998, 10, 2870.
2 T. Suzuki, I. Kosacki and H. U. Anderson, Solid State Ionics, 2002,
151, 111.
3 I. Kosacki, C. M. Rouleau, P. J. Becher, J. Bentley and
D. H. Lowndes, Solid State Ionics, 2005, 176, 1319.
4 J. H. Shim, C.-C. Chao, H. Huang and F. B. Printz, Chem. Mater.,
2007, 19, 3850.
5 J. Maier, Nat. Mater., 2005, 4, 805.
6 N. Sata, K. Eberman, K. Eberl and J. Maier, Nature, 2000, 408, 946.
7 X. X. Guo, I. Matei, J.-S. Lee and J. Maier, Appl. Phys. Lett., 2007,
91, 103102.
8 J. G. Barriocanal, A. R. Calzada, M. Varela, Z. Sefrioui, E. Iborra,
C. Leon, S. J. Pennycook and J. Santamaria, Science, 2008, 321, 676.
9 S. Ramanathan, J. Vac. Sci. Technol., A, 2009, 27, 1126.
10 A. J. Jacobson, Chem. Mater., 2010, 22, 660.
11 X. Guo, Science, 2009, 324, 465; J. G. Barriocanal, A. R. Calzada,
M. Varela, Z. Sefrioui, E. Iborra, C. Leon, S. J. Pennycook and
J. Santamaria, Science, 2009, 324, 465.
12 M. P. O’Callaghan, A. S. Powell, J. J. Titman, G. Z. Chen and
E. J. Cussen, Chem. Mater., 2008, 20, 2360.
13 Y. Hoshina, K. Iwasaki, A. Yamada and M. Konagai, Jpn. J. Appl.
Phys., 2009, 48, 04C125.
14 J. Munguıa, G. Bremond, J. M. Bluet, J. M. Hartmann and
M. Mermoux, Appl. Phys. Lett., 2008, 93, 102101.
15 V. R. D’Costa, Y.-Y. Fang, J. Tolle, J. Kouvetakis and J. Menendez,
Phys. Rev. Lett., 2009, 102, 107403.
16 M. Mavrikakis, P. Stoltze and J. K. Nørskov, Catal. Lett., 2000, 64,
101.
17 J. Kilner, Nat. Mater., 2008, 7, 838.
18 J. Maier, Prog. Solid State Chem., 1995, 23, 171.
19 A. K. Ivanov-****z, Crystallogr. Rep., 2007, 52, 129.
20 D. C. Sayle, J. A. Doig, S. C. Parker and G. W. Watson, Chem.
Commun., 2003, 1804.
21 C. Korte, A. Peters, J. Janek, D. Hesse and N. Zakharov, Phys. Chem.
Chem. Phys., 2008, 10, 4623.
22 N. Schichtel, C. Korte, D. Hesse and J. Janek, Phys. Chem. Chem.
Phys., 2009, 11, 3043.
23 K. Suzuki, M. Kubo, Y. Oumi, R. Miura, H. Takaba, A. Fahmi,
A. Chatterjee, K. Teraishi and A. Miyamoto, Appl. Phys. Lett.,
1998, 73, 1502.
24 W. Araki, Y. Imai and T. Adachi, J. Eur. Ceram. Soc., 2009, 29, 2275.
25 R. Pornprasertsuk, P. Ramanarayanan, C. B. Musgrave and
F. B. J. Prinz, J. Appl. Phys., 2005, 98, 103513.
26 W. L. Roth, R. Wong, A. I. Goldman, E. Canova, Y. H. Kao and
B. Dunn, Solid State Ionics, 1986, 18–19, 1115; H. Morikawa,
Y. Shimizugawa, F. Marumo, T. Harasawa, H. Ikawa, K. Tohji
and Y. Udagawa, J. Ceram. Soc. Jpn, 1988, 96, 253; M. H. Tuilier,
J. Dexpert-Ghys, H. Dexpert and P. Lagarde, J. Solid State Chem.,
1987, 69, 153.
27 M. Weller, Z. Metallkd., 1993, 84, 6.
28 D. Steele and B. E. F. Fender, J. Phys. C: Solid State Phys., 1974, 7, 1.
29 B. W. Veal, A. G. McKale, A. P. Paulikas, S. J. Rothman and
L. J. Nowicki, Physica B+C, 1988, 150, 234; D. Komyoji,
A. Yoshiasa, T. Moriga, S. Emura, F. Kanamaru and K. Koto,
Solid State Ionics, 1992, 50, 291.
30 C. R. A. Catlow, A. V. Chadwick, G. N. Greaves and L. M. Moroney,
J. Am. Ceram. Soc., 1986, 69, 272.
31 P. Li and J. E. Penner-Hahn, Phys. Rev. B: Condens. Matter, 1993, 48,
10063; P. Li and J. E. Penner-Hahn, Phys. Rev. B: Condens. Matter,
1993, 48, 10074; P. Li and J. E. Penner-Hahn, Phys. Rev. B:
Condens. Matter, 1993, 48, 10082.
32 J. P. Goff, W. Hayes, S. Hull, M. T. Hutchings and K. N. Clausen,
Phys. Rev. B, 1999, 59, 14202.
33 K. Kawata, H. Maekawa, T. Nemoto and T. Yamamura, Solid State
Ionics, 2006, 177, 1687.
34 G. Stapper, M. Bernasconi, N. Nicoloso and M. Parrinello, Phys.
Rev. B: Condens. Matter Mater. Phys., 1999, 59, 797.
35 A. Bogicevic, C. Wolverton, G. M. Crosbie and E. B. Stechel, Phys.
Rev. B: Condens. Matter Mater. Phys., 2001, 64, 014106.
36 A. Bogicevic and C. Wolverton, Phys. Rev. B: Condens. Matter
Mater. Phys., 2003, 67, 024106.
37 A. Bogicevic and C. Wolverton, Europhys. Lett., 2001, 56, 393.
38 X. Xia, R. Oldman and R. Catlow, Chem. Mater., 2009, 21, 3576.
39 M. Sakib Khan, M. Saiful Islam and D. R. Bates, J. Mater. Chem.,
1998, 8, 2299.
40 M. O. Zacate, L. Minervini, D. J. Bradfield, R. W. Grimes and
K. E. Sickafus, Solid State Ionics, 2000, 128, 243.
41 D. A. Andersson, S. I. Simak, N. V. Skorodumova, I. A. Abrikosov
and B. Johansson, Proc. Natl. Acad. Sci. U. S. A., 2006, 103, 3518.
42 S. Garcıa-Martın, D. P. Fagg and J. T. S. Irvine, Chem. Mater., 2008,
20, 5933.
43 A. Predith, G. Ceder, C. Wolverton, K. Persson and T. Mueller, Phys.
Rev. B: Condens. Matter Mater. Phys., 2008, 77, 144104.
44 F. Pietrucci, M. Bernasconi, A. Laio and M. Parrinello, Phys. Rev. B:
Condens. Matter Mater. Phys., 2008, 78, 094301.
45 G. Kresse and J. Hafner, Phys. Rev.B: Condens. Matter, 1993, 47, 558.
46 G. Kresse and J. Furthm€uller, Phys. Rev. B: Condens. Matter, 1996,
54, 11169.
47 P. E. Bl€ochl, Phys. Rev. B: Condens. Matter, 1994, 50, 17953.
48 G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater.
Phys., 1999, 59, 1758.
49 J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77,
3865.
50 R. Krishnamurthy, Y.-G. Yoon, D. J. Srolovitz and R. Car, J. Am.
Ceram. Soc., 2004, 87, 1821.
51 G. Mills and H. Jonsson, Phys. Rev. Lett., 1994, 72, 1124.
52 G. Henkelman, B. P. Uberuaga and H. Jonsson, J. Chem. Phys., 2000,
113, 9901.
53 A. Van der Ven, G. Ceder, M. Asta and P. D. Tepesch, Phys. Rev. B:
Condens. Matter Mater. Phys., 2001, 64, 184307.
54 M. Martin, J. Electroceram., 2006, 17, 765.
55 F. Frey, H. Boysen and I. K- Bischoff, Z. Kristallogr., 2005, 220, 1017.
56 D. N. Argyriou, M. M. Elcombe and A. C. Larson, J. Phys. Chem.
Solids, 1996, 57, 183.
57 N. Ishizawa, Y. Matsushima, M. Hayashi and M. Ueki, Acta
Crystallogr., Sect. B: Struct. Sci., 1999, 55, 726.
58 D. Baither, M. Bartsch, B. Baufeld, A. Tikhonovsky, M. R€uhle and
U. Messerschmidta, J. Am. Ceram. Soc., 2001, 84, 1755.
59 A. Van der Ven, J. C. Thomas, Q. C. Xu, B. Swoboda and
D. Morgan, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 78,
104306.
60 A. Kushima and B. Yildiz, ECS Trans., 2009, 25, 1599.
This journal is ª The Royal Society of Chemistry 2010 J. Mater. Chem., 2010, 20, 4809–4819 | 4819
heres the first one..more to come..
