Starting from pure ceria, let's consider first the formation of anti-Frenkel (AF) defects: an oxygen ion jumps to an interstitial site, leaving an V_O at the O site (Mamontov and Egami, 2000). Following the Kröger-Vink notation (Kröger, 1977), the defect equation is:

AF defects are supposed to be involved in the oxygen storage capacity of these materials, especially in nanocrystalline form (Mamontov et al., 2000).

V_Os can be also formed by the interaction of CeO_2−δ with the atmosphere, following the equilibrium:

Oxygen deficiency δ values vs. temperature and oxygen partial pressure pO₂ from different authors were reviewed by Mogensen et al. (2000).

For each V_O, two electrons are injected in the conduction band (Tuller and Nowick, 1979). The pertinent equilibrium constant is:

where pO₂ is the oxygen partial pressure. Electrons produced in equilibrium (2) form adiabatic small polarons and electronic conduction is achieved via an activated diffusive hopping mechanism involving Ce⁴⁺ to Ce³⁺ reduction (Tuller and Nowick, 1977; Chiang et al., 1996; Oliva et al., 2004; Farra et al., 2013). Electronic conductivity σ_e can be expressed as:

In Equation (4), n_e and q are the electron concentration in the conduction band and electric charge, respectively, μ₀ is the T-independent part of the electron mobility term, E_H is the hopping activation energy of the polaron and k is the Boltzmann's constant.

Considering the electro-neutrality condition 2[VO••]=[e′], electron concentration in the conduction band and oxygen partial pressure are related by:

Doubly ionized V_Os are dominant at low δ, while defects interactions and single ionized V_Os become important at larger δ (Tuller and Nowick, 1979), bringing to different equilibria (Mogensen et al., 2000). In any case, an oxygen partial pressure dependent conductivity is an efficient way to detect electronic contribution to charge transport in non-stoichiometric oxides (Scavini et al., 1994).

Upon doping ceria with trivalent RE ions, RE³⁺ substitute Ce⁴⁺ and V_Os are introduced in the structure for charge compensation. Doping occurs during the synthesis process and the reaction can be described using the following equation:

Equation (6) assumes trivalent cations only. If not, V_O concentration differs from x/2.

Ceria doped materials display high oxygen mobility via a vacancy diffusion mechanism and, as a consequence, high ionic conductivity σ_i. σ_i can be expressed as Tuller and Nowick (1975):

where ni=[VO••], is proportional to the probability to find a vacant site around the jumping oxygen ion, E_i and ΔS_i are the activation energy (Enthalpy) and Entropy for oxygen diffusion, respectively, q is the charge (=2e), k is the Boltzmann constant, υ is a frequency factor and a/2 (i.e., half of the cell parameter) is the jump distance for an V_O, along the <100> crystallographic direction (Mogensen et al., 2000; Koettgen et al., 2018). According to equilibrium (6), RE doping introduces additional V_Os; also, the presence of the charged RE′_Ce defects changes the electro-neutrality condition into 2[VO••]=[RECe′]+[e′] and, as a consequence, the dependence of n_e vs. pO₂ described in Equation (5) is turned into ne∝pO2-14 for [RECe′]≫[e′] (Tuller and Nowick, 1975). Finally, for large x values, equilibrium (2) is pushed toward its left side and the concentration of conducting electrons is negligible but at very low oxygen partial pressures (see below).

Electrochemical Impedance Spectroscopy (EIS) is usually adopted to measure conductivity in materials in which the ionic conduction is prevalent on the electronic one (Sacco, 2017). A small sinusoidal voltage V = V₀sin(ωt) is applied and the response current I = I₀sin(ωt+α) is measured at the same frequency. As a consequence of the (possible) phase shift α, the impedance Z calculated through Ohm's law is a complex number (Z = Z′ + iZ″ = V/I). A wide range of frequencies ω are sampled and Z data are typically plotted using the Nyquist representation. Data are then fitted against an equivalent circuit that is a combination of resistive R and capacitive C terms considering the transport across the bulk (R_b and C_b) and the grain boundary (R_gb and C_gb) of a polycrystal. Further elements are usually added for the electrode-electrolyte surface contribution and/or electronic transport in mixed ionic/electronic conductors.

In Figure 2 conductivity data (from Eguchi et al., 1992) on (Y, Gd, Sm)-doped CeO₂ are plotted against T. Activation energies for conduction can be calculated fitting the Arrhenius plots (logσT vs T⁻¹) of Figure 2. The plots of the insets of Figure 2 (logσ vs. pO₂) assess the electronic and ionic contributions to conductivity: at high pO₂, σ is constant and fully ionic; conversely σ increases lowering pO₂ by reason of additional electronic contribution. For application as a solid electrolyte in SOFC/SOEC it is important to assess the zone of x, T and pO₂ where the equilibrium with the atmosphere brings to the formation of additional V_Os according to equilibrium (2) because electronic conduction provokes undesired shortcuts between the electrodes.

A huge number of experimental and review papers focus on the transport properties of RE-doped ceria solid solutions as a function of temperature, oxygen partial pressure, composition and microstructure. Limiting to reviews (Jacobson, 2010) and (Goodenough, 2003), discussed the suitability of Ce_1−xRE_xO_2−x/2 materials for SOFC applications while Inaba and Tagawa (1996) and Kilner (2008) and, more recently, Koettgen et al. (2018) reviewed their electrical conductivity.

Although data from different groups are highly scattered, as pointed out by Koettgen et al. (2018), some general trends emerge. As suggested by Equations (6, 7), isothermal ionic conductivity increases with x. However, this applies to light doping only. σ_i decreases above a critical concentration x_c depending on RE nature and temperature, as observed for Gd (Faber et al., 1989; Zhang et al., 2002, 2004; Zha et al., 2003), Sm (Zhan et al., 2001; Jung et al., 2002; Zha et al., 2003), Y (Wang et al., 1981; Faber et al., 1989; Zhang et al., 2004; Sato et al., 2009) and other dopants like La, Nd, Yb (Faber et al., 1989; Dikmen et al., 1999). The critical concentrations x_c reported in the literature are scattered and usually occur in the range x_c = 0.06-0.2, depending on RE and temperature considered (Koettgen et al., 2018). For example, maxima for Gd were observed by different authors at x_c~0.10 (Steele, 2000), 0.15 (Zha et al., 2003), and 0.20 (Zhang et al., 2002).

As an example, σ_i(x) vs. T registered for Ce_1−xY_xO_2−x/2 are reported in Figure 3A (data from Faber et al., 1989). Isothermal conductivity maxima appear around x ≈ 0.10, at which only 2.5% of the oxygen sites are vacant. The measured activation energies also depend on x and show minima around the same x values (Kilner, 2008). In Figure 3B the E_i values as a function of x are displayed (data from Faber et al., 1989).

E_i and σ_i also change by moving across the lanthanide series: the highest σ_i is observed for Gd and Sm (Eguchi et al., 1992; Chen et al., 1997; Omar et al., 2006; Zajac and Molenda, 2008), while σ_i reduces for lighter and heavier RE elements. The inset of Figure 3B shows E_i and σ_i at T = 400 K and x ≈ 0.10 as a function of the RE ionic radius (data from Faber et al., 1989).

To rationalize the bell curves of Ei and σ_i vs. the ionic radii, several authors evidenced the role of ionic radii (ir) mismatch between Ce⁺⁴ and RE⁺³. It should be noted that, for plotting purposes, in Figure 3B we used the ionic radii of the 8-fold-coordinated Re ions (Shannon, 1976). However, the correct coordination to calculate the “mismatch” is matter of debate and it will be discussed in section Lattice Parameters.

Following the “mismatch” idea, too small and too big RE³⁺ cations cause large structural distortions depressing conductivity. Codoping with a suitable combination of small and big cations has been proposed to minimize elastic strain and improve conductivity (Yoshida et al., 2001; Wang et al., 2004; Omar et al., 2006, 2007; Kilner, 2008; Zajac and Molenda, 2008). For example, it has been observed that Sm and Nd-co-doped ceria is more conducting than single Gd-doped ceria at 550°C (Omar et al., 2008); on the other hand, Ce_0.90Lu_xNd_yO_1.95 is less conductive than Ce_0.90RE_0.10O_1.95 with RE = Y, Sm, Gd (Omar et al., 2006).

Conductivity measurements display complex behaviors moving along the x and r_i coordinates. Actually, ionic conductivity is a macroscopic quantity that sums up (and average) a plethora of elementary migration processes each one involving a jump of one oxygen ion from a lattice site to an empty vacant site around. Thus, the local structure rather than the average one does influence these processes. As pointed out by several authors, part of the V_Os could be excluded from the diffusion process (at least within a certain temperature range) because of RE-V_Os defect clustering. V_Os involved in defect clusters should be thermodynamically more stable than the remaining ones thus hindering or reducing the jump frequency toward them. On the other hand, the activation energy for diffusion could depend on the local structure around the jumping coordinate and also include a term for the interaction of the vacancy with other point defects (Kilner, 2008). Koettgen and coworkers identify two phenomena: “trapping” and blocking. In the former the RE distribution affects the initial and the final oxygen sites energies differently; migration energy is therefore different for forward and backward jumps. In the latter, forward and backward energies are equal but dopants affect the transition state energy (Koettgen et al., 2018).

To apply these ideas to real ceria-based phases it is fundamental to map the RE and V_Os distribution as a function of RE nature, concentration and temperature. This implies knowledge of the structure at different length scales using different experimental techniques spanning from diffraction to spectroscopy, microscopy and magnetic resonance. This is the core of the present review and some findings are presented in the next section. Then, coupling theoretical calculations to experimental findings should allow estimating also the energies at work on defect clusters formation and on oxygen migration paths (see section Atomistic Modeling Methods of Doped Ceria for some details).

Cerium oxide exhibits fluorite structure: cubic, space group Fm-3m, with Ce and O in special positions, namely Ce in 4a (0, 0, 0) and O in 8c (¼, ¼, ¼). Ceria is used as reference standard material by NIST. Fluorite is stable over a wide range of T and oxygen non stoichiometry δ. No phase transformation was found down to 2 K by neutron powder diffraction (Coduri et al., 2012a), while high temperature promotes the splitting of O into two different sites consistent with a disordering along <111> direction (Yashima et al., 2006).

When CeO₂ is reduced to CeO_~1.7−1.8, a disordered, non-stoichiometric, fluorite-related phase (α) forms (Bevan and Kordis, 1964; Trovarelli, 2002). By further increasing δ, a number of fluorite related superstructures were observed, lowering the cell symmetry owing to vacancy orderings along the fluorite <111> direction, consistently with high temperature studies (Yashima et al., 2006). A comprehensive description of the superstructures in CeO_2−δ, down to CeO_1.66, is given by a single crystal neutron diffraction study (Kuemmerle and Heger, 1999).

Though still fluorite, the atomic scale structure can be different. An example for CeO₂ is given in Mamontov and Egami (2000). Using Pair Distribution Function (PDF), they probed Frenkel-type defects in octahedral sites. These defects, which can be removed by high-temperature treatment, are related to the oxygen storage capacity of ceria. Hence, the structure is fluorite, but the application is driven by defects that alter, locally, the fluorite structure.

More in general, fluorites look simple, but they are not. Fluorites accommodate high concentration of lattice defects, especially in terms of V_Os (Kim, 1989; Malavasi et al., 2010). Their wide temperature and stoichiometric stability can hinder the real local atomic arrangement, which in complex materials is often the one at the basis of the application (Egami and Billinge, 2003) This makes fluorite more interesting, but more challenging, to control and investigate.

One of the first local scale investigations came with a pioneering neutron powder diffraction study (Anderson and Cox, 1983), where different dopants-vacancy clusters were tested to model diffuse scattering. Though limited by instrumental setup, they demonstrated the existence of local ordering breaking the symmetry of fluorite, proposing the formation of clusters of dopant and V_Os. Yet the structural information was rather qualitative. In the last decades a number of local structure investigations were carried out. This section will provide a survey, technique by technique, for unveiling the defect structure in doped ceria.

Nowadays different notations are used to describe the compositional region corresponding to the transition from fluorite to C-type. In Coduri et al. (2013a, 2018), Checchia et al. (2015), and Scavini et al. (2015) the authors defined as C^* the compositional region (x~0.3–0.5 depending on composition) characterized by the percolation of C-type domains. From a powder diffraction point of view, the C^* region is characterized by (i) non negligible concentration of APBs, i.e., anomalous broadening of superstructure peaks; (ii) a different linear dependence of atomic coordinate x(M2) with doping than the x(M2) trend in the “mature” C-type phase. This means that the structure is actually long range C-type, resulting from the percolation of C-type nanodomains in a fluorite matrix while large fluorite regions still exist.

Artini et al. (2015, 2016, 2017) described intermediate compositions for RE = Gd, Sm, Lu as a hybrid (H) phase between fluorite and C-type. The structure is generally C-type, but the coexistence of Raman modes typical of fluorite and C-type leads tothe picture of a hybrid structure, i.e., a sort biphasic system occurring at the local, or intermediate scale. This is similar to the C^* formalism, even though the H structure can be found in long range fluorite samples—the case of uncorrelated nanodomains in a long range fluorite structure does not apply to C^*—and extend to a larger region in the C-type compositional range. The existence of H was found to depend on the size mismatch between dopant and host. If the size mismatch is too big, a single hybrid (C-type) phase is not stable, i.e., it cannot accommodate for both Ce and RE, and a long range phase separation into F and C-type is observed. Producing samples following the same protocol, Artini et al. (2016) observed long range F-C separation only for Tm, Yb, and Lu. Phase separations are observed for other dopants as well (Bevan and Kordis, 1964; Shuk et al., 2000). The reaction procedure does play a role in that, as either single fluorite or biphasic systems were obtained on the same material, respectively, after oxidizing and reducing annealings (Małecka et al., 2009).

Doped ceria can be successfully produced through a number of different synthetic routes. Most of early studies used solid state synthesis, involving reaction of CeO₂ and RE₂O₃ mixed powders at high temperature with intermediate regrinding steps. Still, the observed local ordering with dopant oxide structure can either be an intrinsic characteristic or a consequence of an incomplete reaction.

The interest for surface effects called for the development of different synthetic routes to tune the morphology of the products. Manifold wet procedures were proposed (Van Herle et al., 1998; Reddy et al., 2009; Rezaei et al., 2009; Wang et al., 2010). These lead to nanoparticles that can be sintered afterwards to increase grain size and density. Moreover, wet methods are considered to lead to more homogeneous materials (Horlait et al., 2011).

The role of the synthesis is generally studied with respect to transport properties. Its relationship with defect structures or structure in general received less attention, with the exception of defects in nanometric samples, where surface effects become not negligible.

Tsunekawa et al. (2000) attributed the inverse relation of lattice expansion of CeO₂ and crystal size to the stabilization of Ce⁺³ ions. In doped samples Ce⁺³ is often found to segregate at the surface of nanoparticles. Lee et al. (2014) observed that Y-doping for x > 0.09 induced ferromagnetic ordering of surface Ce⁺³ clustered together with Y⁺³ ions. Hence, Ce⁺³ ions stabilized at the surface replace some Y⁺³ ions in the M⁺³-V_O clusters already observed in bulk. Similarly, Sm³⁺-V_O-Ce³⁺ complexes were observed for x > 0.07 on the surface of Sm-doped nanoparticles (Chen et al., 2014), suggesting a core-shell defect structure, already proposed in Malecka et al. (2008); Małecka et al. (2009) in terms of surface RE segregation.

Recently, 3D electron microscopy revealed the 3D surface of La-doped ceria with subnanometric resolution (Collins et al., 2017). Additional V_Os were seen within the first 1.5 nm of the particle surface together with associated changes in Ce(4f) hybridization as well as surface enrichment in La. Acharya et al. (2014) observed through Raman and EXAFS that in nanoparticles Gd induces more intrinsic V_Os than Sm. For a fixed doping amount, reducing the particle size damped the Raman signal at 370 cm⁻¹, which is the fingerprint of C-type ordering (6-fold coordination). Accordingly, the C-type domains observed in bulk samples disappear, or reduce in intensity, when decreasing particle size. No PDF peak of C-ordering is observed for x = 0.313 Gd-doped samples (Coduri et al., 2017) when particle size is below 10–15 nm. Similarly, lower sintering temperatures in Y-doped samples (x = 0.10 and 0.25) led to smaller nanodomains (Ou et al., 2006a). Eventually, Sen et al. (2008) observed that nanometric particle size doubles the population of Y(VIII), increasing the probability of V_Os to be NN as Ce, which can thus be reduced to +3. This can enhance the migration of V_Os, which are less strongly bound to Ce⁺³ than to Y⁺³ (Mogensen et al., 2000).

In conclusion, reducing particle size to the nanometric scale increases the fraction of Ce⁺³, mostly located at the surface and strongly associated with V_Os. Ce⁺³ acts as a dopant even larger than Nd⁺³ and La⁺³, which are known to shift the long range F to C phase transformation to higher dopant concentrations. Ce³⁺ thus stabilizes fluorite. XRD, PDF and Raman provided evidence in this direction.

The above investigations were performed at low temperature, which is the best condition for structural analysis. Yet, this is far from real operating conditions, i.e., high temperature and controlled (oxidizing and reducing) atmosphere. Only a few investigations reported in the literature were not performed at ambient condition, and generally they are carried out under air. This is due firstly to instrumental difficulties, such as the case of electron microscopy. Secondly, increasing atomic vibrations broaden the signals of local probes, undermining the significance of local scale investigations. High temperature does not affect diffraction studies, other than making intensity decaying faster with momentum transfer. Since doped ceria retains its fluorite structure on a wide temperature range, no structural discontinuity has ever been reported. Yashima used high temperature powder diffraction to probe static disorder by analyzing APDs (Yashima and Takizawa, 2010). No deviation from linearity with temperature of lattice parameters and APDs using neutron diffraction on x(La) = 0.25 up to 800°C (Coduri et al., 2013b) was observed, nor it was for different compositions of Sm (Artini et al., 2018) and Gd (Artini et al., 2014). This evidences once more the need for a local probe. To the authors' knowledge, only Wang et al. (2006) used EXAFS on x(Y) = 0.10 doped ceria at high temperature (600°C). Although Y- and Ce-environments appeared more homogeneous at high temperature, they claimed that structural features become ambiguous because of excessive signal broadening.

Raman suffers as well from high temperature effects. Still, the Raman band assigned to dopant-V_O clusters was found to vanish at 450°C for x(Gd) = 0.15, consistently with the dissociation of V_Os from clusters involving dopant ions. These clusters are preserved for higher dopant amount. A similar behavior was found for other dopants (Shirbhate et al., 2016). PDF showed that local C-type orderings are retained up to 750°C, even though to a lesser extent than at RT. Further distortions, compatible with Ce reduction, were observed under reducing atmosphere (Coduri et al., 2013b).

HRTEM was used to compare materials before and after being subject to operating conditions, revealing that reducing atmosphere induces further V_Os ordering consequent to Ce reduction (Li et al., 2013). Increasing sintering temperature promotes a melting, at least a partial one, of the clusters/domains rich in V_Os. Surprisingly, among the vast set of reports in the literature, very little effort has been dedicated so far to move defect structures investigations toward real operating conditions.

Density Functional Theory (DFT) is currently the computational workhorse for the first principles modeling of doped ceria, surfaces and catalytic reactivity. However, solving the Schrödinger equation with standard DFT methods presents one problem, that is, the spurious self-interaction between electrons. The substitution of Ce⁴⁺ with a trivalent ion leaves an unpaired hole in the system. It turns out that in many oxides, such as MgO, TiO₂, SiO₂, CeO₂, both local (LDA) and semilocal (GGA: BLYP, PW91, PBE, PBEsol) functionals yield a hole wavefunction which is delocalized on several ionic sites. Indeed, reduced CeO₂ is predicted to be metallic by GGA functionals, in stark contrast with experiments. Moreover, EPR experiments, through the hyperfine couplings, reveal that the spin-density is very much localized and UPS experiments show that new states, associated with Ce 4f orbitals, appear at ~1.2 eV above the top of the valence band. This drawback carries several additional consequences and leads to wrong predictions about lattice spacing, vacancy formation and migration energies, surface energies, catalytic activity (Pacchioni, 2008).

Two practical approaches emerged to partially correct for the self-interaction issue: DFT+Hubbard U (Dudarev et al., 1998; Cococcioni and de Gironcoli, 2005) and hybrid functionals (Dovesi et al., 2005, 2014). Both methods are rooted in the principle that Hartree-Fock (HF) is almost free of self-interaction. In fact, DFT+U can be viewed as a local-HF correction, acting on the subspace of localized d and f electrons. Similarly, hybrid functionals evaluate the exchange potential by mixing pure orbital-dependent HF and density-based exchange. Hybrid functionals like B3LYP and PBE0 retain the long range Coulomb interaction while the HSE (Heyd et al., 2003) uses a screened Coulomb interaction, which appears to be more appropriate to deal with periodic solids. However, both methods depend on empirical parameters: the value of U and the fraction of HF, respectively, which makes them not fully ab-initio.

There are procedures to determine the U-values (Cococcioni and de Gironcoli, 2005) and the fraction of HF from the dielectric constant (Skone et al., 2014) but they are not widespread in calculations. Rather, it is common practice to adjust these parameters to reproduce some experimental quantities, such as lattice spacing or vacancy formation energy (Fabris et al., 2005a; Pacchioni, 2008). Lu and Liu (2014), presented a rationalization of the Hubbard U parameter for Ce oxides, calculated within different approaches (linear response, constrained RPA). They found that the U parameter ranges from ~4.3 to 6.7 eV and depends strongly on the Ce-O coordination number and bond length. A similar conclusion was reported in Loschen et al. (2007). Recently it has been found that adding a Hubbard U term on oxygen 2p orbital (despite the fact that is not as localized like the 4f orbitals), does improve the electronic structure and the energetics of defective ceria (Yeriskin and Nolan, 2010; Plata et al., 2012). In particular, Figure 9 shows the spin density around a V_O. By employing DFT+U spin density remains strictly localized on two Ce³⁺ ions. Without U corrections, it would spread over many more sites.

Hybrid functionals constitute a valid alternative to DFT+U, at the price of a higher computational cost, especially in the plane-waves. Among all hybrid functionals, the screened-exchange HSE06 functional is the most used to describe bulk CeO₂ and Ce₂O₃ (Hay et al., 2006; Da Silva et al., 2007; Ganduglia-Pirovano et al., 2007; Du et al., 2018), oxygen vacancies at surfaces (Nolan, 2010, 2011; Han et al., 2016), polarons (Sun et al., 2017) and dopants (Shi et al., 2016). Other non-screened functionals, such as B3LYP and B3PW91 were employed successfully to calculate the optical properties of bulk CeO₂ (El Khalifi et al., 2016). In general, hybrid functionals provide a better description of both Ce³⁺ and Ce⁴⁺ ions, larger band gaps, larger vacancy formation energies and more localized hole wavefunctions, with respect to pure DFT and DFT+U. Finally, there have been only few but promising investigations on bulk CeO₂ employing meta-GGA functionals (Tran et al., 2006).

Most calculations reported in literature employ the plane-wave pseudopotential method (Kresse and Joubert, 1999; Giannozzi et al., 2009) and all-electron local-basis calculations (Dovesi et al.) are less frequent due to the large number of electrons of the lanthanides. Pseudopotentials instead treat explicitly only valence electrons and they often discard the 4f electrons (putting them in the core) in order to reduce the computational cost and complexity (i.e., discarding magnetism; Dholabhai et al., 2010).

To describe different doping concentrations, one typically builds periodic supercells consisting of 96 sites (2 × 2 × 2 of the conventional cubic CeO₂). In the presence of multiple vacancies, one must pay attention to spin multiplicity and to the fact that DFT+U can present multiple local minima depending on the starting guess of the orbital occupations. It was found that the energy difference between parallel and anti-parallel spins is very small (hence, a ferromagnetic solution can be imposed, Murgida et al., 2014) and that the full randomization of the starting guess can confidently provide the true electronic ground state.

To determine the optimal dopants and vacancies configuration is a formidable task even for a 96-atoms supercell, where the number of combinations can be very large. To solve this problem, genetic optimization algorithms proved to be efficient (Jung et al., 2018). However, such periodic supercells are too small to calculate accurately the thermodynamics, vacancy ordering and clustering effects. To overcome the finite-size effects, several authors employed cluster expansion techniques (van de Walle et al., 2002). In short, the idea is to map the local atomic configurations of dopants and vacancies on a Heisenberg lattice model. The “J” parameters of Heisenberg Hamiltionian are obtained by total energy differences between different configurations. Next, Monte Carlo methods are used to include temperature effects and to provide useful insights on vacancy order/disorder and on the microscopic structure at the nanometer scale (Gopal and Van De Walle, 2010; Murgida et al., 2014; Žguns et al., 2017, 2018).

Despite the tremendous progress of first-principles techniques and code in the past few years, both DFT+U and hybrid functionals still need to be improved in their formulation or in their parametrization. This would reduce the dependence on experimental inputs and provide more accurate, ab-initio predictions. The development of DFT-based tight-binding methods (i.e., DFTB) allows to simulate larger supercells that can represent closely the local structure, the clustering and ordering of defects, minimizing the spurious interaction between their periodic replicas (Kullgren et al., 2017).

As shown above, the combination of DFT methods and cluster expansion techniques is extremely powerful but also computationally involved. In addition to first-principle methods, doped ceria have also been studied by Molecular Dynamics (MD) and Monte Carlo (MC) methods, using empirical potentials. The most common potentials are pairwise Born and Mayer (1932) rigid-ion potentials. In particular Vives and Meunier (2015) studied the defect structure and defect migration paths with six different parametrizations of the BM potentials. By comparing to experimental data they found that only two parameter sets can be used confidently to study Gd-doped CeO₂. For a recent review of atomistic simulations of oxide interfaces, surfaces and nanoparticles (see Sayle and Sayle, 2007).

In this section we selected from the vast literature recent results on ab-initio modeling of doped ceria. Vanpoucke et al. (2014) studied extensively the aliovalent doping of CeO₂ with lanthanide, transition metal, and alkali atoms with DFT+U and determined the parameters of the Vegard's law, the defect formation energy, bulk modulii, thermal expansion parameters as a function of dopant concentration. Apostolov et al. (2018) and Jung et al. (2018) report phonon calculations and Raman shifts of the F_2g mode of pure CeO₂, as a function of doping. In particular, doping with lanthanide +3 ions (whose Shannon radius is larger than that of Ce⁴⁺) causes a decrease in the frequency of the F_2g Raman mode and considerable local structure relaxations. Dholabhai et al. (2010) and Ahn et al. (2012) and their coworkers studied the oxygen vacancy migration barriers in Pr-doped ceria. Surfaces of pure and doped ceria have been studied in Yeriskin and Nolan (2010) and Farra et al. (2013) where they have found the appearance of Ce³⁺ and one compensated oxygen vacancy. The electronic structure of planar and stepped surface was studied by Fabris et al. (2005b) and Esch et al. (2005).

The ionic conductivity of doped and co-doped ceria has been studied by MD with empirical potential by Burbano et al. (2014) using a sophisticated polarizable ion model fitted to ab-initio results. Their main finding is that oxygen conductivity in co-doped ceria can be well approximated by the weighted average of the conductivity of the single-doped parent compounds. According to their analysis the oxygen conductivity is determined by the local lattice strain generated by the single defect rather than by synergistic effects between different dopants. This confirms EXAFS investigations (Yoshida et al., 2001).

Purton et al. (2017) simulated the ionic conductivity of Ca- and Gd-doped CeO₂ using large simulation cells. They employed a hybrid-MC technique consisting of a short (1 fs) MD run followed by a MC exchange move between a Ce and a dopant ion. The purpose of the MC move is to provide an escape path from local-minima configurations in order to approach the equilibrium in a faster way. They found that for x(Gd)>0.2, Gd-rich domains are formed. They also calculated the oxygen conductivity for randomized Gd positions and Gd-nanodomains. They found, in accordance with percolation theory, that Gd-rich domains limit oxygen mobility. They also studied the segregation of the dopants at grain boundaries of CeO₂.

Koettgen et al. (2018) reported an ab-initio study of the oxygen conductivity, employing DFT and Kinetic Monte Carlo (KMC). As already mentioned above, they considered only jumps in the <100> direction by half of the unit cell, and only three local configurations (Ce-Ce, Ce-RE, RE-RE). Despite the small database of moves, their KMC results are in very good agreement with experiments. Contrary to observations of Burbano et al. (2014), they found that the ionic conductivity is influenced by trapping, blocking and vacancy–vacancy synergistic interactions: blocking limits the dopant fraction at the ionic conductivity maximum while trapping limits the maximum ionic conductivity. They also found a non-linear Arrhenius behavior of the conductivity, with a reduced activation energy at high temperature, which they ascribe to the “association” between an oxygen ion and a rare-earth dopant.

The thermodynamics of defective ceria was studied by cluster expansion methods by Gopal and Van De Walle (2010); Žguns et al. (2017, 2018). In particular, Gopal et al. addressed the thermodynamics of intrinsic oxygen vacancies. In addition to configurational entropy they also included the lattice vibrational contribution to the free energy. Their lattice Monte Carlo simulations showed that vacancies have a tendency to cluster and to order along preferred directions. DFT calculations showed for instance that for small radius dopants (La to Nd), the oxygen vacancy tend to occupy the nearest-neighbor (NN) position, whereas large radius dopants (Sm to Er) the oxygen vacancy occupies the next-nearest-neighbor (NNN) site, in order to minimize the elastic energy (Andersson et al., 2006; Gupta et al., 2010). This is good agreement with EXAFS experiments (see Section X-Ray Absorption Spectroscopy). Using empirical potentials, Hayashi et al. (2000) showed that by doping with La, Gd, and Y, the local structure of the dopant-vacancy complex is characterized by oxygen relaxation along the [100] direction. As a consequence the Ce-O distance decreases even if the average volume of the fluorite cell is increased.

The ordering of Gd dopants was studied extensively by Žguns and coworkers in two recent papers (Žguns et al., 2017, 2018). Their Monte Carlo simulations of both thermal equilibrium and rapid quenching showed a critical temperature just below 1,000 K. Their simulation showed that below T_c Gd tends to phase-separate from CeO₂, forming nano-domains of C-type Gd₂O₃ (see Figure 8 of Žguns et al., 2017). At high temperature the Gd ions take a completely random distribution but oxygen vacancies tend to cluster in the coordination shell of the dopant. They further studied the entire CeO₂-Gd₂O₃ phase diagram, in the full concentration range, as a function of temperature (Žguns et al., 2018).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The handling editor declared a past co-authorship with several of the authors (MC and MS).