Background

Most of the past interest in pure ceria (CeO$_2$) and doped ceria materials have involved their applications in the automotive industry. Ceria has been studied extensively for its use in catalysis and oxygen sensing in automobile exhaust/emission systems [1]. It has also gained recent attention as a candidate for fuel cell electrodes [1] because of its high electrical conductivity and high thermal stability. These studies have provided significant insight into the understanding of material properties and defect chemistry of ceria and doped ceria materials.

Material Properties and Defect Chemistry of Ceria

Cerium(IV) oxide (CeO$_2$) crystallizes in the fluorite crystal structure with lattice constant $a=5.41134(12)$ Å. The fluorite structure consists of a face-centered cubic (f.c.c.) unit cell of cations with anions occupying the octahedral interstitial sites. This can also be seen as a superposition of an f.c.c. lattice of cations (Ce$^{4+}$) with lattice constant $a$, and a simple cubic (s.c.) lattice of anions (O$^{2-}$) with lattice constant $a/2$. In this structure (shown schematically in Figure 2.1), each cerium cation is coordinated by eight nearest-neighbor oxygen anions, while each oxygen anion is coordinated by four nearest-neighbor cerium cations.

Figure 2.1: Schematic of CeO$_2$ unit cell.

Fluorite structure oxides exhibit similar material properties, such as high radiation tolerance and high thermal stability. Cerium oxide is attractive as a surrogate for uranium oxide because it has the same fluorite crystal structure and many similar material properties, including melting temperature (UO $_2:\sim2870 \ensuremath{\,^{\circ}\mathrm{C}}$, CeO $_2: \sim 2600 \ensuremath{\,^{\circ}\mathrm{C}}$) and thermal diffusivity, and has been well characterized experimentally up to $700 \ensuremath{\,^{\circ}\mathrm{C}}$ [2][3][4][5][6].

While ionic conductivity is believed to be negligible in pure ceria, it increases significantly when ceria is doped with an aliovalent oxide like Y$_2$O$_3$ and La$_2$O$_3$. The open structure of the fluorite lattice is able to tolerate the high level of atomic disorder that would be introduced by this type of doping. When ceria is doped with a trivalent ion like Lanthanum which forms La$_2$O$_3$, a local charge imbalance is created. The lattice must compensate for this excess negative charge using one of three mechanisms: vacancy compensation, dopant interstitial compensation, and cerium interstitial compensation [7]. The mechanisms can be represented in Kröger-Vink notation:

$$M_2O_3 \overset{CeO_2}{\longleftrightarrow} 2M_{Ce}' + V_O^{\cdot\cdot} + 3 O_O^\times$$ (2.1)

$$2M_2O_3 \overset{CeO_2}{\longleftrightarrow} 3M_{Ce}' + M_i^{\cdot\cdot\cdot} + 6 O_O^\times$$ (2.2)

$$Ce_{Ce}^\times + 2 O_O^\times + 2M_2O_3 \overset{CeO_2}{\longleftrightarrow} 4M_{Ce}' + Ce_i^{\cdot\cdot\cdot\cdot} + 8 O_O^\times$$ (2.3)

In vacancy compensation (Equation 2.1), an anion vacancy is produced for every two dopant ions placed on the host cation sites. In dopant interstitial compensation (Equation 2.2), one dopant cation is placed in on interstitial site for every three dopant cations that are placed on the host cation sites. In cerium interstitial compensation (Equation 2.3), a cerium Frenkel pair is produced and the displaced cerium cation is placed on an interstitial site for every four dopant cations placed on the host cation sites. Empirical calculations performed by Minervini et al. show that for large dopant cations (cation radius $> 0.8$ Å), vacancy compensation is the preferred charge compensation mechanism [7]. Blumenthal et al. also ruled out the formation of interstitial as the compensation mechanism by measuring true density and comparing it with calculated values [8]. In terms of dopant solute concentration, the vacancy compensation mechanism can be represented by:

$$xMO_{1.5} + (1-x)CeO_2 \longleftrightarrow xM_{Ce}' + 0.5xV_{O}^{\cdot\cdot} + (1-x)Ce_{Ce} + (2-0.5x)O_O$$ (2.4)

This reactions implies that when $x/2$ moles of dopant oxide (M$_2$O$_3$) are added, cation sites are filled with $x$ moles of dopant cation (M$^{3+}$) and $(1-x)$ moles of host cation (Ce$^{4+}$), while anion sites are filled with $(2-x/2)$ moles of host anion (O$^{2-}$) and $x/2$ moles of anion vacancies (V $_O^{\cdot\cdot}$). Therefore, a predictable concentration of oxygen vacancies can be introduced into the crystal by controlling the concentration of the lanthanum dopant added. This process can be used to create an oxygen vacancy environment that is similar hypostoichiometric configuration of pure ceria.

Previous Computational Work on Ceria related Materials

Simulation techniques have been used to study a wide range of material properties over a wide range of length and time scales, including thermodynamic properties, defect structure and clustering, defect clustering, and transport phenomena. Techniques ranging from Density Functional Theory (DFT), a quantum mechanical theory used to study electron structure over a couple of atoms for a matter of picoseconds, to Finite Element Method (FEM), which can used to study structures as large as a reactor core over a matter of years, each have their own applications and advantages/disadvantages. Figure shows the classic illustration of the application of various modeling techniques of the vast time and length scales of interest.

Figure 2.2: Visualization of computational modeling techniques shown over the length and time scales that they apply to.

This study links the results from several different time/length scales. At the lowest level, the starting point, are the potentials for the La-doped ceria system that have been developed, either through DFT calculations or from fitting to experimental data. These potentials are fed into Molecular Dynamics (MD) simulations to calculate local configuration-dependent oxygen vacancy migration energies. These migration energies are then fed into KMC simulations to calculate oxygen diffusivity.

Since its development in the 1970s, Molecular Dynamics has been widely used to predict various material properties. However, while it has been used extensively in metal and metal alloy systems, its use in ceramic oxide systems has been comparatively lacking. This is likely due to a lack of confidence in available interatomic potentials. While the potentials in metal and metal alloy systems have been well developed and validated by experimental results, the potentials for ceramic alloys generate results that are not always consistent with experimental data.

Since the the understanding of atomic level interactions is very important in nuclear fuel research, a series of studies have been done to model and understand the thermodynamics and defect chemistry of uranium oxide (UO${_2}$) [9][10][11][12][13][14][15][16][17][18]. An extensive literature review was provided by Grovers et al. [19][20], where all of the available interatomic potentials on UO$_2$ were compared. They found that none of the potentials reviewed could predict all of the thermodynamic, mechanical, and defect transport properties. For example, they found that oxygen vacancy migration energies obtained from simulations ranged from $0.1$ eV to $0.7$ eV, while the experimentally measured energy was found to be $0.5$ eV. Similarly, oxygen interstitial migration energies obtained from simulations ranged from $0.1$ eV to $3.6$ eV, while the experimentally measured energies ranged from $0.9$ eV to $1.3$ eV. These discrepancies could be attributed to the fact that all of the potentials reviewed were obtained by fitting parameters to experimentally measured properties, such as lattice constant, lattice energy, dielectric constants, etc. As most of these potentials were fitted without considering defect energetics, their inability to accurately predict these values is possibly unavoidable. This survey also demonstrated the importance of testing and validating potentials before using them for predictive purposes.

While there are nineteen interatomic potentials available for UO$_2$, there are only a few potentials available for CeO$_2$, and even fewer with dopant parameters. Through an extensive literature survey conducted by Dr. Di Yun, only eight potentials for CeO$_2$ were found. From those eight, three were selected for comparison in this study. Similar to what was found Grovers' literature review for UO$_2$ [19][20], all of the potentials found for CeO$_2$ were in two forms.

The first potential form consists of the addition of a Buckingham term to the basic Coulomb potential. This form can be described by:

$$V_{ij}\left(r\right) = \frac{q_i q_j e^2}{4\pi\epsilon_0} + A_{ij}\exp{\left(-\frac{r}{\rho_{ij}}\right)} - \frac{C_{ij}}{r^6}$$ (2.5)

where $V_{ij}$ is the pair potential between atoms $i$ and $j$, $r$ is the distance between atoms $i$ and $j$, and $q_{i}$,$q_{j}$ are the charges of atoms $i$ and $j$, respectively. The pair parameters $A_{ij}$, $\rho_{ij}$, and $C_{ij}$ are free parameters that are obtained by fitting to material properties as discussed earlier. This potential can be extended to take the polarization effects of the nucleus and electron shell into account by using the shell-core model developed by Dick and Overhauser [21]. In this model, the charged nucleus and electron cloud are treated as a massless, negatively charged shell bound by a spring to a massive, positively charged core. The spring constant for this model is determined by the polarizability of the atoms being modeled.

The second potential form consists of the addition of a Morse potential to the Buckingham potential form in Equation 2.5, which is used to describe the covalent bonding between the anions and cations. This form can be described by:

$$V_{ij}\left(r\right) = \frac{q_i q_j e^2}{4\pi\epsilon_0} + f_0\l... ...\right) \exp{\left(\frac{a_i + a_j - r}{b_i + b_j}\right)} - \frac{C_iC_j}{r^6}$$ (2.6)

where $a_i$,$a_j$,$b_i$,$b_j$,$c_i$, and $c_j$ are the free parameters that are obtained by fitting to material properties as discussed earlier. Table 2.1 lists the fitted parameters for the potentials considered in this study.

Table 2.1: Shell-core parameters for Buckingham form pair potentials in CeO$_2$.

Parameters	Units	Potentials
		Gotte [22]	Minervini [7]	Sayle [23]
O shell charge	e	-6.5667	-2.04	-6.1
O core charge	e	4.5667	0.04	4.1
O spring constant	eV/Å$^2$	1759.8	6.3	419.9
Ce shell charge	e	4.6475	4.2	7.7
Ce core charge	e	-0.6475	-0.2	-3.7
Ce spring constant	eV/Å$^2$	43.451	177.84	291.75
O - O interactions
A	eV	9533.421	9547.96	22764.3
$\rho$	Å	0.234	0.2192	0.149
C	Å$^6$	224.88	32	43.83
O - Ce interactions
A	eV	755.1311	1809.68	1986.83
$\rho$	Å	0.429	0.3547	0.35107
C	Å$^6$	0	20.4	20.4
O - La interactions
A	eV		2088.79
$\rho$	Å		0.346
C	Å$^6$		23.25

Previous work has been done using both modeling and experimental techniques to investigate the clustering of oxygen vacancies around dopant ions in ceria doped with trivalent ions (in this system, the lanthanum trapping effect). Wang et al. showed that in the dilute range, charged dimers $(M_{Ce}':V_O^{\cdot\cdot})$ form [24]. Gerhardt-Anderson and Nowick extended this work on other $CeO_2:M_2O_3$, and suggested that conductivity and dimer binding energy vary inversely with dopant cation radius [25]. Kilner and Brook found that due to the size mismatch between the host and dopant cations, elastic strain energy makes a large contribution to the binding energy of the dimer [26]. Subsequent theoretical studies also showed the importance of defect clustering in the determination of free charge carrier concentration in fluorite oxides [27][28][29].

Previous work has also been carried out to validate this type of local configuration dependent kinetic Monte Carlo technique. Murray et al. carried out one of the earliest investigations modeling oxygen vacancy conductivity in yttria-doped cerium oxide [30]. They showed that the oxygen vacancy barrier is sensitive to the local dopant environment and that a first nearest neighbor approximation could generate ionic conductivity results that are somewhat consistent with experimental results. Pornprasertsuk et al. used a similar KMC approach by calculating the binding energies of the oxygen vacancies and dopant ions in Yttria-stabilized zirconia [31]. They showed that the association energy of oxygen vacancies and dopant ions was big enough compared to the oxygen vacancy migration barrier to make the vacancy migration energy depend on the local dopant environment. These studies helped to validate the use of local configuration-dependent migration energies in kinetic Monte Carlo simulations, but these codes were written specifically for the systems in question. The goal of this study was to create a generalized code that could be used to simulate arbitrary systems with arbitrarily complex local configurations, so long as the migration energies for these configurations are known.