Molecular dynamics investigation of the structural and energetic properties of CeO₂–MO_x (M = Gd, La, Ce, Zr) nanoparticles

Abstract

CeO₂-based materials have played a critical role in catalysis, where the substrate particles have reduced in size year by year due to experimental achievements in synthesis control. Thus, there is increasing interest to improve our atomistic understanding of the structural and energetic properties of mixed CeO₂-based nanoparticles of 1 nm to 5 nm. Here, we employed classical molecular dynamics (MD) simulations to study the following solid solutions, CeO₂–Gd₂O₃, CeO₂–La₂O₃, CeO₂–Ce₂O₃, and CeO₂–ZrO₂, using 5 compositions (0, 25, 50, 75, and 100%). The amorphization and re-cystallization process via MD simulations was employed to generate the nanoparticles, which were characterized by several analyses. We found that even in small CeO₂ concentrations, the systems maintain the cubic fluorite structure and the truncated octahedron shape found in pure ceria, evidencing the strong influence of Ce⁴⁺ on the nanoparticle morphology due to the higher phase transition temperature of the CeO₂ compound. On the other hand, the addition of the 3+ species leads to the spontaneous appearance of higher concentrations of solvated cations and vacancies near to the surface of the CeO₂-based solutions. Beyond that, the 3+ species also influence the electrostatic potential in the nanoparticle surface, and hence controlling the 3+–4+ ratio may be an interesting approach to control the nanoparticle physicochemical properties for catalytic purposes.

---

1 Introduction

Cerium-based oxides have played an important role in several industrial applications,^1,2 which spread from catalytic^3,4 to biomedical applications.^5,6 Recently, G. Lu et al. highlighted their importance in electrocatalysis, namely, in the hydrogen and oxygen evolution reactions.⁷ Those applications have been possible due to the high figures of merit of cerium oxides (CeO_x, ), in particular, their oxygen storage capacity, which is interconnected with the formation of oxygen vacancies and changes in the oxidation state of Ce cations, i.e., from Ce⁴⁺ (CeO₂) to Ce³⁺ (Ce₂O₃).^8,9

CeO₂ (ceria) has a fluorite crystal structure (space group Fm3̄m), in which each cerium atom is bound to 8 oxygen atoms, while each oxygen is bound to 4 Ce atoms,¹⁰ and its reduced form (Ce₂O₃) crystallizes in the hexagonal La₂O₃-type structure.¹¹ The CeO₂(111) surface has the lowest surface energy, and hence it is the most abundant surface on ceria particles, while the (100) and (110) surfaces are less stable.¹² Those less stable surfaces have been found in large truncated octahedron-like particles with diameters from nano up to micro size.^13,14 Several studies have indicated a direct correlation between the morphology of cerium oxide particles and their applications, especially in the field of catalysis, in which the surface termination and defects play a crucial role in the presence of oxygen vacancy formation, oxygen diffusion, low-coordination sites, oxidation state of the Ce cations, etc.^9,12

Although ceria has high figures of merit, its applications in particular fields require higher stability, easiness to create oxygen vacancies, higher oxygen diffusion, etc. In the literature, we can see that some recent applications of ceria are built taking advantage of the tunability and engineering of those properties in nanoparticles. For example, it is already known that ceria can be used to build hydrophobic nanoparticles for coating, although powder or bulk ceria is slightly hydrophilic.¹⁵ The CeO₂(111) surface is more hydrophobic and the edges are the main factor in the surface–water interaction,¹³ so this application is only possible because of the particle edge abundance. However, our atomistic understanding of the CeO₂ surface–edge relation is far from satisfactory when compared with surface–edge effects in metallic nanoparticles.^16–19 Additional to that, Tarnuzzer et al.²⁰ proposed the use of ceria nanoparticles with diameters from 2 up to 5 nm to repair cell damage induced by radiation, which explores the antioxidant properties of ceria particles with a similar ratio of Ce³⁺ and Ce⁴⁺ cations on the particle surfaces.

Doping or mixing of CeO_x with different oxides is an alternative strategy to tune its physical–chemical properties.^21–25 For instance, it is already established that the formation of a solid solution of ZrO₂ with CeO₂ can enhance the reducibility of Ce⁴⁺ towards Ce³⁺ cations;²² at the same time, solid solutions with different lanthanide atoms can enhance the oxygen conductibility.^21,26 Furthermore, the formation of solid solutions of CeO₂ with Gd₂O₃, La₂O₃, and ZrO₂ has attracted attention for catalysis applications. For example, Gd₂O₃ doped Ce₂O₃ has already been reported as a good ionic conductor for catalytic applications²⁷ and its properties can be improved with co-doping with different rare-earth elements.²⁸ Since these materials are applied typically at high temperatures, the phase diagram is also of interest and some important insights have been obtained in this direction. For example, Andrievskaya et al. reported the La₂O₃–CeO₂ phase equilibrium on large systems for temperatures of 1100 °C to 1500 °C,²⁹ while Zguns et al. theoreticaly predicted similar results for Gd₂O₃–CeO₂ mixtures.³⁰

In theoretical studies, ceria-based mixed oxides can be approached in many different methodologies and scales. To cite a few examples, Balducci et al. studied the surface and reduction energies of CeO₂–ZrO₂ using slab models,³¹ while Bulgakov et al. analysed the solution stability and oxygen diffusion in CeO₂–ZrO₂ bulk models³² and Koleva et al. studied the impact of Zr on the reducibility of small CeO₂ clusters.³³ On the other hand, little is known about the behavior of those mixtures in large solid solution particles with more than 1000 atoms.³⁴ Based on density functional theory (DFT) calculations, Felício-Sousa et al. studied the energetic stability of small Ce_nZr_15−nO₃₀ nanoalloy particles, where they reported a negative excess energy (stable) for all compositions,²⁵ while Mucelini et al. addressed the formation of ZrO₂-like structures on CeO₂(111).²⁴ Although DFT calculations yield an accurate description of the electronic and structural properties, those calculations have an expensive computational cost, and hence the study of temperature effects via molecular dynamics is forbidden for large particles.

A less expensive alternative, based on force field molecular dynamics (MD), was employed by Sayle et al. to study the amorphization and recrystallization of ceria particles, in which it was possible to address successfully CeO₂ particles with diameters up to 8 nm with different shapes, such as nanorods, polyhedral nanoparticles and mesoporous ceria, with a good correlation with the experimental findings.^35–37 Using MD simulations, Sayle and coworkers were also able to correlate the reactivity, morphology, and composition of CeO₂ nanoparticles in different applications such as supercapacitors and the oxi-reduction process.^36,37 Beyond that, the amorphization and recrystallization methodology can be easily extended to mixed particles and other systems with an adequate force-field.^38,39

Therefore, all those experimental^{5,6,10,11,14,15,21,26–29} and theoretical^{8,9,13,22,25,30,31,33,35} studies have provided a great contribution to improve our atom-level understanding of mixed CeO₂–MO_x nanoparticles; however, our atomistic understanding is far from complete. To improve our atomistic understanding, we employed classical MD simulations to investigate the crystallization process of mixed ceria-based nanoparticles, namely, CeO₂–Gd₂O₃, CeO₂–La₂O₃, CeO₂–Ce₂O₃, and CeO₂–ZrO₂. For that, we employed nanoparticles formed by 500 or 864 cations combined with the appropriate number of oxygen anions, which generates particles with sizes from 1 nm to 5 nm.

We found from our analyses that the CeO₂ morphology and structure (fluorite crystal truncated octahedron-like particles) are preserved for all mixed nanoparticles, which can be explained as an equilibrium mechanism between the crystal lattice and the surface energy of each pure material. For CeO₂–M₂O₃, M³⁺ cations and oxygen vacancies (defects) are abundant on the surface region, showing that Coulomb interactions and the ionic radii have a key role as a stabilising mechanism on oxide alloys. In addition, the final alloyed structures have a complex combination of positive and negative electrostatic potentials, modulated by the concentration of the metal atoms, which is advantageous in the catalysis context.

2 Theoretical approach and computational details

We have the aim to obtain a deep atomistic understanding of the crystallization process of finite-size CeO₂–MO_x solid solutions, which is crucial knowledge to enhance the applications of CeO₂-based compounds in catalysis.^40,41 Thus, it requires a description of temperature effects, time-dependent phenomena, point defect formation (e.g., oxygen vacancies), shape (surface) evolution along the crystallization process, etc., where particle size effects play a crucial role. Therefore, it is a challenging problem that requires a strategical balance between an atom-level description and computational cost.⁴² To achieve our goals, we considered two particle sizes, namely, 500 and 864 cations with the appropriate number of oxygen anions, which yields oxide particles with diameters similar to experimental observations,^43,44 which are large enough to show ordered core regions, internal nucleation, and domain walls. Systems with that size are not affordable by ab initio MD calculations, and hence all MD simulations were based on atomistic force field simulations, as described below.

2.1 Molecular dynamics and force field

Our MD simulations are based on the Born-like force field, which is based on the combination of the long-range Coulomb potential and the short-range Buckingham potential.⁴⁶ It has been used for several CeO₂-based studies for a wide range of applications,^47,48 including amorphization/recrystallization simulations of finite-size oxide systems.^45,49 The total charge of each system was set to zero by adjusting the number of oxygen atoms to balance the positive charge (cations), e.g., by building the system directly as the combination of Ce⁴⁺O₂²⁻ and M₂³⁺O₃²⁻. This can be seen as an advantage of the approach proposed here, since the system can be optimized for a given configuration that preserves a specific 3+/4+ ion ratio, which is desirable to understand the oxygen vacancy distribution and its effect on the structural and energetic properties in relation to the system composition.

From several initial tests, we selected the force field parameters developed by Grimes and coworkers⁴⁵ instead of Sayle's original ones.⁴⁹ From our knowledge, the parameters proposed by Sayle were built for different applications, in particular, with a focus on unalloyed oxide systems, and they do not include parameters for all the desired cations in our simulations. On the other hand, Grimes's parameters were derived from mixed oxides of 3+ and 4+ cations, and hence are more suitable for our applications.

The potential term for the total energy of the system is written as follows,

where ε₀ is the vacuum dielectric constant. A, ρ, and C are the Buckingham potential parameters, which are listed in Table 1. It is important to highlight that, because of the stronger nature of the Coulomb repulsion between cations in those systems, it is commonly assumed that the A and C parameters are set as zero and ρ to infinity, as in the original studies reported by Sayle⁴⁹ and Grimes.⁴⁵r_cut1 and r_cut2 are the cutoff radius for the Coulomb interactions and Buckingham potential, which were set as 100.0 and 6.0 Å, respectively. Furthermore, we considered the q_i parameters as fully determined by the atomic species, i.e., q_Ce⁴⁺ = +4.0, q_Ce³⁺ = +3.0, q_La³⁺ = +3.0, q_Ga³⁺ = +3.0, q_Zr⁴⁺ = +4.0 and q_O²⁻ = −2.0.

Table 1 Force field parameters for the short-range Buckingham potential for the cation–anion, anion–anion, and cation–cation interactions obtained from ref. 45 (M = Ce⁴⁺, Ce³⁺, Gd³⁺, La³⁺, Zr⁴⁺). Here, “∞” can be set as any large number in the calculations without a direct impact on the potential values since C = 0 for the M–M interactions

Interaction	A (eV)	ρ (Å)	C (eV Å)
Ce^4+–O^2−	1809.68	0.35470	20.40
Ce^3+–O^2−	2010.18	0.34490	23.11
Gd^3+–O^2−	1885.75	0.33990	20.34
La^3+–O^2−	2088.89	0.34600	23.25
Zr^4+–O^2−	1502.11	0.34770	5.10
O^2−–O^2−	9547.96	0.21920	32.00
M–M′	0.00	∞	0.00

---

All the molecular dynamics simulations were carried out using the LAMMPS package⁵⁰ with a time step of 1 fs and using the standard implementation of the Verlet integration for the canonical (NVT) ensemble. Cubic simulation boxes of 20 nm were set for all simulations, which are approximately five times larger than the largest particles. Periodic boundary conditions were employed for all simulations. Thus, in order to avoid nanoparticle translation and anomalous rotations being counted as contributions to the thermal energy, the particles were re-centered at the origin coordination system at each time step and the angular and linear momentum were set to zero for each 500 time steps. For analysis, the atomic positions and energies were saved at each 500 steps.

2.2 Amorphization and annealing

The framework for the amorphization of ordered structures and their recrystallization was adapted from the procedures reported by Sayle for CeO₂ nanoparticles.³⁵ Based on several MD simulations for smaller and larger particles, our amorphization process includes seven separate steps of NVT simulations:

(1) 200 ps at 4500 K using 2 layers free and 4 layers frozen.

(2) 300 ps at 4500 K using 4 layers free and 2 layers frozen.

(3) 100 ps at 4500 K with all atoms free.

(4) 300 ps at 4500 K with external tension.

(5) Simulation with a temperature ramp from 4500 K to 5500 K for 400 ps.

(6) Simulation at 5500 K for 100 ps.

(7) Simulation over 1000 ps with a temperature ramp from 5500 K to 8500 K, and the amorphous CeO₂ samples were obtained.

To reduce the computational cost, we transformed the amorphous CeO₂ samples into the initial amorphous structures for the recrystallization process for all mixed-oxides systems, namely, CeO₂–Ce₂O₃, CeO₂–Gd₂O₃, CeO₂–La₂O₃, and CeO₂–ZrO₂ with 0, 25, 50, 75, and 100% CeO₂ concentration, by random cation substitution and creation of oxygen vacancies to maintain the total charge of the particles as zero. The recrystallization procedure performs 1 ns lowering of the temperature from 8500 K to 5500 K. Then, we ran subsequent 1 ns production simulations at a constant temperature and 1 ns ramps that decreased the total temperature 100 K at a time. It finishes once the temperature simulation reaches 10⁻³ K (replaces 0 K), i.e., a total simulation time of 112 ns for each system.

2.3 Analyses

Several analyses were performed along the present study, which are defined below, while additional technical details can be found in the ESI.†

2.3.1 Density. To calculate the density of the particles (i.e., the ratio between the mass and volume of the particle), we implemented a dedicated in-house subroutine based on FORTRAN 90. As expected, the nanoparticle mass is computed by the sum of the relative isotopic mass for all atoms in the system. The total volumes of the nanoparticles were estimated by numerical integration of the volume inside the nanoparticle covalent surface via Monte Carlo integration. The randomly generated points for integration were tested using a combination of two criteria: a point is considered inside the nanoparticle volume if (i) its distance to any atom of the system is smaller than the atom covalent radii or (ii) a ray-marching test says that the point “cannot see the simulation box”. Test (ii) is a sophisticated procedure that ensures that the nanoparticle volume estimation includes voids, cavities, and vacancies inside it as part of the nanoparticle volume. For example, the volume calculated for a buckminsterfullerene must include the part inside it that is not filled with atoms, so the final volume of C₆₀ will be closer to the volume of a sphere of the same diameter. Therefore, the points for which there exists at least one line that goes from the point to infinity without crossing any atomic covalent volume are the ones considered outside of the nanoparticle volume. For ray-marching, we employed 1000 directions, while the volume was converged up to 0.1 ų.

2.3.2 X-ray diffraction and radial distribution functions. For each chemical species with different oxidation states, e.g., Ce³⁺ or Ce⁴⁺, we obtained the radial distribution function (RDF) with respect to the geometric center of the nanoparticle for each species and the pair correlation function (g(r)) for each atomic pair. Coordination and tension analysis implemented in OVITO,⁵¹ version 3.1.3, were employed to identify crystalline phases. Furthermore, the LAMMPS X-ray diffraction (XRD) tool⁵² was employed with 100 ps sampling at 300 K and a 1.541838 Å wavelength.

2.3.3 Oxygen vacancy distribution. The search for oxygen vacancies induced by temperature effects in ordered oxides such as CeO₂ nanoparticles can be easily done by the analysis of the cation coordination; however, it is not the case in disordered or partially disordered mixed oxide particles obtained by simulating the annealing process. Thus, in the present work, the location of the oxygen vacancies was obtained using an in-house implementation based on graph theory, which is composed of the following steps: (i) We construct a graph that has the same adjacency matrix (an equivalent form of the first neighbourhood list) for the cations in the nanoparticle. (ii) We search for the full list of K₄ subgraphs in this graph. (iii) For each K₄, we check the existence of an oxygen atom within the volume of the tetrahedron formed by the cations correspondent to the graph vertexes. (iv) If none is found, a vacancy is reported and its position is listed as the center of the tetrahedron of cations corespondent to the K₄ subgraph.

2.3.4 Diffusion barriers via umbrella sampling. To improve our atomistic understanding of the oxygen mobility in pure CeO₂ nanoparticles, we performed potential of mean force calculations through the umbrella sampling (US) method.^53,54 US makes possible the analysis of high energy configurations where the ergodicity of the system is not accessible in a feasible computational time. Initially, we split the trajectory of the oxygen atom into 6 parts, where the first part is the starting position and the last is the final position. Each one of these positions is subjected to the US molecular dynamics independently, where an initial energy minimization stage was performed using the steepest descent algorithm.

After the minimization step, we performed a short NVT equilibration simulation over 10 ps at 2000 K. This temperature is justified, as it is still far from the melting point of the material⁵⁵ and it is enough to enable oxygen diffusion. All the US simulations were performed over 15 ns using the GROMACS package,⁵⁶ version 2018.3, and its in-built tools for data analysis.⁵⁷ Because of the large number of possible atomic arrangements in the CeO₂-based mixtures,⁵⁸ we carried out US simulations focusing only on the pure CeO₂ nanoparticles, which was combined with the oxygen vacancy analysis to improve our atomistic understanding of the oxygen mobility.

2.3.5 Electrostatic potential mapping. To obtain insights into the charge distribution on the nanoparticle surfaces and its dependence on the composition, we computed the electrostatic potential using the smooth particle-mesh Ewald method (PME),⁵⁹ in which the point charges are approximated by Gaussian functions, as implemented in the visual molecular dynamics (VMD) software.⁶⁰ This approach evaluates the long-range electrostatic interactions to obtain the map of electrostatic potential (MEP) by solving numerically the Poisson equation on a numerical grid. For the PME analysis in VMD, the point charge and the van der Waals radius are given for each atom type, where the atoms are treated as a small sphere with the respective charge distributed on the surface. For this analysis, we considered the covalent radius, namely, 2.04 (Ce³⁺), 2.04 (Ce⁴⁺), 1.69 (La³⁺), 1.96 (Gd³⁺), 1.75 (Zr⁴⁺) and 1.52 Å (O²⁻).

3 Results and discussion

Based on the structures obtained by MD calculations, we characterized the nanoparticle morphology, phase transition, stability, radial distribution functions, oxygen vacancy distribution, oxygen mobility, electrostatic potential, etc. For simplicity, our analysis and discussion here will be based only on the largest particles (i.e., the ones with 864 cations). Similar results were found for the smaller particles and are fully reported within the ESI.†

3.1 Nanoparticle morphologies

In the present MD simulations, the liquid to solid phase transition occurs at the temperature at which the shape of the nanoparticle consolidates itself, i.e., there are no significant structural changes for smaller temperatures, except for local atomic rearrangements. Fig. 1a shows the structures of the pure oxides for all systems at 3 selected temperatures, in which the structures obtained below the phase transition are visually similar for all the pure oxides. Thus, for simplicity, we will focus our discussion only on the obtained 0 K final structures.

Fig. 1 Molecular representation of the mixed oxide nanoparticles: (a) Sample geometries for pure nanoparticles at different temperatures during the annealing process. (b) Graphical representation of all the mixed nanoparticle geometries with respect to the M/Ce⁴⁺ ratio (M = Ce³⁺, Gd³⁺, La³⁺ and Zr⁴⁺) at 0 K obtained from the annealing simulations for the structures based on 2592 atoms.

The pure CeO₂ nanoparticle shows a truncated octahedron-like structure, which is expected from previous experimental^61,62 and theoretical studies,^13,35 showing that the combination of the procedure and force-field parameters is working well and presenting proper results. The octahedron-like structure is characterized by the presence of the close-packed (111) surface, which covers most of the nanoparticle facets, which are interconnected by small (100) facets. Those trends can be explained by the highest stability of the oxygen terminated close-packed (111) surface (lowest surface energy) compared with the less stable and open (100) and (110) surfaces. For example, the DFT-PBE surface energies are 0.45 J m⁻² for (111), 1.28 J m⁻² for (100), and 0.86 J m⁻² for (110).⁶³

The pure ZrO₂ and Gd₂O₃ systems also presented the truncated octahedron-like shape with dominance of the (111) surface. In this case, comparatively larger (100) and smaller (111) facets are observed, which is indicative of a different (111)/(100) surface energy ratio compared to the CeO₂ system.⁶⁴ The edges in the polyhedron created by the ZrO₂ nanoparticles are also more irregular, with substantial oxygen rearrangements. In particular, Gd₂O₃ is found in a configuration that has a crystal grain boundary, which originated from a phase transition as a consequence of a double nucleation in the critical temperature. As a consequence, the Gd₂O₃ nanoparticles do not show a perfect octahedron shape, but the combination of two pieces of different truncated octahedra, Fig. 1a. Ce₂O₃ and La₂O₃ present an irregular scheme of faces in the surface, with a more complex final geometry with no clear face definition. Thus, this evidence indicates that these two systems are in different crystal structures than the cubic fluorite structure reported for CeO₂. Furthermore, those structures have weaker non-local ordering, which will be proved by X-ray analysis in the section below.

All nanoalloys ended the annealing process with a nanoparticle shape similar to the pure CeO₂, i.e., in a truncated octahedron-like structure, Fig. 1b. Our results show that the behavior of the mixed systems in the case of ZrO₂ and Gd₂O₃ alloys with CeO₂ can be explained by the combination of the pure system behavior. A decrease in the concentration of cerium species results in an increase in the surface area of (100) faces and a decrease in the (111) ones (even with this decrease, the (111) faces are the largest ones in all cases). Ce₂O₃ and La₂O₃ alloys have the same geometry and face structure as the pure CeO₂, which implies that the more amorphous systems just follow the crystal structure and arrangement imposed by CeO₂. In particular, (CeO₂)_0.25(La₂O₃)_0.75/2 also presents a grain boundary like the one observed in pure Gd₂O₃. Grain boundary formation was observed in two simulations with different materials; however, we cannot conclude that this phenomenon is related to the atomic species itself. In our observation, it is associated with multiple nucleation sites along the phase transition, which is more likely in larger particles. In agreement, this behavior was not observed for any of the 500 cation particles.

3.2 Atomic local ordering and arrangement

From Fig. 2, we can conclude that the CeO₂, Gd₂O₃ and ZrO₂ particles have truncated octahedra-like structures and display almost the same crystalline XRD signature. In the case of the CeO₂ and ZrO₂ systems, XRD confirms the fluorite cubic structure with space group Fm3̄m, associated with the well ordered bulk in the nanoparticle core, and different lattice parameters. The small peaks near 20° and 45° also show that the Gd₂O₃ nanoparticle is, in fact, a more intricate cubic structure with space group Ia3̄, as expected for this material.⁶⁵

Fig. 2 X-ray diffraction pattern obtained for the nanoparticles based on our 864 cation model in 500 ps molecular dynamics at 300 K starting from the 0 K geometry. The simulation aims to capture the signature with a Cu X-ray source (λ = 1.541838 Å).

Pure Ce₂O₃ and La₂O₃ have a common behavior that is not the traditional one shown by any other sample, including the mixtures. Their XRD patterns, Fig. 2, display broader peaks and no strong signal for larger angles. This by itself is a signal that the particle does not show long range order. Fig. 3(a) shows the partial Ce³⁺–O²⁻g(r) for our obtained Ce₂O₃ nanoparticle, comparing it to the same three crystals. In this case, the structure does not show exactly any of the patterns, but it has a strong resemblance with the cubic phase, e.g., the same Ia3̄ of the Gd₂O₃ nanoparticle. But, in any of the cases, it is impossible to attribute a specific crystal structure to those two particles.

Fig. 3 (a) Comparative plot for partial Ce³⁺–O²⁻g(r) for our final Ce₂O₃ nanoparticles obtained from molecular dynamics. The points indicated in the plot are the position and intensities for the pure infinite crystal in three different phases. (b) Visualization of low strain cations obtained for our final Ce₂O₃ and La₂O₃ models, using OVITO.⁵¹ The atoms shown as green are classified as cubic crystal atoms, while the red ones are trigonal.

A direct inspection of those two systems reveals an intricate structure with many local crystallographic domains. In Fig. 3(b), we used the OVITO⁵¹ analysis tool to exclude all oxygen atoms of the structure and analyse only the cation sublattice. The atoms with low strain, except for the surface atoms, were classified as FCC (green atoms) or HCP (red atoms) sub-lattices, which correspond directly to the cubic Ia3̄ and trigonal P3̄m1 phases of Ce₂O₃, respectively. Additionally, this analysis shows that 88.4% of the cerium atoms and 77.3% of the lanthanum atoms are in the surface or in non-crystalline sites. As seen in Fig. 1, the surface shows local order and it is responsible for 36.23% of cerium and 34.95% of lanthanum in those systems (computed using the cation FCC sub-lattice coordination number, considering as part of the surface the atoms with coordination smaller or equal to 10). That points us to the conclusion that the grain boundaries and amorphous internal regions are responsible for the larger part of the cations in the structure. In conclusion to that, the pure Ce₂O₃ and La₂O₃ nanoparticles obtained in this work have a complex topology, with grains of cubic and trigonal structures and amorphous grain boundaries.

For the mixed systems, we see the same behavior as for the pure CeO₂ system. Even with only 25% of the cations as Ce⁴⁺, the nanoparticle displays the truncated octahedron shape, Fig. 1. Also, both g(r) and the XRD (Fig. 2) of those systems indicate that even locally the system behaves as CeO₂, regarding the local and global ordering. The major changes are the lattice parameters, which adjust to better suit the mixed species, and the loss of precision in the partial g(r) for the O²⁻–O²⁻ interaction, which indicates a less ordered oxygen sub-lattice.

3.3 Phase transition

The annealing process simulated here covers a wide range of temperatures, starting from 8500 K and going up to 0 K. As discussed before, those systems start in completely amorphous configurations, which led to different crystal arrangements by the simulations. In order to do so, we found that the system undergoes a phase transition that can be associated with the liquid/amorphous to solid/crystal transition, Fig. 1. In this section, we discuss the process of the phase transition of the nanoparticles.

Fig. 4 presents the overall phase transition observed in all samples, which is illustrated by the pure CeO₂ nanoparticles. The phase transition starts by a nucleation process that takes place in the atomic positions near to the surface of the nanoparticle. Those atoms are in a region that provides more disorder, since the surface yields additional freedom and mobility to the atoms. This mobility, as the temperature decreases, helps the atoms to find the optimal positions, e.g., the lowest energy coordinates, in relation to the others and it starts the nucleation. Once the initial crystal seed is formed, the structure is completely converted to its final crystalline form in less than 100 ps. In the already discussed cases in which the final crystal has interfaces or grain boundaries, the geometry is defined at this time, and the grain boundaries consolidate in the same time window.

Fig. 4 Cations of the pure CeO₂ nanoparticle, colored by the neighborhood crystalline structure. Green atoms are located in the face-centered cubic sub-lattice sites, which are expected in the bulk fluorite CeO₂ crystal, while white ones are in an amorphous or incomplete neighbourhood.

The phase transition temperature observed here is not enough to estimate a specific value for the systems, since there is only one experiment and the compositions of the crystal in the nucleation seed can influence its value. Hence, we studied the phase transition temperature in two ways, both exploring the fluctuation of extensive variables. As both ways describe the same phenomena, we present here the energy evolution over temperature, Fig. 5, and in the ESI† we present the density evolution over time in the simulation, Section 6.3. As shown in Fig. 5, the pure CeO₂ has the highest phase transition temperature in Fig. 5, and hence the decrease in the Ce⁴⁺ ratio in the nanoparticles leads to a decrease in the temperatures at which the transitions happen. The observed phase transition temperatures are not linear with the nanoparticle composition, fluctuating around this behavior. This pattern is seen for all nanoparticles, except for the CeO₂–ZrO₂ mixtures. While the other systems show a correlation between the phase transition temperature and the relative Ce⁴⁺ concentrations in the mixture, in the Zr⁴⁺-based systems, the 50% nanoparticle shows a phase transition at higher temperature than the 25% and 75%.

Fig. 5 Potential energy of the nanoparticle (relative to the 0 K configuration) as a function of instantaneous temperature for the CeO₂–Ce₂O₃ 864 cation nanoparticles. Analogous plots for all systems are provided as ESI.†

3.4 Relative stability via excess energy

To address the relative stability of those systems, we calculated the excess energy for the mixed systems, E^A+B_exc, which is defined as follows,where E^A+B_pot is the calculated total potential energy for the mixed system A + B, while the E^A_pot and E^B_pot terms are the total potential energies obtained for the same size systems containing only the species A and B, respectively. C^A and C^B are the relative concentrations of species A and B in the mixed system A + B. Thus, for pure system A (or B), we obtained E_exc = 0. As a single geometry was obtained here for each system, the values presented in Fig. 6 are an approximation of the real behavior, which can show smaller excess energies if more stable structures were found.

Fig. 6 Excess energy (E_exc) calculated from the structures at zero temparature.

The apparent anomaly observed in the phase transition of the CeO₂–ZrO₂ mixtures is here explained by the lower excess energy at the 50% composition, when compared with the 25% and 75% mixtures. This behaviour seems to be a characteristic of the scale, since previous DFT studies have shown that CeO₂–ZrO₂ mixtures are stable on smaller scales and have a non trivial dependency between concentration and excess energy.²⁵ In the negative excess energy cases, the minimal excess energy observed for each composition is −75.9 eV in (CeO₂)_0.25(Ce₂O₃)_0.75/2, −10.4 eV in (CeO₂)_0.5(Gd₂O₃)_0.5/2, and −78.8 eV in (CeO₂)_0.5(La₂O₃)_0.5/2. So, the observed stability of these solid mixtures follows the order La³⁺, Ce³⁺, Gd³⁺, and Zr⁴⁺. In this same order, the Buckingham potential loses strength (i.e., A and C decrease) and the Coulomb interaction contribution to the total energy of the system increases. Thus, in our results, the relative strength of the Coulomb term in the atomic interactions has a negative correlation with the solution stability.

3.5 Oxygen vacancy distribution

So far, we have analyzed these oxide structures from the point of view of the morphology and the involved chemical species. However, if we are interested in applying this material in catalysis, we also have to focus on its ability to accept oxygen atoms in its structure. To analyze this, we studied the vacancies from both the point of view of their positions and the temporal evolution of their population in the nanoparticles. First, we found a very small number of oxygen vacancies in the pure CeO₂ structures (17 in the final structure at 0 K), which are related to surface edges and defects on the (100) surfaces, as shown in Fig. 8. The remaining systems show more vacancies in general, and this can be seen in Fig. 7a. The Ce₂O₃ and La₂O₃ systems are the ones that present the largest number of vacancies, followed by the cubic Gd₂O₃ and ZrO₂ systems. The mixtures display a behavior that is intermediate between the pure systems with small deviations, which indicates that the mixing of the materials induces the formation of vacancies in the nanoparticles.

Fig. 7 (a) Dependency of the number of oxygen vacancies in the final structure on the relative concentration of the cations in mixed nanoparticles. (b) Evolution of the vacancy percentile in the oxygen sites of the pure structures. The total number of oxygen sites is also evolving, and is presented in Fig. S36 (ESI†).

Fig. 8 Vacancy (represented as green spheres) distribution on CeO₂ pure systems at three different temperatures. The Ce and O atoms are represented by gray and red dots, respectively, for ease of visualization of the vacancy positions.

The same behavior can be seen in the evolution of the ratio between the number of vacancies and the number of oxygen sites in the structures seen in Fig. 7b. Crystallization generates a discontinuity in the occupation of vacancies for most of the systems, with the exception of ZrO₂, which is indicative of the low or no mobility of the oxygen atoms in the system. The M₂O₃ systems displayed a small discontinuity, which showcases the oxygen mobility and the local rearrangement to fill the vacancies in the phase transition. Comparing this graph and the evolution of the number of oxygen sites, we found that Ce₂O₃ and La₂O₃ do not show discontinuities in the number of oxygen sites in the phase transition, which is a sign of a smooth transition from a liquid to a more amorphous solid in those materials.

CeO₂ and Gd₂O₃ have transitions that lower the number of oxygen sites, while in the ZrO₂ nanoparticle this number increases. In both cases, there is a strong reorganization of the system that indicates crystallization. The slow increase in the vacancy site count even in the liquid phase is again a sign of low oxygen mobility in ZrO₂. The increase seen in the phase transition shows that the system undergoes strong reorganization that gives the oxygen and cations the “chance” to rearrange themselves and go to a number of sites similar to the ones obtained continuously by the other systems in the liquid phase. This strong rearrangement is in agreement with the energy fluctuation in the phase transition, e.g., Fig. 5, where ZrO₂ is the largest one after CeO₂.

3.6 Oxygen mobility

As described above, we analyzed the energy barriers when moving an oxygen atom to a vacancy in different regions in the pure CeO₂ nanoparticle structure as in the schematic representation of Fig. 9b. In the first case, when moving an oxygen atom located in the first layer to a vacancy located at the second layer of the nanoparticle, this motion is energetically unfavorable. Hence, the displacement of anions from the second layer of oxygen to fill vacancies on the surface of the particle is energetically favorable, which can be seen in Fig. 9a(i).

Fig. 9 (a) Energy barriers for oxygen atom migration to a vacancy. The US was calculated in three different regions of the particle, where the oxygen atoms were moved from their initial position to vacant positions in which they were the first neighbors. All initial positions coincide with the origin of the plot. (b) Schematic representation of the moving oxygen atom trajectories for the US simulation. (c) Radial distribution function obtained for the oxygen atom migration in the core indicating the Ce atom and O atom distribution relative to the moving O atom. The coordinate ξ is the distance traveled by the moving O atom in nm, which represents different windows of the US simulations.

We also considered the energy barrier to move an oxygen atom located in the third oxygen layer to a vacancy site in the second layer, Fig. 9a(ii). The energy difference between the start position, i.e., the origin of the graph, and the end position at ξ ≈ 2.15 Å is close to 1.5 eV. Although it may be overestimated due to the limitations of the US technique, this value allows a qualitative interpretation that shows that the migration of the O anions from the third oxygen layer to a vacancy site of the particle is also energetically unfavorable. In this scenario, O atoms in the second oxygen layer are strongly susceptible to local rearrangement if located near an oxygen vacancy.

Finally, we considered the displacement of an oxygen atom to a vacancy site that was artificially created by removing the oxygen atom nearest to the geometric center of the nanoparticle. Even in this case, we notice an energy barrier of the order of 1.5 eV, but, unlike the previous situations, there is a peak of energy during the displacement, Fig. 9a(iii). This result shows that even at the center of the nanoparticle there is a potential difference that could induce atomic rearrangements. But, the presence of a strong potential barrier that originates from the crystalline order itself lets us conclude that this displacement is improbable at low temperatures.

An additional concern when using US to determine an energy barrier is to check the integrity of the crystalline structure. In this case, we employed radial distribution functions to check for possible changes in the structure of the nanoparticle during the displacement of the moving oxygen atom in the core, as we see in Fig. 9c. Because it is a periodic structure, the g(r)s show that there is not much difference between the starting position and the ending position. In Fig. 9c, we notice for the O–Ce g(r) that the start and end positions (ξ = 0.000 and 0.233) have 3 peaks at almost equal distances. During the trajectory there is an almost perfect symmetry in the first peaks. For the O–O g(r), we notice a symmetry in the first 2 peaks, with these maxima occurring in the initial and final positions. This result shows that the integrity of the crystalline structure remains cohesive during the umbrella sampling simulations.

3.7 Electrostatic potential

We analysed the map of the electrostatic potential in the nanoparticles in relation to the M/Ce⁴⁺ molar ratio (M = Ce³⁺, Gd³⁺, La³⁺, Zr⁴⁺) using the electrostatic potential in the range from −0.52 to 2.08, Fig. 10. The red regions in the nanoparticle structures represent low electrostatic potential regions, e.g., large concentrations of negative charges in the structure, while the blue-colored atoms indicate the most positive regions.

Fig. 10 Electrostatic potential map obtained with the particle-mesh Ewald method analysis for all nanoparticle structures at 0 K in relation to the M/Ce⁴⁺ ratio.

Migani and coworkers investigated the vacancy formation energy in small (CeO₂)_n nanoclusters where n = 21, 30, 40 and 80 using DFT calculations. They found that the energy formation of an O vacancy is related to the electrostatic potential of the material, and then the reactivity can be rationalized using the electrostatic potential, in particular, at the surface region.⁶⁶ Along this line, Sayle and coworkers have employed the map of the electrostatic potential as a descriptor to identify reactivity “fingerprints” on Ce³⁺/Ce⁴⁺ nanoparticles obtained by MD simulations. They found that the Ce³⁺/Ce⁴⁺ ratio and morphology play crucial roles in the oxygen vacancy formation energy and consequently in the oxygen storage capability and reactivity of the material.^13,67

The concentration and predominance of blue or red regions on the nanoparticle surface, Fig. 10, is related not only to the point charge of each atom, q, but also to the van der Waals radius, r. As mentioned before, both parameters, q and r, are required to calculate the reciprocal sum of the PME and the charge q is distributed on the surface of the sphere with radius r. Thus, atoms with small radius contribute with higher charge density and larger homogeneous regions. We would like to highlight that the electrostatic potential in Fig. 10 is the contribution of all atoms in the nanoparticle and not only the surface atoms.

The presence of M³⁺ cations in the nanoparticles decreases the electrostatic potential, which is related to the larger populations of M³⁺ cations in the surface and agrees with the greater reactivity of Ce³⁺/Ce⁴⁺, e.g., lower O²⁻ vacancy energy and higher oxygen storage capacity, in relation to the neat Ce⁴⁺ particles.^13,67 Beyond the influence of the M³⁺ cations, we also see a larger occurrence of negative regions in the edges and defects of the nanoparticles, which can be related to the lower bound of the oxygen in these regions.⁶⁶

On the other hand, the CeO₂–ZrO₂ nanoparticles present the largest and homogeneous blue regions on the surface, e.g., positive regions, and it can be understood by relating the Zr⁴⁺ high charge q = +4 with the small zirconium ionic radius, e.g., r = 1.75 Å. Moreover, as the Zr⁴⁺/Ce⁴⁺ molar ratio increases, the blue regions get more dominant over the negative potential, which is also the consequence of the high density of the positive charge generated by Zr⁴⁺ in comparison with q = +4 and r = 2.04 Å of the Ce⁴⁺ cation. Hence, the MEP in Fig. 10 indicates that the M³⁺/Ce⁴⁺ nanoparticles may present interesting activity in reactions involving C–H bond activation in hydrocarbons, a key process in conversion of CH species to CO species.⁶⁸

On the other hand, the M⁴⁺/Ce⁴⁺ nanoparticles may be promising materials for conversion of O₂-species to molecular O₂. This can in principle be done by adsorbing ions on the particle surface (in the positive potential region) leading to oxidation, as in O²⁻ → O⁰ + 2e⁻. As an example, this is an already explored application of CeO₂ in the context of water splitting.^69,70 Thus, with the classical point of view of the electrostatic potential, by combining red and blue regions nearby, a nanoparticle could be for instance an interesting candidate as a catalyst for oxi-reduction reactions.^68,71,72

4 Insights in the role of temperature and composition in mixed CeO₂-based nanoalloys

Within this section, we discuss the general remarks and the knowledge that arises from our results. First of all, the shapes of the particles are unchanged by temperature effects below the phase transition point, and only change by the effect of the particle composition. The truncated octahedron shape seen in the pure ceria is inherited by the mixtures, even with the small 25% CeO₂ concentration. In the particular case of the Ce₂O₃ and La₂O₃ systems, the final particles have a less regular morphology that follows the more irregular (locally trigonal) crystal arrangement seen in this work.

In mixed particles, the presence of Ce⁴⁺ and the lower melting point of CeO₂ lead to the nucleation and growth of a crystal that follows the cubic fluorite structure of pure ceria, with vacancies appearing both by temperature effects and the presence of 3+ species. The pure systems present locally the same morphology expected from the bulk configurations, but for Ce₂O₃ and La₂O₃ the balance between bulk and surface stress favored a less globally ordered structure with irregular faceting. For pure La₂O₃, this can be associated with the weaker bond energy expected from lanthanum to oxygen due to its small electronegativity, which is the smallest among all the metals studied here (1.10 on the Pauling scale). Ce has a similar electronegativity (1.12 on the Pauling scale), which also results in the weakening of the Ce³⁺–O²⁻ bond in comparison to Ce⁴⁺–O²⁻.

Even without having an impact on the particle shape, the temperature effects below the phase transition have a strong connection with the oxygen mobility and distribution. From the evolution of the number of vacancies on the particles, we have shown that different oxides show different resistance for oxygen mobility, in which CeO₂ presents the largest mobility and ZrO₂ the smallest. The intermediate mobility seems to be similar in all the three 3+ cations studied here, and more detailed experiments are necessary to evaluate their exact relation. When combined with the relative concentration distribution of the cations on the nanoparticle, we can conclude that the CeO₂–ZrO₂ mixture will behave more linearly with respect to mobility, since the cations are more homogeneously distributed in the particle. The other mixtures will display a more nonlinear behavior since the surface is richer in 3+ cations in relation to the bulk.

This last point connects directly with our observations on the excess energy of those systems. The high excess energy shown in the ZrO₂ mixtures indicates that, if the systems were annealed at a slower rate, the cations could segregate even more. In complete opposition, we observed a small segregation of the 3+ cations towards the surface, even with negative excess energy. This is expected, since the 3+ cations do not need the same coordination as the 4+ ones to minimize the system energy. This surface preference for 3+ cations can be useful since it creates the possibility of a high dopant concentration in the chemically active surface using a small M³⁺/Ce⁴⁺ molar ratio.

The Ce⁴⁺ abundance in the bulk of CeO₂–M₂O₃ systems is an explanation for why those mixtures can maintain the fluorite crystal structure even at a lower relative Ce⁴⁺ concentration. The other explanation is that CeO₂ has a higher phase transition temperature in comparison to the other oxides and can nucleate first in the simulations. The surface of the final CeO₂ nanoparticle has an interesting mixture of positive and negative electrostatic potential, in which the (111) surface is more positive and the (100) more negative. As the mixtures tend to maintain the geometry of the CeO₂ particle, we observe that the addition of 3+ cations lowers the potential on the (111) face, i.e., makes the surface in general less positive.

The addition of ZrO₂ up to 50% changes the surface only slightly, but higher concentrations generate a particle that has a completely positive surface. This result agrees qualitatively with the observations from Pengpanich et al., in which they observed that the addition of Zr⁴⁺ up to 25% leads to an “optimal” methane conversion performance for the Ce_1−xZr_xO₂ nanoparticle, while a higher concentration leads to the decay of the catalytic performance.⁷³

The positive profile of the surface in Ce⁴⁺/Zr⁴⁺ is strongly associated with the oxygen mobility we observed in this work using the umbrella sampling technique. We observed that to move the O²⁻ anions from the first layer of the particles to the surface is more favorable than moving in the opposite direction. On the other hand, the combination of these results points to a general path of conduction that is only reliable for neutral particles. If the systems are unbalanced in charge, more complex behaviors can appear.

5 Conclusion

We employed classical molecular dynamics simulations and the process of amorphization and recrystallization to investigate the structural formation of mixed CeO₂-based nanoparticles with different concentrations of Ce₂O₃, Gd₂O₃, La₂O₃, and ZrO₂. The studied systems covered two different nanoparticle sizes and 5 compositions for each CeO₂–MO_x system. Our results regarding the nanoparticle morphology revealed that the process of annealing makes the CeO₂ morphology overcome the other oxide ones, which gives to the mixture the fluorite crystal structure in a truncated octahedron-like particle. The CeO₂ concentration is also detrimental to the phase transition in those systems, both in the temperature and in the fluctuation of extensive variables such as density and potential energy. For these three properties, the relation with concentration is closely approximated by a linear combination of isolated species. CeO₂–ZrO₂ is the only solid solution found to have positive excess energy and more abrupt fluctuations in the phase transition temperature as a function of the concentration, which is the opposite of the DFT results obtained for small CeO₂–ZrO₂ clusters.²⁵

Our findings showed that the nanoparticles have a good combination of positive and negative electrostatic potential regions on the surface, with a good modulation of this value with the 3+ or 4+ species. Also, oxygen mobility and vacancies were investigated, revealing the strong potential barriers that exist for oxygen transport in those solids. CeO₂ crystals eliminate most of the vacancies in the annealing process, but the 3+ mixtures retain the vacancies even at low temperature. The sites near the surface are the more energetic ones, but vacancies in the bulk generate stronger energy fluctuations by virtue of atomic rearrangement. These characteristics observed and explained here are signs that particles in this size range may present catalytic properties such as observed in experiments and can be used as more realistic models for computational simulation than few atoms clusters. Therefore, we expect these results to motivate future research on the direct observation of oxi-reduction reactions over those nanoparticles.

Author contributions

João Paulo A. de Mendonça – methodology, software, formal analysis, investigation, data curation, writing – original draft, visualization. Tuanan C. Lourenço – methodology, software, validation, formal analysis, investigation, writing – review & editing. Luis Paulo M. Freitas – formal analysis, writing – review & editing. Anderson A. E. Santo – formal analysis. Gustavo T. Feliciano – formal analysis, writing – review & editing, supervision. Juarez L. F. Da Silva – conceptualization, methodology, resources, data curation, writing – review & editing, visualization, supervision, project administration, funding acquisition.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors gratefully acknowledge support from FAPESP (São Paulo Research Foundation, Grant Numbers 2017/11631-2, 2018/21401-7 and 2019/23681-0), Shell and the strategic importance of the support given by ANP (Brazils National Oil, Natural Gas and Biofuels Agency) through the R&D levy regulation and also to the University of São Paulo for providing the computational resources for the project execution. The authors acknowledge also the Laboratory of Advanced Scientific Computing (University of São Paulo) and the Department of Information Technology – Campus São Carlos, for hosting our cluster. The authors acknowledge also the National Laboratory for Scientific Computing (LNCC/MCTI, Brazil) for providing HPC resources of the SDumont supercomputer, which have contributed to the research results reported within this paper. URL: http://sdumont.lncc.br.

References

1. R. Schmitt, A. Nenning, O. Kraynis, R. Korobko, A. I. Frenkel, I. Lubomirsky, S. M. Haile and J. L. M. Rupp, A review of defect structure and chemistry in ceria and its solid solutions, Chem. Soc. Rev., 2020, 49, 554–592,  10.1039/C9CS00588A.
2. Y. Shan, Y. Liu, Y. Li and W. Yang, A review on application of cerium-based oxides in gaseous pollutant purification, Sep. Purif. Technol., 2020, 250, 117181,  DOI:10.1016/j.seppur.2020.117181.
3. P. Li, X. Chen, Y. Li and J. W. Schwank, A Review on Oxygen Storage Capacity of CeO2-Based Materials: Influence Factors, Measurement Techniques, and Applications in Reactions Related to Catalytic Automotive Emissions Control, Catal. Today, 2019, 327, 90–115,  DOI:10.1016/j.cattod.2018.05.059 , New Insight of Environmental Catalysis (18th NCC).
4. A. Vita, Catalytic Applications of CeO₂-Based Materials, Catalysts, 2020, 10, 576,  DOI:10.3390/catal10050576.
5. P. Bellio, C. Luzi, A. Mancini, S. Cracchiolo, M. Passacantando, L. Di Pietro, M. Perilli, G. Amicosante, S. Santucci and G. Celenza, Cerium oxide nanoparticles as potential antibiotic adjuvant. Effects of CeO₂ nanoparticles on bacterial outer membrane permeability, Biochim. Biophys. Acta, Biomembr., 2018, 1860, 2428–2435,  DOI:10.1016/j.bbamem.2018.07.002.
6. S. Rajeshkumar and P. Naik, Synthesis and biomedical applications of Cerium oxide nanoparticles - A Review, Biotechnol. Rep., 2018, 17, 1–5,  DOI:10.1016/j.btre.2017.11.008.
7. G. Lu, H. Zheng, J. Lv, G. Wang and X. Huang, Review of recent research work on CeO₂-based electrocatalysts in liquid-phase electrolytes, J. Power Sources, 2020, 480, 229091,  DOI:10.1016/j.jpowsour.2020.229091.
8. J. L. F. Da Silva, M. V. Ganduglia-Pirovano, J. Sauer, V. Bayer and G. Kresse, Hybrid Functionals Applied to Rare-Earth Oxides: The Example of Ceria, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 75, 045121,  DOI:10.1103/physrevb.75.045121.
9. M. V. Ganduglia-Pirovano, J. L. F. Da Silva and J. Sauer, Density-Functional Calculations of the Structure of Near-Surface Oxygen Vacancies and Electron Localization on CeO₂(111), Phys. Rev. Lett., 2009, 102, 026101,  DOI:10.1103/physrevlett.102.026101.
10. V. Goldschmidt and L. Thomassen, Crystal structure of natural and synthetic oxides of U, Th and Ce. Skrifter utgitt av det Norske Videnskaps-Akademi i Oslo 1: Matematisk-Naturvidenskapelig Klasse, 1923, 1923, 1–48.
11. H. Whitfield, D. Roman and A. Palmer, X-ray study of the system ThO₂-CeO₂-Ce₂O₃, J. Inorg. Nucl. Chem., 1966, 28, 2817–2825.
12. D. R. Mullins, The Surface Chemistry of Cerium Oxide, Surf. Sci. Rep., 2015, 70, 42–85,  DOI:10.1016/j.surfrep.2014.12.001.
13. T. X. Sayle, M. Molinari, S. Das, U. M. Bhatta, G. Möbus, S. C. Parker, S. Seal and D. C. Sayle, Environment-mediated structure, surface redox activity and reactivity of ceria nanoparticles, Nanoscale, 2013, 5, 6063–6073,  10.1039/c3nr00917c.
14. P. Li, Y. Zhou, Z. Zhao, Q. Xu, X. Wang, M. Xiao and Z. Zou, Hexahedron Prism-Anchored Octahedronal CeO2: Crystal Facet-Based Homojunction Promoting Efficient Solar Fuel Synthesis, J. Am. Chem. Soc., 2015, 137, 9547–9550,  DOI:10.1021/jacs.5b05926.
15. S. Yasmeen, M. R. Khan, K. Park, Y. Cho, J. W. Choi, H.-S. Moon and H.-B.-R. Lee, Preparation of a Hydrophobic Cerium Oxide Nanoparticle Coating with Polymer Binder via a Facile solution route, Ceram. Int., 2020, 46, 12209–12215,  DOI:10.1016/j.ceramint.2020.01.268.
16. J. C. Hamilton, Edge Energies: Atomistic Calculations of a Continuum Quantity, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 125447,  DOI:10.1103/PhysRevB.73.125447.
17. S. F. Li, X. J. Zhao, X. S. Xu, Y. F. Gao and Z. Zhang, Stacking Principle and Magic Sizes of Transition Metal Nanoclusters Based on Generalized Wulff Construction, Phys. Rev. Lett., 2013, 111, 115501,  DOI:10.1103/physrevlett.111.115501.
18. X. J. Zhao, X. L. Xue, Z. X. Guo and S. F. Li, Relative Edge Energy in the Stability of Transition Metal Nanoclusters of Different Motifs, Nanoscale, 2016, 8, 12834–12842,  10.1039/C6NR00486E.
19. J. P. A. de Mendonça, M. G. Monteiro, S. de Andrade Leonel and F. Sato, Modeling Surface Energy in Porous metallic Nanostructures, Eur. Phys. J. B, 2018, 91, 1–12,  DOI:10.1140/epjb/e2018-80655-2.
20. R. W. Tarnuzzer, J. Colon, S. Patil and S. Seal, Vacancy Engineered Ceria Nanostructures for Protection from Radiation-Induced Cellular Damage, Nano Lett., 2005, 5, 2573–2577,  DOI:10.1021/nl052024f.
21. H. Yahiro, Y. Eguchi, K. Eguchi and H. Arai, Oxygen Ion Conductivity of the Ceria-Samarium Oxide System with Fluorite Structure, J. Appl. Electrochem., 1988, 18, 527–531,  DOI:10.1007/BF01022246.
22. J. J. Plata, A. M. Márquez and J. F. Sanz, Understanding the Interplay of Dopants, Interfaces, and Anionic Conductivity in Doped Ceria/Zirconia Heteroepitaxial Structures, Chem. Mater., 2014, 26, 3385–3390,  DOI:10.1021/cm500415h.
23. T. Montini, M. Melchionna, M. Monai and P. Fornasiero, Fundamentals and Catalytic Applications of CeO2-Based Materials, Chem. Rev., 2016, 116, 5987–6041,  DOI:10.1021/acs.chemrev.5b00603.
24. J. Mucelini, R. Costa-Amaral, Y. Seminovski and J. L. F. Da Silva, Ab initio investigation of the formation of ZrO₂-like structures upon the adsorption of Zr_n on the CeO₂(111) surface, J. Chem. Phys., 2018, 149, 244702,  DOI:10.1063/1.5063732.
25. P. Felício-Sousa, J. Mucelini, L. Zibordi-Besse, K. F. Andriani, Y. Seminovski, R. C. Prati and J. L. F. Da Silva, Ab initio Insights Into the Structural, Energetic, Electronic, and Stability Properties of Mixed Ce_nZr₁₅−nO₃₀ Nanoclusters, Phys. Chem. Chem. Phys., 2019, 21, 26637–26646,  10.1039/c9cp04762j.
26. K. Jiang, S. Liu and X. Wang, Low-Thermal-Conductivity and High-Toughness CeO₂-Gd₂O₃ Co-Stabilized Zirconia Ceramic for Potential Thermal Barrier Coating Applications, J. Eur. Ceram. Soc., 2018, 38, 3986–3993,  DOI:10.1016/j.jeurceramsoc.2018.04.065.
27. Z. Tianshu, Ionic Conductivity in the CeO₂–Gd₂O₃ System (0.05 ≤ Gd/Ce ≤ 0.4) Prepared by Oxalate Coprecipitation, Solid State Ionics, 2002, 148, 567–573,  DOI:10.1016/s0167-2738(02)00121-2.
28. V. P. Kumar, Y. Reddy, P. Kistaiah, G. Prasad and C. V. Reddy, Thermal and Electrical Properties of Rare-Earth Co-Doped Ceria Ceramics, Mater. Chem. Phys., 2008, 112, 711–718,  DOI:10.1016/j.matchemphys.2008.06.030.
29. E. Andrievskaya, O. Kornienko, A. Sameljuk and A. Sayir, Phase Relation Studies in the CeO₂-La₂O₃ System at 1100-1500 °C, J. Eur. Ceram. Soc., 2011, 31, 1277–1283,  DOI:10.1016/j.jeurceramsoc.2010.05.024 , special issue: Directionally solidified eutectic and advanced ceramics.
30. P. A. Žguns, A. V. Ruban and N. V. Skorodumova, Phase Diagram and Oxygen–Vacancy Ordering in the CeO₂–Gd₂O₃ System: a Theoretical Study, Phys. Chem. Chem. Phys., 2018, 20, 11805–11818,  10.1039/c8cp01029c.
31. G. Balducci, J. Kašpar, P. Fornasiero, M. Graziani and M. S. Islam, Surface and Reduction Energetics of the CeO₂-ZrO₂ Catalysts, J. Phys. Chem. B, 1998, 102, 557–561,  DOI:10.1021/jp972400n.
32. N. N. Bulgakov, V. Sadykov, V. V. Lunin and E. Kemnitz, React. Kinet. Catal. Lett., 2002, 76, 103–110,  DOI:10.1023/a:1015621613212.
33. I. Z. Koleva, H. A. Aleksandrov, K. M. Neyman and G. N. Vayssilov, Preferential Location of Zirconium Dopants in Cerium Dioxide Nanoparticles and the Effects of Doping on Their Reducibility: A DFT Study, Phys. Chem. Chem. Phys., 2020, 22, 26568–26582,  10.1039/D0CP05456A.
34. K. Reed, A. Cormack, A. Kulkarni, M. Mayton, D. Sayle, F. Klaessig and B. Stadler, Exploring the properties and applications of nanoceria: Is there still plenty of room at the bottom?, Environ. Sci.: Nano, 2014, 1, 390–405,  10.1039/c4en00079j.
35. T. X. T. Sayle, S. C. Parker and D. C. Sayle, Shape of CeO₂ nanoparticles using simulated amorphisation and recrystallisation, Chem. Commun., 2004, 2438–2439,  10.1039/b408752f.
36. T. X. Sayle, F. Caddeo, X. Zhang, T. Sakthivel, S. Das, S. Seal, S. Ptasinska and D. C. Sayle, Structure-Activity Map of Ceria Nanoparticles, Nanocubes, and Mesoporous Architectures, Chem. Mater., 2016, 28, 7287–7295,  DOI:10.1021/acs.chemmater.6b02536.
37. A. Jeyaranjan, T. S. Sakthivel, M. Molinari, D. C. Sayle and S. Seal, Morphology and Crystal Planes Effects on Supercapacitance of CeO₂ Nanostructures: Electrochemical and Molecular Dynamics Studies, Part. Part. Syst. Charact., 2018, 35, 1–9,  DOI:10.1002/ppsc.201800176.
38. T. X. T. Sayle, P. E. Ngoepe and D. C. Sayle, Generating structural distributions of atomistic models of Li2O nanoparticles using simulated crystallisation, J. Mater. Chem., 2010, 20, 10452,  10.1039/c0jm01580f.
39. R. S. Ledwaba, D. C. Sayle and P. E. Ngoepe, Atomistic Simulation and Characterization of Spinel Li_1+xMn₂O₄ (0 ≤ x ≤ 1) Nanoparticles, ACS Appl. Energy Mater., 2020, 3, 1429–1438,  DOI:10.1021/acsaem.9b01870.
40. A. H. Tamboli, A. A. Chaugule, S. W. Gosavi and H. Kim, Ce_xZr_1-xO₂ Solid Solutions for Catalytic Synthesis of Dimethyl Carbonate from CO₂: Reaction Mechanism and the Effect of Catalyst Morphology on Catalytic Activity, Fuel, 2018, 216, 245–254,  DOI:10.1016/j.fuel.2017.12.008.
41. Y. Shu, J. Bao, S. Yang, X. Duan and P. Zhang, Entropy-Stabilized Metal-CeO_x Solid Solutions for Catalytic Combustion of Volatile Organic Compounds, AIChE J., 2020, 67, 1–13,  DOI:10.1002/aic.17046.
42. B.-H. Kim, J. Kullgren, M. J. Wolf, K. Hermansson and P. Broqvist, Multiscale Modeling of Agglomerated Ceria Nanoparticles: Interface Stability and Oxygen Vacancy Formation, Front. Chem., 2019, 7, 203,  DOI:10.3389/fchem.2019.00203.
43. R. Álvarez Asencio, R. W. Corkery and A. Ahniyaz, Solventless Synthesis of Cerium Oxide Nanoparticles and Their Application in UV Protective Clear Coatings, RSC Adv., 2020, 10, 14818–14825,  10.1039/D0RA01710H.
44. S. Carvajal, M. Perramón, G. Casals, D. Oró, J. Ribera, M. Morales-Ruiz, E. Casals, P. Casado, P. Melgar-Lesmes, G. Fernández-Varo, P. Cutillas, V. Puntes and W. Jiménez, Cerium Oxide Nanoparticles Protect against Oxidant Injury and Interfere with Oxidative Mediated Kinase Signaling in Human-Derived Hepatocytes, Int. J. Mol. Sci., 2019, 20,  DOI:10.3390/ijms20235959.
45. L. Minervini, M. O. Zacate and R. W. Grimes, Defect cluster formation in M₂O₃-doped CeO₂, Solid State Ionics, 1999, 116, 339–349,  DOI:10.1016/S0167-2738(98)00359-2.
46. R. A. Buckingham, The Classical Equation of State of Gaseous Helium, Neon and Argon, Proc. R. Soc. A, 1938, 168, 264–283.
47. H. Hojo, T. Mizoguchi, H. Ohta, S. D. Findlay, N. Shibata, T. Yamamoto and Y. Ikuhara, Atomic Structure of a CeO₂ Grain Boundary: The Role of Oxygen Vacancies, Nano Lett., 2010, 10, 4668–4672,  DOI:10.1021/nl1029336.
48. G. Balducci, M. S. Islam, J. Kašpar, P. Fornasiero and M. Graziani, Reduction Process in CeO₂-MO and CeO₂-M₂O₃ Mixed Oxides: A Computer Simulation Study, Chem. Mater., 2003, 15, 3781–3785.
49. T. X. Sayle, S. C. Parker and C. R. Catlow, The role of oxygen vacancies on ceria surfaces in the oxidation of carbon monoxide, Surf. Sci., 1994, 316, 329–336,  DOI:10.1016/0039-6028(94)91225-4.
50. S. Plimpton, Fast Parallel Algorithms for Short-Range Molecular Dynamics, J. Comput. Phys., 1995, 117, 1–19,  DOI:10.1006/jcph.1995.1039.
51. A. Stukowski, Visualization and Analysis of Atomistic Simulation Data With OVITO–the Open Visualization Tool, Modell. Simul. Mater. Sci. Eng., 2019, 18, 015012,  DOI:10.1088/0965-0393/18/1/015012.
52. S. P. Coleman, D. E. Spearot and L. Capolungo, Virtual diffraction analysis of Ni [010] symmetric tilt grain boundaries, Modell. Simul. Mater. Sci. Eng., 2013, 21, 055020,  DOI:10.1088/0965-0393/21/5/055020.
53. G. Patey and J. Valleau, The Free Energy of Spheres with Dipoles: Monte Carlo with Multistage Sampling, Chem. Phys. Lett., 1973, 21, 297–300,  DOI:10.1016/0009-2614(73)80139-3.
54. G. Torrie and J. Valleau, Nonphysical Sampling Distributions in Monte Carlo Free-Energy Estimation: Umbrella Sampling, J. Comput. Phys., 1977, 23, 187–199,  DOI:10.1016/0021-9991(77)90121-8.
55. A. Blackman and L. R. Gahan, Aylward and Findlay's SI chemical data, John Wiley & Sons, 2014.
56. M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess and E. Lindahl, GROMACS: High Performance Molecular Simulations Through Multi-Level Parallelism from Laptops to Supercomputers, SoftwareX, 2015, 1–2, 19–25,  DOI:10.1016/j.softx.2015.06.001.
57. J. S. Hub, B. L. de Groot and D. van der Spoel, g_wham—A Free Weighted Histogram Analysis Implementation Including Robust Error and Autocorrelation Estimates, J. Chem. Theory Comput., 2010, 6, 3713–3720,  DOI:10.1021/ct100494z.
58. J. Koettgen, S. Grieshammer, P. Hein, B. O. H. Grope, M. Nakayama and M. Martin, Understanding the Ionic Conductivity Maximum in Doped Ceria: Trapping and Blocking, Phys. Chem. Chem. Phys., 2018, 20, 14291–14321,  10.1039/c7cp08535d.
59. A. Aksimentiev and K. Schulten, Imaging alpha-Hemolysin with Molecular Dynamics: Ionic Conductance, Osmotic Permeability, and the Electrostatic Potential Map, Biophys. J., 2005, 88, 3745–3761,  DOI:10.1529/biophysj.104.058727.
60. W. Humphrey, A. Dalke and K. Schulten, VMD: Visual molecular dynamics, J. Mol. Graphics, 1996, 14, 33–38,  DOI:10.1016/0263-7855(96)00018-5.
61. M. Baalousha, P. L. Coustumer, I. Jones and J. R. Lead, Characterisation of structural and surface speciation of representative commercially available cerium oxide nanoparticles, Environ. Chem., 2010, 7, 377,  DOI:10.1071/en10003.
62. S. Takami, S. Ohara, T. Adschiri, Y. Wakayama and T. Chikyow, Continuous Synthesis of Organic–Inorganic Hybridized Cubic Nanoassemblies of Octahedral Cerium Oxide Nanocrystals and Hexanedioic Acid, Dalton Trans., 2008, 5442,  10.1039/b801099d.
63. Y. Jiang, J. B. Adams and M. van Schilfgaarde, Density-Functional Calculation of CeO₂ Surfaces and Prediction of Effects of Oxygen Partial Pressure and Temperature on Stabilities, J. Chem. Phys., 2005, 123, 064701.
64. U. Castanet, C. Feral-Martin, A. Demourgues, R. L. Neale, D. C. Sayle, F. Caddeo, J. M. Flitcroft, R. Caygill, B. J. Pointon, M. Molinari and J. Majimel, Controlling the {111}/{110} Surface Ratio of Cuboidal Ceria Nanoparticles, ACS Appl. Mater. Interfaces, 2019, 11, 11384–11390,  DOI:10.1021/acsami.8b21667.
65. M. Scavini, M. Coduri, M. Allieta, M. Brunelli and C. Ferrero, Probing Complex Disorder in Ce_1-xGd_xO_2-x/2 Using the Pair Distribution Function Analysis, Chem. Mater., 2012, 24, 1338–1345,  DOI:10.1021/cm203819u.
66. A. Migani, G. N. Vayssilov, S. T. Bromley, F. Illas and K. M. Neyman, Dramatic Reduction of the Oxygen Vacancy Formation Energy in Ceria Particles: A Possible Key to their Remarkable Reactivity at the Nanoscale, J. Mater. Chem., 2010, 20, 10535–10546,  10.1039/c0jm01908a.
67. T. X. T. Sayle, M. Cantoni, U. M. Bhatta, S. C. Parker, S. R. Hall, G. Möbus, M. Molinari, D. Reid, S. Seal and D. C. Sayle, Strain and architecture-tuned reactivity in ceria nanostructures; Enhanced catalytic oxidation of CO to CO₂, Chem. Mater., 2012, 24, 1811–1821,  DOI:10.1021/cm3003436.
68. M. M. Nair and S. Abanades, Solid-State Redox Kinetics of CeO2 in Two-Step Solar CH4 Partial Oxidation and Thermochemical CO2 Conversion, Catalysts, 2021, 11, 723.
69. J. R. Scheffe and A. Steinfeld, Thermodynamic Analysis of Cerium-Based Oxides for Solar Thermochemical Fuel Production, Energy Fuels, 2012, 26, 1928–1936,  DOI:10.1021/ef201875v.
70. R. R. Bhosale, G. Takalkar, P. Sutar, A. Kumar, F. AlMomani and M. Khraisheh, A decade of ceria based solar thermochemical H2O/CO2 splitting cycle, Int. J. Hydrogen Energy, 2019, 44, 34–60,  DOI:10.1016/j.ijhydene.2018.04.080.
71. M. Capdevila-Cortada, M. García-Melchor and N. López, Unraveling the Structure Sensitivity in Methanol Conversion on CeO₂: A DFT+U Study, J. Catal., 2015, 327, 58–64,  DOI:10.1016/j.jcat.2015.04.016.
72. G. Righi, R. Magri and A. Selloni, Methane Activation on Metal-Doped (111) and (100) Ceria Surfaces with Charge-Compensating Oxygen Vacancies, J. Phys. Chem. C, 2020, 124, 17578–17585,  DOI:10.1021/acs.jpcc.0c03320.
73. S. Pengpanich, V. Meeyoo, T. Rirksomboon and K. Bunyakiat, Catalytic oxidation of methane over CeO2-ZrO2 mixed oxide solid solution catalysts prepared via urea hydrolysis, Appl. Catal., A, 2002, 234, 221–233,  DOI:10.1016/s0926-860x(02)00230-2.

---

Footnote

† Electronic supplementary information (ESI) available: We included several additional details, namely, generation of the initial mixed oxide structures, the amorphization protocol, recrystallization tests, relaxation at zero temperature, etc. Furthermore, all the results obtained for the smallest particles (500 cations) are also reported. See DOI: 10.1039/d1ma00543j

---