Structure and reducibility of yttrium-doped cerium dioxide nanoparticles and (111) surface

Abstract

Using periodic density functional calculations, we studied the local structure and preferred locations of yttrium cations and oxygen vacancies in Y-doped cerium dioxide. We employed three kinds of models – a slab of the CeO₂(111) surface and two ceria nanoparticles of different sizes and shapes. In the slab models, which represent the (111) surface of ceria and the corresponding extended terraces on the facets of its nanoparticles, Y³⁺ cation dopants were calculated to be preferentially located close to each other. They tend to surround a subsurface oxygen vacancy that forms to maintain the charge balance. Such general behavior was not found for the nanoparticle models, in which structural flexibility and the presence of various low-coordinated surface centers seem to be crucial and suppress most of the trends. Configurations with four Y³⁺ cations were calculated to be particularly stable when they combined two of the most stable configurations with two Y³⁺ cations. However, no clear trend was found regarding the preferential spatial distribution of the Y³⁺ pairs – they can be stable both in isolation and close to each other. In general, doping by yttrium does not notably change the reducibility of ceria systems but selectively facilitates the formation of oxygen vacancies at the ceria surface in comparison with pristine ceria. Yttrium cations also slightly increase the basicity of the nearby oxygen centers with respect to a stoichiometric ceria surface.

---

1. Introduction

In recent years, ceria-based materials have excited tremendous interest from both fundamental scientific and commercial points of view. In particular, ceria has attracted much attention because of its crucial role as a component of catalysts for numerous catalytic applications. The ever-growing interest is demonstrated by the vast numbers of publications starting from the review by Trovarelli in 1996,¹ which was followed by two books in 2002 (ref. 2) and 2013 (ref. 3) and many reviews,^4–6 as well as a special issue of Catalysis Today in 2015.⁷

Many studies have revealed the applicability of ceria as a support in different catalysts owing to its well-known high oxygen storage capacity (OSC).^1–3 It is a reducible oxide, and O vacancies can be created in ceria at elevated temperatures. This process is accompanied by the reduction of two Ce⁴⁺ cations to the Ce³⁺ state for each O vacancy and can be reversible depending on the conditions.^1–4

Doping of ceria is a commonly exploited approach for tailoring its redox properties and concomitantly affecting its catalytic behavior.⁸ The presence of aliovalent dopants⁹ causes the appearance of extrinsic defects that play a key role in the performance of many catalytic processes. Of particular importance for both technological applications and fundamental research are CeO₂ materials doped with rare-earth elements, including yttrium. When CeO₂ is doped by an oxide of a trivalent metal, such as Y₂O₃, the initial structure contains O vacancies that are created in order to compensate the effective negative charge produced by the substitution of Ce⁴⁺ with Y³⁺. Thus, a distortion is induced in the cubic matrix of the resulting Y₂O₃–CeO₂ mixed oxide. The beneficial effect of yttrium doping on the activity of ceria in soot combustion has been demonstrated by Atribak et al.¹⁰ Considering the abundance of yttrium resources, She et al. studied the effect of the addition of Y (0–5 wt%) to a CuO/CeO₂ catalyst on its activity in the water–gas shift.¹¹ The highest catalytic activity was exhibited by the catalyst doped with 2 wt% Y. The effects of different synthesis routes for the preparation of Y-doped ceria and of various amounts of Y on the performance of the preferential oxidation of CO (PROX) in the presence of excess hydrogen on supported gold catalysts have been discussed by Ilieva et al.¹² Irrespective of the preparation method, very similar degrees of conversion of CO were recorded for samples containing 1–5 wt% Y₂O₃ in the temperature range of 80–120 °C, which is of interest for the operation of fuel cells. The relatively low activity that was observed for doped ceria samples with a Y₂O₃ content of 7.5 wt% was attributed to the restricted supply of oxygen due to the assumed ordering of surface oxygen vacancies around segregated Y³⁺. Similarly, no improvement in the PROX activity and stability of gold catalysts supported on Y-doped ceria was reported by Jardim et al.,¹³ who explained this result by the surface segregation of the large amount of Y₂O₃ used (30 wt%).

Several theoretical investigations have focused on the understanding of the doping of ceria by yttrium and lanthanide metals. For instance, Dholabhai et al. studied the formation and migration of oxygen vacancies in Pr- and Gd-doped ceria with the help of methods based on density functional theory (DFT).^14,15 The local ordering of oxygen vacancies in Ce_1−xY_xO_2−x/2 where x = 0.04, 0.08, 0.11, 0.18, and 0.26 (corresponding to 2, 4, 6, 10, and 15 mol% Y₂O₃) has been investigated by analyzing data from experimental measurements (impedance spectroscopy and neutron diffraction) and ab initio molecular dynamics simulations.¹⁶ It was concluded that the ionic conductivity of Y-doped ceria decreased with an increase in the amount of Y₂O₃ because of ordering of anion vacancies. The distribution of rare-earth oxide dopants in ceria and the oxygen ion conductivity of these systems were also investigated by Grieshammer et al., who combined DFT and Monte Carlo studies.¹⁷ Hu and Metiu¹⁸ modelled the doping of the ceria(111) surface by a single Y³⁺ cation with no charge-compensating O vacancy, which would require the introduction of a second dopant cation. This single dopant lowered the formation energy of an O vacancy by more than 1.2 eV. Surprisingly, in another DFT-based study of an yttrium-doped CeO₂(110) surface¹⁹ it was reported that a peculiar arrangement of the charge distribution in the form of a structure comprising two Y³⁺ cations, one O vacancy, one Ce³⁺ cation and one O⁻ anion was found. This charge distribution was reported to be the most stable but did not provide the energetics of the conventional arrangement containing only Ce⁴⁺ and O²⁻ ions. In any case, it will be interesting for catalytic applications and further studies to understand the influence of Y³⁺ cations on the reducibility of stoichiometric CeO₂/Y₂O₃ systems and to clarify whether nanostructuring changes the behavior of these systems with respect to the surface models.

In the current study, we modeled various structures in which Y³⁺ cations replace Ce⁴⁺ cations on/in the CeO₂(111) slab and two models of ceria nanoparticles with different sizes (Ce₂₁O₄₂ and Ce₄₀O₈₀) with low Y : Ce atomic ratios of 3.9%/8.0% for the slab model and 10.5%/23.5% and 5.3%/11.1% for the two nanoparticle models, respectively. The goals of the present study were to clarify whether yttrium dopants prefer to segregate or be distantly distributed among cerium cations and to estimate changes in the reducibility of ceria in the presence of Y³⁺ cations. It is also interesting to observe the differences between the behavior of the slab, which can represent larger particles, and the nanoparticle models, where irregularities and low-coordinated surface positions (edges and corners) are more abundant.

2. Computational details and models

The calculations were performed using a periodic plane-wave DFT method with the PW91 (ref. 20) exchange-correlation functional of the generalized gradient approximation (GGA) type as implemented in the VASP program.^21–23 Stoichiometric and reduced ceria were described within the so-called GGA+U approach,^24,25 in which an on-site coulombic interaction with U = 4 eV (ref. 26) was included. Our test calculations without applying the U correction for systems that do not contain Ce³⁺ cations showed that the relative stability of the structures remains similar to that determined with U = 4 eV. A plane-wave basis set with a cut-off of 415 eV for the kinetic energy and the projector-augmented wave²³ description of core-valence electron interactions were employed. During geometry optimization all atoms were allowed to relax locally until the maximum forces that acted on each atom became less than 0.02 eV Å⁻¹. We performed only Γ-point calculations. Spin polarization was taken into account when O vacancies were modeled because this leads to the formation of Ce³⁺ cations, which have unpaired electrons. The formation energy of an oxygen vacancy E_vac was calculated with respect to half the energy of an O₂ molecule in its triplet state as follows:

E_vac = E(Ce_xO_2x−1) + 1/2 × E(O₂) − E(Ce_xO_2x)

E_vac = E(Ce_x−2Y₂O_2x−2) + 1/2 × E(O₂) − E(Ce_x−2Y₂O_2x−1)

2.1 CeO₂(111) surface slab model

The CeO₂(111) surface was modelled with a Ce₅₄O₁₀₈ slab comprising three CeO₂ (nine atomic) layers and a 3 × 6 surface cell with cell parameters of a = 13.224 Å, b = 19.836 Å, c = 19.353 Å, α = β = 90°, and γ = 60°. The lateral dimensions of the cell were fixed according to the experimental bulk geometry (a = 541 pm).²⁷ The notation used for the modelled structures is given in the sections below, where the corresponding structures are discussed.

2.2 Ce₂₁O₄₂ and Ce₄₀O₈₀ nanoparticle models

We employed a small Ce₂₁O₄₂ ceria nanoparticle with a diameter of about 1 nm,^28–34 which was located in a cubic unit cell with dimensions of 2.0 × 2.0 × 2.0 nm, as well as a larger Ce₄₀O₈₀ nanoparticle with a diameter of about 1.5 nm,^30,35–38 which was positioned in a unit cell with a rectangular parallelepiped shape and dimensions of 2.2 × 1.9 × 1.9 nm. In these models, the distance between neighboring images in the periodic unit cells is at least 0.9 nm. The notation used for the modelled structures is given in the sections below, where the corresponding structures are discussed.

3. Results and discussion

3.1 Yttrium-doped CeO₂(111) surface

For the regular CeO₂(111) surface, we modeled various structures of two types of stoichiometric system in which two and four cerium ions, respectively, were substituted by yttrium ions. In addition, for the former type of model we computed the formation energies of an oxygen vacancy.

3.1.1 Stoichiometric systems containing two yttrium ions. We modeled 27 structures in which two Y³⁺ cations replace two Ce⁴⁺ cations. In order to balance the charge on the system, one O atom was removed from the structure shown in Fig. 1 (see also Table 1). We considered structures in which an O vacancy was created in three different positions: (i) in the surface layer (SL, O_SL); (ii) in the first subsurface layer (FS, O_FS); and (iii) in the second subsurface layer (SS, O_SS). The positions of cerium cations replaced by Y³⁺ cations are specified in Fig. 1. To denote the resulting structures, the position of the created O vacancy is followed (in parentheses) by the positions of the replaced Ce⁴⁺ ions.

Fig. 1 Top (upper panel) and side (lower panel) views of CeO₂(111) surface slab model displaying one unit cell. Each position where a Ce cation is replaced by a Y³⁺ cation is denoted by a number and a capital letter. The number indicates the vertical position of the Ce cation: 1 – in the first (surface) cationic layer and 2 – in the second (subsurface) cationic layer, which is not shown in the upper panel, i.e., below an O anion in the first subsurface anionic layer. Capital letters are used to distinguish various surface and subsurface positions. An O vacancy created in the surface O layer or in the first or second subsurface O layers is marked by a green, blue or yellow circle, respectively. Color coding: Ce – yellow, O in the surface layer – dark red, O in the first subsurface layer – red, O in the second subsurface layer (located exactly below a surface Ce atom) – light red.

Table 1 Energetic (in eV), structural (in pm) and coordination characteristics of the modeled CeO₂(111) structures containing two Y³⁺ cations and one O vacancy either in the surface (O_SL), first subsurface (O_FS) or second subsurface (O_SS) layer

Structure^a	Erel^b	Y–Y^c	Y–Ovac^d	N^e	M^f
a Position of Y^3+ cations as shown in Fig. 1. Structure names in italics indicate that the Y^3+ cations are located in nearest-neighboring positions.b Relative energy with respect to the most stable structure, i.e., OFS(1B,2C).c Distance between two nearest-neighboring Y cations.d Distances between Y atoms and the O vacancy.e Coordination numbers of Y^3+ with respect to O atoms.f Number of O atoms bound simultaneously to both Y^3+ cations.
OSL(1A,1B)	0.32	421	260, 263	6, 6	1
OSL(1B,2C)	0.33	384	263, 456	6, 8	2
OSL(2C,2D)	0.40	375	458, 460	8, 8	2
OSL(2E,2F)	0.67	378	893, 894	8, 8	2
OSL(2E,2G)	0.68	376	895, 1041	8, 8	2
OSL(1H,2G)	0.79	662	890, 1042	7, 8	0
OSL(2F,1H)	0.83	381	889, 893	8, 7	2
OSL(1I,1J)	0.86	1137	596, 609	7, 7	0
OSL(1H,1J)	1.00	375	600, 887	7, 7	2
OFS(1B,2C)	0.00	410	252, 256	6, 7	1
OFS(2C,2D)	0.02	378	257, 448	7, 8	2
OFS(1A,1B)	0.26	409	250, 253	6, 6	1
OFS(2E,2F)	0.37	378	706, 801	8, 8	2
OFS(2E,2G)	0.40	378	803, 1038	8, 8	2
OFS(1H,2G)	0.51	661	799, 801	7, 8	0
OFS(2F,1H)	0.52	380	706, 800	8, 7	2
OFS(1I,1J)	0.63	1156	591, 596	7, 7	0
OFS(1H,1J)	0.73	376	595, 799	7, 7	2
OSS(2C,2D)	0.02	406	249, 252	7, 7	1
OSS(1B,2C)	0.38	420	249, 265	6, 7	1
OSS(2E,2F)	0.49	379	593, 594	8, 8	2
OSS(2E,2G)	0.55	376	595, 798	8, 8	2
OSS(2F,1H)	0.72	377	593, 701	8, 7	2
OSS(1H,2G)	0.73	660	703, 799	7, 8	0
OSS(1A,1B)	0.83	377	266, 442	7, 6	2
OSS(1I,1J)	0.91	1147	448, 797	7, 7	0
OSS(1H,1J)	0.93	381	448, 700	7, 7	2

---

The stability of the modeled structures (Table 1) depends on the positions of both the two yttrium centers and the O vacancy and varies by up to 1.0 eV. The three most stable and essentially isoenergetic structures, namely, O_FS(1B,2C), O_FS(2C,2D) and O_SS(2C,2D), contain at least one subsurface Y³⁺ cation, and the O vacancy is in either the subsurface or the second subsurface layer. The next two structures, namely, O_FS(1A,1B) and O_SL(1A,1B), in which both yttrium cations are on the surface, are less stable by ca. 0.3 eV. In all the most stable structures destabilized by up to 0.5 eV the two Y³⁺ cations are located at neighboring sites separated by up to ca. 420 pm. Note, however, that there are structures that completely or partially have the above features but are notably less stable (see the brown triangles in Fig. 2A). Fig. 2B (brown triangles), which includes all the modeled relative locations of the yttrium centers, also shows the lack of a clear relation between the relative stability and the location of the two dopant cations on the surface, in the subsurface layer or in both positions.

Fig. 2 Relations between the relative stability of the various calculated models and the nearest Y–Y distances (A) or the positions occupied by the two Y centers (B). Ceria surface – brown triangles, smaller nanoparticles – green bullets and larger nanoparticles – blue squares.

Although the O vacancies in the most stable structures are located at sites that neighbor yttrium centers, a general correlation between the Y–O_vac distances (accounted for in different ways, see Fig. S4 in ESI†) and the relative stability of all the modeled structures cannot be found.

In order to determine whether the space of the missing O atom can be refilled by an O atom, we modeled structures in which an added O atom was put in the position of the missing O atom in the Y-containing models discussed above. The systems were optimized to be spin-polarized and neutral. All the structures with an additional O atom were calculated to be less stable by 0.83, 1.27 and 1.12 eV than the initial structures, namely, O_SL(1A,1B), O_FS(1B,2C) and O_SS(2C,2D), plus half a gas-phase O₂ molecule, respectively. Hence, the filling by an oxygen atom of the charge-compensating O vacancy in such surface and subsurface positions is, not unexpectedly, energetically unfavorable up to rather high oxygen pressures.

In order to determine the relative stability of the electronic states of a system reported elsewhere,¹⁹ we performed geometry optimization of the O_SL(1I,1J) and O_SL(1B,2C) structures in the fixed triplet state. In this way, we obtained structures with an electron density distribution analogous to that in ref. 19, with one Ce³⁺ cation and partially electron-deficient oxygen anions due to the forced transfer of one O 2p electron to a Ce⁴⁺ cation. However, in the cases of both O_SL(1I,1J) and O_SL(1B,2C) this triplet state was less stable by more than 1.4 eV than the aforementioned singlet state. Hence, the spontaneous transfer of an electron from an O²⁻ ion to a Ce⁴⁺ ion in the yttrium-doped ceria systems appears to result in an excited state rather than the ground state.

3.1.2 Oxygen vacancy in systems containing two yttrium ions. In order to estimate the reducibility of ceria in the presence of two Y³⁺ cations, we modeled the creation of an oxygen vacancy in the most stable structures in which the oxygen vacancy was located in the surface layer, i.e., O_SL(1A,1B), first subsurface layer, i.e., O_FS(1B,2C), and second subsurface layer, i.e., O_SS(2C,2D). The O center was removed from various surface, subsurface and second subsurface (below surface Ce or Y) positions (see Fig. 3 and Table 2). In most cases, the locations of the two Ce⁴⁺ cations that were reduced to Ce³⁺ were calculated to be nearby the created O vacancy. We note in passing that both the nearest-neighbor and the next-nearest-neighbor locations of the two Ce³⁺ cations with respect to the O vacancies resulted from DFT-based calculations for slab models of pristine ceria surfaces.^39,40

Fig. 3 Top view of the initial structures of the doped ceria surface used to model the formation of a second oxygen vacancy: A – O_SL(1A,1B) structure, B – O_FS(1B,2C) structure, and C – O_SS(2C,2D) structure. The first O vacancy is marked as a green, blue or yellow circle; see the caption of Fig. 1. Each position of the second removed O atom is denoted by a number and a lower-case letter. The number indicates the vertical position of the removed O atom: 1 – from the first (surface) anionic layer, 2 – from the second (subsurface) anionic layer and 3 – from the third (directly below a surface cation) anionic layer. Lower-case letters are used to distinguish various surface and subsurface positions. Color coding: Ce – yellow, O in the surface layer – dark red, O in the first subsurface layer – red, O in the second subsurface layer – light red, Y – light blue. For clarity the Y–O bonds are shown.

Table 2 Energetic (in eV) and structural (in pm) characteristics of the modeled Y-doped CeO₂(111) surface structures containing two Y cations and two O vacancies

Structure^a	Erel^b	Evac^c	Ovac–Ovac^d	#Ce^e
a The positions of the second O vacancy in the most stable structures, i.e., OSL(1A,1B), OFS(1B,2C), and OSS(2C,2D) (see Fig. 1 and Table 1), are denoted by lower-case letters, as shown in Fig. 2–4, and numbers. The numbers specify the positions of the removed O atom: 1 – in the first (surface) layer, 2 – in the second (subsurface) layer, and 3 – in the third (below a surface Ce or Y cation) layer.b Relative energy with respect to the most stable structure, i.e., OFS(1B,2C)-2f.c Formation energy of an O vacancy.d Distance between two O vacancies.e Positions of the Ce^3+ cations (see numbers in Fig. S1).
OSL(1A,1B)-1a	0.41	2.61	389	39, 40
OSL(1A,1B)-2g	0.55	2.75	1179	13, 48
OSL(1A,1B)-1b	0.60	2.80	656	47, 48
OSL(1A,1B)-1c	0.60	2.80	389	15, 16
OSL(1A,1B)-3j	0.63	2.83	1014	47, 48
OSL(1A,1B)-2f	0.75	2.94	482	27, 31
OSL(1A,1B)-3h	0.90	3.10	388	16, 20
OSL(1A,1B)-2e	0.92	3.12	242	15, 16
OSL(1A,1B)-2d	1.02	3.21	233	31, 40
OSL(1A,1B)-3i	1.05	3.25	386	20, 31
OFS(1B,2C)-2f	0.00	2.52	391	27, 39
OFS(1B,2C)-1c	0.14	2.65	606	6, 16
OFS(1B,2C)-2e	0.15	2.67	388	4, 16
OFS(1B,2C)-1a	0.25	2.76	247	39, 40
OFS(1B,2C)-3j	0.29	2.80	811	47, 48
OFS(1B,2C)-2g	0.30	2.82	1009	13, 47
OFS(1B,2C)-3k	0.72	3.24	255	1, 30
OFS(1B,2C)-3i	0.76	3.28	248	4, 32
OSS(2C,2D)-2g	0.10	2.61	974	13, 48
OSS(2C,2D)-1c	0.18	2.70	662	6, 15
OSS(2C,2D)-2e	0.23	2.75	463	4, 16
OSS(2C,2D)-3i	0.38	2.90	385	4, 20
OSS(2C,2D)-3j	0.45	2.97	764	42, 48
OSS(2C,2D)-3l	0.47	2.98	386	18, 27
OSS(2C,2D)-2d	0.53	3.05	241	3, 4
OSS(2C,2D)-1a	0.53	3.05	397	39, 40
OSS(2C,2D)-2f	0.71	3.23	243	26, 27
OSS(2C,2D)-4a	0.78	3.30	253	30, 35

---

As a reference, we considered the formation of an oxygen vacancy on the surface and in the first and second subsurface layers of the undoped surface (Fig. 1). The most stable structure is that in which the O vacancy is in the first subsurface layer, for which E_vac = 2.55 eV, which is close to the value of 2.48 eV from PW91+4 calculations.⁴⁰ The structures with the vacancy in the second subsurface layer or surface layer are higher in energy by 0.22 and 0.33 eV, respectively.

The lowest formation energy of an O vacancy among the Y-doped structures corresponds to O_FS(1B,2C)-2f, in which both missing O atoms are in the first subsurface layer with E_vac = 2.52 eV, which is only slightly lower than the value calculated for the pristine CeO₂(111) surface. The next four structures have E_vac values that are higher by 0.1–0.2 eV and feature O vacancies in different locations – a combination of the surface and first and second subsurface layers. This similarity in the stability of structures with the O vacancy in different locations suggests that various spatial distributions of O vacancies may exist in Y-doped ceria. One should note the O_SL(1A,1B)-1a structure, for which E_vac = 2.61 eV, where both O vacancies are on the surface. This value is lower by 0.27 eV than the E_vac calculated for the removal of an equivalent surface O atom from the pristine CeO₂(111) surface. This result implies that the presence of Y³⁺ can facilitate the removal of surface oxygen atoms in comparison with those on the undoped ceria surface. Thus, the difference between the formation energies of an oxygen vacancy on the surface and in the subsurface layer of the oxide is reduced to only 0.10 eV from the notably larger value of 0.33 eV for the bare ceria surface.

A comparison of the locations of vacancies in different structures suggests that two O vacancies may be either close to or far from each other, but it is not favorable to remove an O atom located between two Y³⁺ cations.

In summary, doping by Y does not influence the formation energy of an O vacancy in the subsurface region. However, the formation energy of a surface vacancy is reduced in comparison with that for the undoped surface and becomes close to the formation energy of a subsurface vacancy.

This finding agrees with the results of a recent experimental study of relationships between the structural features and activity in the water-gas shift (WGS) of gold catalysts supported on Y-doped ceria (1.0, 2.5, 5.0 and 7.5 wt% Y₂O₃).⁴¹ Temperature-programmed reduction measurements revealed that the doping of ceria with Y₂O₃ affected the reduction behavior. The effect was stronger when the doped ceria supports were prepared by impregnation rather than co-precipitation. The improved reducibility of the gold catalysts on ceria with predominant surface modification (impregnated samples) was closely correlated with better WGS performance in comparison with the gold samples on supports prepared by co-precipitation with prevailing bulk modification by Y. A nanoscale 3D electron microscopy study of La-doped ceria by Collins et al. also showed that reduction occurred predominantly at the surface of nanoparticles, as estimated from the locations of Ce³⁺ ions.⁹

3.1.3 Stoichiometric systems with four yttrium ions. We obtained structures with four Y³⁺ cations and two O vacancies by adding a pair of Y³⁺ cations to the most stable structures with two Y³⁺ cations and one O vacancy, i.e., O_FS(1B,2C) and O_SS(2C,2D), and creating a new O vacancy. These added Y³⁺ cations and the new O vacancy were positioned according to the finding (see above) that pairs of Y³⁺ cations tend to stay close to each other near an oxygen vacancy. We identified five lowest-energy structures of this kind with stabilities within a range of 0.10 eV (Table 3). Some of these structures have a combination of Y ions in the surface and subsurface layers, whereas in the other structures both Y pairs are in the subsurface layer and the O vacancies are in the subsurface or second subsurface layer. In O_FS(1B,2C)/O_FS(1H,2F) and O_SS(2C,2D)/O_SS(2E,2F) the Y pairs are far apart, whereas in the other three structures (given in Table 3 in italics) the pairs are close to each other. In summary, one can see from Table 3 that structures that contain O vacancies on the surface are less stable than those with vacancies in the first and second subsurface layers. It can also be seen that the combinations of two stable configurations of a Y³⁺ pair and an O vacancy lead to stable structures with four Y³⁺ cations and two O vacancies.

Table 3 Energetic (in eV) and structural (in pm) characteristics of CeO₂(111) structures containing four Y³⁺ cations and two O vacancies

Structure^a	Erel^b	Y–Y^c	ΔY–Y^d
a Positions of Y^3+ cations as shown in Fig. 1. Structure names in italics indicate that both Y^3+ pairs are located close to each other.b Relative energy with respect to the most stable structure, i.e., OFS(1B,2C)/OFS(1H,2F).c Distances between two nearest Y cations.d Difference between the Y–Y distances in the structure with four Y^3+ cations and the structures with two Y^3+ cations, i.e., OFS(1B,2C) (410 pm) and OSS(2C,2D) (406 pm), respectively.
OFS(1B,2C)/OSL(1A,1N)	0.64	414; 421	4
OFS(1B,2C)/OFS(1A,1N)	0.19	407; 399	−3
OFS(1B,2C)/OSS(2D,2Q)	0.51	461; 402	51
OFS(1B,2C)/OSL(1A,1L)	0.80	408; 455	−2
OFS(1B,2C)/OFS(1H,2F)	0.00	408; 410	−2
OFS(1B,2C)/OFS(1R,2M)	0.10	408; 409	−2
OFS(1B,2C)/OFS(1S,2T)	0.14	410; 408	0
OSS(2C,2D)/OSS(1N,1B)	0.64	404; 417	−2
OSS(2C,2D)/OSS(1L,2K)	0.37	396; 421	−10
OSS(2C,2D)/OSS(2K,2M)	0.05	396; 407	−10
OSS(2C,2D)/OSS(2E,2F)	0.01	405; 406	−1
OSS(2C,2D)/OSS(2M,2Q)	0.10	408; 399	2

---

3.1.4 Influence of yttrium on the basicity of surface oxygen centers. We used the most stable Y-doped structures to estimate changes in the basicity of the surface oxygen centers. According to earlier studies, a shift in the O 1s core level is correlated with the basicity of the oxygen centers, as estimated from the proton affinity of these centers – a less stable O 1s core level (with a lower absolute energy E(O 1s)) corresponds to a more basic oxygen center.⁴² It was estimated that a destabilization of the O 1s core level by 1.0 eV corresponds to an increase in the proton affinity of the O center of ca. 80 kJ mol⁻¹.

Oxygen atoms in pristine CeO₂(111) feature a very distinct separation of the O 1s core levels – the surface and subsurface oxygen centers have E(O 1s) values of around −503.7 and −504.0 eV, respectively (see the green bullets in Fig. 4). The presence of yttrium cations and an oxygen vacancy modifies the O 1s energy values (see the brown triangles and blue rhombuses in Fig. 4 and S3A in ESI†), and the E(O 1s) values extend over a somewhat wider interval from −503.5 to −504.2 eV. For most of the oxygen centers, the O 1s core level is stabilized by only up to 0.1 eV (Fig. S3A in ESI†). Oxygen centers close to yttrium ions and oxygen vacancies exhibit a more specific behavior. In the structure with both yttrium cations and the oxygen vacancy in the surface layer, the O 1s level of the oxygen centers between two cations (in the first subsurface layer) is destabilized by 0.4 eV, and the other oxygen centers around each of the dopant cations that do not neighbor an oxygen vacancy are destabilized by 0.2 eV. The most basic oxygen centers in the structure are shown in indigo in Fig. S3B in ESI;† they have E(O 1s) values that are lower in magnitude by 0.2–0.4 eV than those of the surface oxygen centers in pristine ceria and are thus somewhat more basic. On the other hand, the O 1s core levels of the oxygen centers around oxygen vacancies (shown in orange in Fig. S3B†) are stabilized by 0.2 eV with respect to the corresponding centers in the pristine ceria surface.

Fig. 4 Estimated O 1s core-level energies for the pristine CeO₂(111) structure (green bullets) and the two structures with two yttrium ions, i.e., O_SL(1A,1B) (blue rhombuses) and O_SL(2C,2D) (orange triangles).

3.2 Yttrium-doped Ce₂₁O₄₂ nanoparticles

As well as the CeO₂(111) surface slab model, we modeled structures with two and four Y³⁺ cations by replacing two and four Ce⁴⁺ cations, respectively, using a smaller Ce₂₁O₄₂ nanoparticle (Fig. 5), from which one O atom was removed to maintain charge balance for each pair of Y³⁺ cations.

Fig. 5 Side and top views of Ce₂₁O₄₂ nanoparticle model. The positions of the cerium cations that are replaced by Y³⁺ are denoted by capital letters and numbers. The numbers indicate the position of the Ce cations in the surface or subsurface regions (only one Ce cation is in a subsurface position). The created O vacancy is denoted by a lower-case letter, and the number indicates whether the oxygen center is in a surface (1) or subsurface (2) position. Color coding: Ce – yellow, surface O centers – dark red, subsurface O centers – red.

There are two- and three-coordinated O centers on the nanoparticle surface, whereas the subsurface O centers are four-coordinated. The two-coordinated oxygen ions 1a, 1b, 1c and 1d shown in Fig. 5 form a small (100) facet. The corresponding Ce⁴⁺ cations, which form part of this facet, are denoted as 1A, 1B, 1C and 1D. The other two-coordinated O centers (1i, 1j, 1m, 1n, and 1o) are located at the corners of the nanoparticle. The three-coordinated O centers are denoted as 1e, 1f, and 1l, whereas the subsurface centers (four-coordinated) are denoted as 2g, 2h and 2k. Owing to the small size of the nanoparticle, there is only one subsurface Ce⁴⁺ cation, which is denoted as 2I.

3.2.1 Stoichiometric Ce₂₁O₄₂ nanoparticles containing two yttrium ions. We considered 23 structures with two Y³⁺ cations, which can be classified according to the type of oxygen center removed from the nanoparticle. The stability of the modeled structures was determined from their relative energies as shown in Table 4. In the most stable structure, i.e., O_4c(2I,1H)-2h (Fig. 6a), one of the Y³⁺ cations is in the only subsurface position, whereas another Y³⁺ cation is in the (111) surface facet and the vacancy is created by removing one of the two four-coordinated oxygen centers located between the Y³⁺ cations. In total, nine of the most stable structures have energies within a range of 0.20 eV. In these stable structures, both Y³⁺ cations replace Ce⁴⁺ cations from the top two layers, namely, the (100) facet and the layer below it, i.e., from the positions denoted as A–J; in six of the most stable structures the subsurface cation position 2I is occupied by Y³⁺. The O vacancy in these nine structures is formed by the removal of either a two-coordinated or a four-coordinated oxygen center.

Table 4 Energetic (in eV), structural (in pm) and coordination characteristics of structures of the smaller ceria nanoparticle Ce₂₁O₄₂ containing two Y³⁺ cations and one O vacancy

Structure^a	Erel^b	Y–Y^c	Y–Ovac^d	N^e	M^f
a Position of Y^3+ cations as shown in Fig. 5. Structure names in italics indicate that the Y^3+ cations are located in nearest-neighboring positions.b Relative energy with respect to the most stable structure, i.e., O4c(2I,1H)-2h.c Distance between two nearest Y cations.d Distances between Y atoms and the O vacancy.e Coordination numbers of Y^3+ with respect to O atoms.f Number of O atoms bound simultaneously to both Y^3+ cations.
O2c(1A,1B)-1b	0.31	412	215, 214	5, 5	1
O2c(1A,1E)-1a	0.39	345	215, 540	5, 5	2
O2c(1A,1H)-1a	0.16	359	215, 404	5, 6	2
O2c(1A,1C)-1a	0.15	402	215, 215	5, 5	1
O2c(1A,1C)-1m	1.59	339	771, 1005	5, 5	2
O2c(1D,1E)-1d	0.40	817	214, 778	5, 5	0
O2c(1D,1E)-1o	1.07	800	817, 553	6, 5	0
O2c(1E,1F)-1b	0.39	358	533, 373	5, 6	2
O2c(1E,1F)-1i	1.01	357	211, 411	4, 6	2
O2c(1E,1H)-1a	0.37	362	540, 404	5, 6	1
O2c(1K,1L)-1m	0.36	398	215, 215	4, 4	1
O2c(1K,1E)-1i	1.21	436	223, 211	4, 4	1
O2c(1K,1M)-1n	0.82	419	211, 216	4, 5	1
O2c(1J,1H)-1a	0.13	740	714, 404	6, 6	0
O2c(2I,1A)-1a	0.03	397	477, 215	7, 5	1
O2c(2I,1F)-1b	0.08	420	465, 373	7, 6	0
O2c(2I,1H)-1a	0.04	381	477, 404	7, 6	1
O2c(2I,1M)-1a	0.44	361	477, 690	7, 6	1
O2c(2I,1M)-1o	0.62	360	440, 215	8, 5	2
O3c(1A,1E)-1f	1.24	392	239, 228	5, 4	1
O3c(1A,1H)-1f	1.01	394	239, 228	5, 5	1
O3c(1M,1N)-1l	0.55	555	223, 223	5, 5	0
O4c(2I,1H)-2h	0.00	399	244, 228	7, 5	1
O4c(2I,1A)-2h	0.19	410	244, 230	7, 5	1
O4c(1D,1E)-2h	0.52	793	421, 433	6, 5	0
O4c(2I,1F)-2g	0.08	420	236, 221	7, 6	0

---

Fig. 6 Side (left column) and top (right column) views of optimized models of the most stable structures of the Ce₂₁O₄₂ nanoparticle with two Y³⁺ cations and one O vacancy: (a) O_4c(2I,1H)-2h; (b) O_2c(2I,1A)-1a; and (c) O_2c(2I,1H)-1a. Color coding: Ce – yellow, surface O centers – dark red, subsurface O centers – red, Y – light blue. For clarity the Y–O bonds are shown.

The calculated energies of the other 14 modeled structures are higher by 0.36–1.21 eV with respect to the most stable structure, i.e., O_4c(2I,1H)-2h. In these cases, yttrium dopants are present in all three possible cation layers of the smaller nanoparticle. In the least stable structures, of which the relative energies exceed 0.80 eV, both Y³⁺ cations are in surface positions (see the green bullets in Fig. 2B).

In general, the energetic stability of the doped structures under consideration is not related to the coordination number of the removed O center. The mutual location of two Y³⁺ cations does not determine the stability either. Hence, there is no clear correlation between the stability of a given structure and the corresponding Y–Y distance (see the green bullets in Fig. 2A).

The distances between the oxygen center that is removed and the Y³⁺ cations are not crucial for the stability of a given Ce₁₉Y₂O₄₁ structure (Table 4), which is partly due to the geometric flexibility of the nanoparticle model. For instance, in structures that are almost the most stable, namely, O_2c(2I,1A)-1a and O_2c(2I,1H)-1a, the oxygen vacancy neighbors none or only one of the Y³⁺ cations, whereas in the O_2c(1A,1B)-1b, O_2c(1K,1M)-1n, and O_2c(1K,1E)-1i structures, which are destabilized in comparison with the most stable structure, i.e., O_4c(2I,1H)-2h, by 0.31, 0.82 and 1.21 eV, respectively, the oxygen vacancy neighbors both Y³⁺ ions.

3.2.2 Oxygen vacancy in nanoparticles containing two yttrium ions. From the most stable structures with two Y³⁺ cations, i.e., O_4c(2I,1H)-2h and O_2c(2I,1A)-1a, we made models with a second O vacancy by removing one of the two-coordinated O centers (Table 5). Six models derived from O_4c(2I,1H)-2h were considered: in four of them an O center was removed from the (100) facet, whereas in the remaining models an O center was removed from an edge site. The coordination number and location of the second removed O center are specified at the beginning and end, respectively, of the structure notation. For instance, the notation O_2c(O_2c(2I,1A)-1a)-1b means that an O vacancy was created in the two-coordinated 1b position of the O_2c(2I,1A)-1a structure. The removal of an O center generates two Ce³⁺ cations. In all structures, both Ce⁴⁺ cations that were reduced were positioned directly adjacent to the site of the removed O center, except for the O_2c(O_2c(2I,1A)-1a)-1b structure, in which the second removed O atom was initially bound to surface Y³⁺ cations.

Table 5 Energetic (in eV) and structural (in pm) characteristics of the modeled structures of the smaller ceria nanoparticle Ce₂₁O₄₂ containing two Y³⁺ cations and two O vacancies

Structure^a	Erel^b	Evac^c	Ovac–Ovac^d	#Ce^e
a The letters indicate the positions of the Y^3+ cations and removed O centers shown in Fig. 5, whereas the numbers indicate whether the Y^3+ cation or removed O center is located in the first (surface) layer, 1, or the second (subsurface) layer, 2. In parentheses is the structure with two Y^3+ cations and one Ovac from which the O center is removed.b Relative energy with respect to the most stable structure, i.e., O2c(O2c(2I,1A)-1a)-1c.c Formation energy of an O vacancy.d Distance between two O vacancies.e Positions of the Ce^3+ cations (see numbers in Fig. S2).
O2c(O2c(2I,1A)-1a)-1c	0.00	1.87	459	1, 4
O2c(O2c(2I,1A)-1a)-1b	0.46	2.33	319	1, 19
O2c(O2c(2I,1A)-1a)-1d	0.42	2.29	319	4, 13
O2c(O2c(2I,1A)-1a)-1i	0.21	2.08	742	2, 10
O2c(O4c(2I,1H)-2h)-1a	0.97	2.85	262	13, 18
O2c(O4c(2I,1H)-2h)-1b	0.38	2.25	379	1, 18
O2c(O4c(2I,1H)-2h)-1c	0.16	2.03	486	1, 4
O2c(O4c(2I,1H)-2h)-1d	0.30	2.18	379	4, 13
O2c(O4c(2I,1H)-2h)-1o	0.85	2.72	527	7, 11
O2c(O4c(2I,1H)-2h)-1i	0.14	2.01	613	2, 10

---

The most stable structure is O_2c(O_2c(2I,1A)-1a)-1c, in which the two vacancies are located opposite each other in the (100) facet with E_vac = 1.87 eV. The latter value is slightly lower than 2.04 eV, which we calculated for the pristine Ce₂₁O₄₂ nanoparticle with the same locations of the O vacancy and Ce³⁺ cations. Note, however, that the lowest formation energy of an O vacancy in the Ce₂₁O₄₂ nanoparticle determined with the same computational setup, namely, 1.67 eV,^28,29 indicates that the presence of Y³⁺ does not facilitate the formation of an O vacancy. The formation energies of an O vacancy in the other positions in the (100) facet of the pristine Ce₂₁O₄₂ nanoparticle are similar, namely, 2.07 eV (1a), 2.07 eV (1b) and 2.20 eV (1d). For the removal of an O center from the 1b, 1d and 1i positions of the O_2c(2I,1A)-1a structure, the E_vac values are at least 0.2 eV higher than the value for O_2c(O_2c(2I,1A)-1a)-1c.

For the removal of an O center from the O_2c(2I,1H)-2h structure, the E_vac values vary from 2.01 eV to 2.85 eV. In a similar way to the removal of an O center from O_2c(2I,1A)-1a, the lowest E_vac values were also calculated for the removal of an O center from the 1i and 1c positions and are 2.01 and 2.03 eV, respectively. Notably, we were unable to find any correlation between the distances between two O vacancies (O_vac–O_vac in Table 5) and the E_vac values.

3.2.3 Stoichiometric Ce₂₁O₄₂ nanoparticles containing four yttrium ions. We also modeled nanoparticle structures with four Y³⁺ cations and two O vacancies. For this purpose, we used the two most stable structures with two Y³⁺ cations and one O vacancy, i.e., O_4c(2I,1H)-2h and O_2c(2I,1A)-1a, together with the O_2c(1A,1C)-1a structure, in which both Y³⁺ ions are located in the (100) facet (see Table 4). In each of these three structures we replaced two Ce⁴⁺ cations by Y³⁺ cations and concomitantly removed one O atom. The new Y³⁺ sites and O vacancy are accordingly reflected in the notation of the structures. For instance, the notation O_2c(1L,1F)-1m/O_4c(2I,1H)-2h means that in the O_4c(2I,1H)-2h structure the two Ce⁴⁺ cations in the 1L and 1F positions are replaced by Y³⁺ cations and the second O vacancy is formed in the 1m position. The calculated results for the Ce₁₇Y₄O₄₀ structures are summarized in Table 6.

Table 6 Energetic (in eV) and structural (in pm) characteristics of structures of the smaller ceria nanoparticle Ce₂₁O₄₂ containing four Y³⁺ cations and two O vacancies (Ce₁₇Y₄O₄₀)

Structure^a	Erel^b	Y–Y^c	ΔY–Y^d
a The positions of Y^3+ and O^2− ions are shown in Fig. 5. Structure names in italics indicate that both Y^3+ pairs are located close to each other.b Relative energy with respect to the most stable structure, i.e., O2c(1L,1F)-1m/O4c(2I,1H)-2h.c Distances between two nearest Y cations.d Difference between the Y–Y distances in the structures with two Y^3+ cations (O4c(2I,1H)-2h: 399 pm; O2c(2I,1A)-1a: 397 pm; and O2c(1A,1C)-1a: 402 pm) and the structure with four Y^3+ cations.
O2c(1B,1A)-1b/O4c(2I,1H)-2h	0.40	397, 408	9
O2c(1B,1D)-1c/O4c(2I,1H)-2h	0.19	400, 395	−4
O2c(1E,1F)-1b/O4c(2I,1H)-2h	0.66	356, 407	8
O2c(1K,1L)-1b/O4c(2I,1H)-2h	1.88	338, 408	9
O2c(1K,1L)-1m/O4c(2I,1H)-2h	0.23	400, 400	1
O2c(1K,1F)-1m/O4c(2I,1H)-2h	0.04	354, 399	0
O2c(1L,1F)-1m/O4c(2I,1H)-2h	0.00	354, 401	2
O2c(1L,1F)-1j/O4c(2I,1H)-2h	0.67	382, 401	2
O2c(1L,1G)-1j/O4c(2I,1H)-2h	0.86	437, 401	2
O2c(1M,1O)-1o/O4c(2I,1H)-2h	0.95	393, 396	−3
O4c(1M,1O)-2k/O4c(2I,1H)-2h	2.60	369, 432	33
O2c(1B,1J)-1b/O2c(2I,1A)-1a	0.88	348, 443	46
O2c(1B,1J)-1c/O2c(2I,1A)-1a	0.39	364, 383	−14
O2c(1G,1J)-1b/O2c(2I,1A)-1a	1.11	355, 449	52
O2c(1G,1J)-1c/O2c(2I,1A)-1a	0.47	361, 376	−21
O2c(1K,1E)-1i/O2c(2I,1A)-1a	1.01	439, 398	1
O2c(1K,1L)-1m/O2c(2I,1A)-1a	0.18	404, 398	1
O2c(1L,1G)-1j/O2c(2I,1A)-1a	0.92	440, 398	1
O2c(1L,1F)-1m/O2c(2I,1A)-1a	0.04	352, 393	−4
O2c(1G,1F)-1b/O2c(2I,1A)-1a	1.06	356, 417	20
O3c(1G,1F)-1e/O2c(2I,1A)-1a	1.26	397, 390	−7
O2c(1B,1D)-1c/O2c(1A,1C)-1a	0.38	403, 403	1
O2c(1K,1L)-1m/O2c(1A,1C)-1a	0.29	402, 403	1
O2c(1G,1L)-1j/O2c(1A,1C)-1a	1.09	442, 410	8
O2c(1G,1J)-1j/O2c(1A,1C)-1a	0.88	346, 410	8

---

There are six modeled structures with relative energies within a range of 0.23 eV, and in all of these one of the yttrium ions is located in the subsurface position denoted as 2I (see Fig. 7). The relative energy of the other studied structures ranges between 0.29 and 2.60 eV. Note that in several of these destabilized structures a Y³⁺ cation also occupies the subsurface 2I position, which hence does not seem to be a decisive factor for the stabilization of the Ce₁₇Y₄O₄₀ structures.

Fig. 7 Side (left column) and top (right column) views of optimized models of the most stable structures of the Ce₂₁O₄₂ nanoparticle with four Y³⁺ cations and two O vacancies: (a) O_2c(1L,1F)-1m/O_4c(2I,1H)-2h; (b) O_2c(1K,1F)-1m/O_4c(2I,1H)-2h; and (c) O_2c(1L,1F)-1m/O_2c(2I,1A)-1a. Color coding: Ce – yellow, surface O centers – dark red, subsurface O centers – red, Y – light blue. For clarity the Y–O bonds are shown.

The least stable structure is O_4c(1M,1O)-2k/O_4c(2I,1H)-2h, in which the four-coordinated O center between the Y³⁺ ions in the 2I and 1H cation positions is removed. As a comparison, the O_2c(1M,1O)-1o/O_4c(2I,1H)-2h structure, which has the same four Y³⁺ positions but instead a two-coordinated O center (1o) is removed, is more stable by ca. 1.7 eV than O_4c(1M,1O)-2k/O_4c(2I,1H)-2h. In the latter, four Y³⁺ cations form a triangular pyramid, inside which the 2k center is located. The distance between a Y³⁺ cation in the pyramid base and the top Y³⁺ cation is greater by 35 pm in comparison with that in O_4c(1M,1O)-1o/O_4c(2I,1H)-2h as a consequence of the removal of the 2k center.

The O_2c(1B,1D)-1c/O_2c(1A,1C)-1a structure, in which all four Y³⁺ cations are in the (100) facet and the two removed O centers were located opposite each other in the same facet, is less stable by 0.38 eV than the most stable structure, namely, O_2c(1L,1F)-1m/O_4c(2I,1H)-2h.

Table 6 also shows the differences between the Y–Y distances in the parent structures (with two Y³⁺ cations and one O vacancy) and the corresponding resulting structures with four Y³⁺ cations and two O vacancies, denoted as ΔY–Y. In some structures (including three of the most stable structures) this distance is essentially identical or only changes by up to 8 pm. However, when all the modeled structures are considered no clear correlation between the stability and the variations in the Y–Y distance can be observed.

In summary, the stability of the Y-doped structures based on the Ce₂₁O₄₂ nanoparticle is not clearly determined by only one factor, such as the coordination of the removed O centers or the Y–Y and Y–O_vac distances. In general, it can be seen that the structures in which one of the Y³⁺ cations is located in a subsurface position are among the most stable structures.

3.3 Yttrium-doped Ce₄₀O₈₀ larger nanoparticles

3.3.1 Stoichiometric nanoparticles containing two yttrium ions. The results for 26 modeled structures based on the larger Ce₄₀O₈₀ nanoparticle with different combinations of two Y³⁺ cations and one O vacancy are shown in Table 7; see Fig. 8 for details of the positions of the Y³⁺ cations and O vacancies. The most stable structure is O_2c(1E,1K)-1a, which has symmetrically situated surface Y³⁺ ions at the apexes of the nanoparticle and an O center removed from a two-coordinated position in the (100) facet. In an earlier study, the Ce⁴⁺ cations in these apex positions were found to be reduced more easily in comparison with the other Ce⁴⁺ cations in a Ce₄₀O₈₀ nanoparticle.^28,29 The next five structures that are higher in energy exhibit similar stability, with relative energies in the range of 0.27–0.32 eV. In these structures, the Y³⁺ cations are in different types of location: (i) in neighboring subsurface positions, e.g., O_2c(2H,2O)-1a; (ii) in distant subsurface positions, e.g., O_2c(2H,2S)-1a; (iii) in neighboring surface positions, e.g., O_2c(1B,1E)-1a; (iv) in distant-corner surface positions, e.g., O_2c(1L,1R)-1a; and (v) in neighboring surface and subsurface positions, e.g., O_2c(1F,2H)-1a. In all these structures, the O vacancy is in the low-coordinated surface position 1a. During geometry optimization of several structures with an O vacancy in the subsurface 2e position, the oxygen center from the 1a position moved to the 2e vacancy site. Hence, the O_2e vacancy became an O_1a vacancy.

Table 7 Energetic (in eV), structural (in pm) and coordination characteristics of structures of the larger Ce₄₀O₈₀ ceria nanoparticle containing two Y³⁺ cations and one O vacancy

Structure^a	Erel^b	Y–Y^c	Y–Ovac^d	N^e	M^f
a Positions of Y^3+ cations as shown in Fig. 8. Structure names in italics indicate that the Y^3+ cations are located in nearest-neighboring positions.b Relative energy with respect to the most stable structure, i.e., O2c(1E,1K)-1a.c Distance between two nearest Y cations.d Distances between Y cations and the O vacancy.e Coordination numbers of Y^3+ with respect to O atoms.f Number of O atoms bound simultaneously to both Y^3+ cations.
O2c(1A,1B)-1a	0.65	419	219; 214	5, 5	1
O2c(1A,2H)-1a	0.44	370	219; 446	5, 8	2
O2c(1A,1L)-1a	0.44	976	219; 900	5, 4	0
O2c(1A,2S)-1a	0.56	829	219; 945	5, 8	0
O2c(1A,1G)-1a	1.05	369	219; 567	5, 5	2
O2c(1E,1K)-1a	0.00	1034	545; 789	6, 6	0
O2c(1B,1E)-1a	0.32	354	214; 545	5, 6	2
O2c(1F,1M)-1g	0.50	383	408; 214	7, 5	2
O2c(1F,2H)-1a	0.28	375	408; 446	7, 8	2
O2c(2H,2O)-1a	0.28	372	446; 680	8, 8	2
O2c(2H,2S)-1a	0.29	540	446; 945	8, 8	0
O2c(1L,1M)-1a	0.85	731	900; 214	6, 4	0
O2c(1L,1R)-1a	0.32	1509	900; 1134	4, 4	0
O2c(1B,1D)-1a	1.87	522	214; 453	5, 6	0
O2c(1B,1D)-1p	0.86	505	1212; 1275	6, 6	0
O3c(1L,1R)-1o	1.75	1522	1437; 213	3, 4	0
O3c(1A,1G)-1f	1.92	429	242; 241	5, 5	2
O4c(1L,1R)-2e	0.32	1509	752; 983	4, 4	0
O4c(1A,1B)-2e	0.64	419	232; 232	5, 5	1
O4c(1F,2H)-2e	0.36	400	232; 230	6, 7	1
O4c(2H,2O)-2h	0.79	395	233; 233	8, 8	1
O4c(2H,2O)-2e	0.27	372	232; 376	8, 9	2
O4c(2H,2S)-2h	0.98	537	233; 446	8, 8	0
O4c(2O,2Q)-2l	0.61	372	233; 229	7, 7	1
O4c(2O,2P)-2l	0.98	537	233; 445	8, 7	0
O4c(2O,2N)-2h	0.61	395	233; 229	7, 7	1

---

Fig. 8 Side view of Ce₄₀O₈₀ nanoparticle model. The positions of the cerium cations that are replaced by Y³⁺ ions are denoted by capital letters and numbers. The numbers indicate the positions of the Ce cations in the surface or subsurface shells. The created O vacancy is denoted by a lower-case letter, and the number indicates whether the oxygen center is located in a surface (1) or a subsurface (2) position. Color coding: Ce – yellow, surface O centers – dark red, internal O centers – red.

As a general trend, structures with O vacancies created by removing three-(1o and 1f) or four-coordinated (2h and 2l) centers are unstable, and such structures become stable only when the vacancy is created in the position of a two-coordinated O center (1a). For instance: (i) the O_2c(2H,2O)-1a structure is more stable by 0.51 eV than O_2c(2H,2O)-2h; (ii) the O_2c(2H,2S)-1a structure is more stable by 0.69 eV than O_2c(2H,2S)-2h; and (iii) the O_2c(1L,1R)-1a structure is more stable by 1.43 eV than O_2c(1L,1R)-1o.

In a similar way to the abovementioned models of the ceria(111) surface and the smaller ceria nanoparticle, no clear dependence of the structural stability on the distances between two yttrium cations (blue squares in Fig. 2A) or the location of the dopant cations in surface or subsurface positions (blue squares in Fig. 2B) was found. We also examined conceivable correlations between the relative stability and other characteristics of the models – the influence of the calculated electrostatic potentials at the cerium centers substituted by yttrium and the oxygen centers that were removed (Fig. S5A–C, ESI†) and the average number of oxygen centers surrounding two Y³⁺ cations (Fig. S5D, ESI†). However, no obvious trends were found.

3.3.2 Oxygen vacancy in nanoparticles containing two yttrium ions. We modeled 17 structures (Table 8) in which an O vacancy was created in the three most stable Ce₃₈Y₂O₇₉ structures that were found, i.e., O_2c(1E,1K)-1a, O_2c(2H,2O)-1a, and O_2c(1F,2H)-1a, which feature different combinations of Y³⁺ cations: two surface, two subsurface and one surface and one subsurface, respectively. Not unexpectedly, it is most energetically favorable to create an O vacancy in a low-coordinated position in the small (100) facets. The most stable structure is O_2c((1E,1K)-1a)-1c, which has two surface Y³⁺ cations and two O vacancies at the two most distant low-coordinated sites of one (100) facet. It is less favorable to remove both O centers from neighbouring positions in a (100) facet: the corresponding structures O_2c((1E,1K)-1a)-1b and O_2c((1E,1K)-1a)-1d are less stable by 0.65–0.66 eV than the O_2c((1E,1K)-1a)-1c structure. When the second vacancy is created in another (100) facet opposite the first vacancy, the structure that is obtained, namely, O_2c((1E,1K)-1a)-1p, is less stable by only 0.08 eV than the most stable structure, i.e., O_2c((1E,1K)-1a)-1c.

Table 8 Energetic (in eV) and structural (in pm) characteristics of the modeled structures of the larger Ce₄₀O₈₀ nanoparticle containing two Y³⁺ cations and two O vacancies as well as of the pristine Ce₄₀O₈₀ nanoparticle

Structure^a	Erel^b	Evac^c	Ovac–Ovac^d	#Ce^e
a The positions of Y^3+ cations and removed O centers are shown in Fig. 8, whereas the numbers indicate whether the Y^3+ cation or removed O center is located in the surface region, 1, or the subsurface region, 2. In parentheses is the structure with two Y^3+ cations and one Ovac from which the O center is removed.b Relative energy with respect to the most stable structure, i.e., O2c((1E,1K)-1a)-1c.c Formation energy of an O vacancy.d Distance between two O vacancies.e Positions of the Ce^3+ cations (see numbers in Fig. S2).
CeO2-1a		0.76		9, 12
CeO2-1b		1.51		2, 3
CeO2-1c		1.52		3, 4
CeO2-1d		1.53		1, 4
CeO2-1k		1.60		7, 10
CeO2-1p		1.53		14, 15
CeO2-2h		1.81		5, 12
O2c((1E,1K)-1a)-1b	0.66	1.61	319	5, 3
O2c((1E,1K)-1a)-1c	0.00	0.95	449	3, 4
O2c((1E,1K)-1a)-1d	0.65	1.60	315	1, 4
O2c((1E,1K)-1a)-1k	0.40	1.35	765	7, 10
O2c((1E,1K)-1a)-1p	0.08	1.03	1304	13, 14
O2c((2H,2O)-1a)-1b	0.94	1.89	319	3, 5
O2c((2H,2O)-1a)-1c	0.36	1.31	449	3, 8
O2c((2H,2O)-1a)-1d	0.99	1.94	315	1, 8
O2c((2H,2O)-1a)-1k	0.29	1.24	765	7, 11
O2c((2H,2O)-1a)-1p	0.06	1.01	1304	9, 14
O4c((2H,2O)-1a)-2h	1.00	1.95	517	5, 6
O2c((1F,2H)-1a)-1b	0.92	1.87	319	2, 3
O2c((1F,2H)-1a)-1c	0.40	1.35	449	3, 8
O2c((1F,2H)-1a)-1d	1.02	1.97	315	1, 8
O2c((1F,2H)-1a)-1k	0.68	1.63	765	7, 10
O2c((1F,2H)-1a)-1p	0.34	1.29	1304	13, 14
O2c((1F,2H)-1a)-2e	1.50	2.45	255	2, 5

---

The lowest formation energy of an O vacancy in the latter structure is 0.9 eV, which is only slightly higher than the value calculated for the pristine Ce₄₀O₈₀ nanoparticle, namely, 0.8 eV.^28,29 The next most stable structure, namely, O_2c((2H,2O)-1a)-1p, is characterized by the fact that both Y³⁺ cations are located in neighboring subsurface positions, whereas the positions of both oxygen vacancies are the same as in O_2c((1E,1K)-1a)-1p. The former structure is less stable by only 0.06 eV than the O_2c((1E,1K)-1a)-1c structure and has a similar formation energy of an O vacancy. The other structures that were considered are less stable by at least ca. 0.3 eV than O_2c((1E,1K)-1a)-1c.

In the O_2c((1E,1K))-1a)-1p, O_2c((2H,2O)-1a)-1k, O_2c((2H,2O)-1a)-1p, O_4c((2H,2O)-1a)-2h, O_2c((1F,2H)-1a)-1p and O_2c((1F,2H)-1a)-2e structures one of the Ce³⁺ cations is located next to the created O vacancy, whereas the other Ce³⁺ cation is far from the vacancy. In all the other structures that were considered both Ce³⁺ cations are in the vicinity of the created O vacancy.

3.3.3 Stoichiometric Ce₄₀O₈₀ nanoparticles containing four yttrium ions. We considered 23 Ce₃₆Y₄O₇₈ structures containing two pairs of Y³⁺ cations and two O vacancies (see Table 9). The most stable structure is O_2c(1E,1K)-1a/O_2c(2H,2S)-1p. Other structures that are built from two stable combinations of two Y³⁺ cations and one O vacancy also result in stable arrangements, such as O_2c(1L,1R)-1a/O_2c(2H,2S)-1p, O_2c(1E,1K)-1a/O_2c(2H,2O)-1p and O_2c(2H,2O)-1a/O_2c(1L,1R)-1p.

Table 9 Energetic (in eV) and structural (in pm) characteristics of the structures of the larger Ce₄₀O₈₀ ceria nanoparticle containing four Y³⁺ cations and two O vacancies (Ce₃₆Y₄O₇₈)

Structure^a	Erel^b	Y–Y^c	ΔY–Y^d
a The positions denoted by the corresponding letters are shown in Fig. 8. Structure names in italics indicate that both Y^3+ pairs are located close to each other.b Relative energy with respect to the most stable structure, i.e., O2c(1E,1K)-1a/O2c(2H,2S)-1p.c Distances between two nearest Y cations.d Difference between the Y–Y distances in the structure with two Y^3+ cations (O4c(2I,1H)-2h: 399 pm; O2c(2I,1A)-1a: 397 pm; and O2c(1A,1C)-1a: 402 pm) and the structure with four Y^3+ cations.
O2c(1E,1K)-1a/O2c(1F,2H)-1p	0.34	1034; 369	0
O2c(1E,1K)-1a/O2c(2H,2O)-1p	0.21	1034; 371	0
O2c(1E,1K)-1a/O2c(2H,2S)-1p	0.00	1034; 531	0
O2c(1E,1K)-1a/O2c(1L,1R)-1p	0.38	1041; 1512	7
O2c(1E,1K)-1a/O2c(1B,1D)-1p	1.17	1032; 529	−2
O2c(1F,2H)-1a/O2c(1L,1R)-1p	0.36	372; 1512	−3
O2c(1F,2H)-1a/O2c(2Q,2N)-1p	0.28	372; 521	−3
O2c(1F,2H)-1a/O2c(2O,2N)-1p	0.52	367; 368	−8
O2c(1L,1R)-1a/O2c(2H,2S)-1p	0.15	1512; 531	3
O2c(2H,2O)-1a/O2c(1G,1M)-1g	0.78	373; 412	1
O2c(2H,2O)-1a/O2c(1L,1R)-1p	0.29	372; 1512	0
O2c(2H,2O)-1a/O2c(1L,1R)-1b	0.44	379; 1513	7
O2c(2H,2O)-1a/O2c(2Q,2N)-1c	0.52	374; 518	2
O2c(2H,2O)-1a/O2c(2Q,2N)-1p	0.32	368; 519	−4
O2c(2H,2O)-1a/O2c(2S,2P)-1n	0.59	370; 370	−2
O2c(2H,2O)-1a/O2c(2S,2P)-1p	0.51	373; 372	1
O2c(1A,1G)-1a/O2c(1I,1F)-1p	1.95	351; 528	−18
O2c(1A,1B)-1a/O2c(1C,1D)-1c	1.30	422; 423	3
O2c(1A,1B)-1a/O2c(1C,1D)-1p	3.23	418; 347	−1
O2c(1A,1B)-1a/O2c(1R,1L)-1c	0.50	426; 1510	7
O2c(1A,1B)-1a/O2c(1R,1L)-1d	1.24	420; 1510	1
O2c(1A,1B)-1a/O2c(1R,1L)-1p	0.80	418; 1510	−1
O2c(1A,1B)-1a/O2c(1U,1T)-1p	1.14	417; 416	−2

---

When both O centers are removed from one small (100) facet, it is not most favorable to remove them from two neighboring O sites. It is more favorable to remove two O centers that are located diagonally with respect to each other or are from two opposite (100) facets.

4. Summary and conclusions

We performed periodic density functional calculations on a CeO₂(111) slab as well as models of one smaller (Ce₂₁O₄₂) and one larger (Ce₄₀O₈₀) nanoparticle doped by two and four Y³⁺ cations. This substitution induces the formation of one and two O vacancies, respectively, for charge compensation.

4.1 CeO₂(111) slab model

The most stable structures on the CeO₂(111) surface are calculated to be those in which Y³⁺ cations are located at neighboring sites close to an oxygen vacancy created in the first or second subsurface layer.

The presence of Y³⁺ does not noticeably change the reducibility of the ceria slab model in the case of the most stable structures when the vacancy is in the first subsurface layer. However, when the oxygen vacancy is created in the surface or second subsurface layer, the formation energies of an O vacancy decrease by 0.2–0.3 eV in the presence of Y³⁺ cations in comparison with the corresponding values calculated for the pristine ceria system. This result is important from a practical point of view, because yttrium can increase the mobility of surface oxygen centers and in this way facilitate the oxidation of adsorbates with their participation.

Our CeO₂(111) slab model with four Y³⁺ cations showed that there is no clear preference regarding whether both Y³⁺ pairs stay close to each other or separate. We also found that by combining two stable configurations of a Y³⁺ pair and an O vacancy one obtains a stable structure.

When the dopant ions are in the surface layer, the basicity of surface oxygen centers around the dopant ions is slightly higher than that of the corresponding O centers on the surface of pristine ceria, as estimated from the shifts in the O 1s core level. In contrast, oxygen centers near the oxygen vacancy experience stabilization of their O 1s levels, i.e., lower basicity.

4.2 Models of smaller Ce₂₁O₄₂ nanoparticle

Our previous calculations showed that in pristine ceria nanoparticles of this kind O centers are most easily removed from low-(two-)coordinated sites.^28,29 However, we were unable to confirm this observation in the case of the doped Ce₁₉Y₂O₄₁ structures, in which one O center was removed for charge compensation upon the introduction of the two Y³⁺ cations. In the most stable structures containing yttrium dopants, oxygen centers are removed from the top two O layers: two-coordinated O centers from the (100) facet and four-coordinated O centers from the layer below. In contrast to the doped ceria slab models, in these models of ceria nanoparticle it does not seem to be obligatory that both Y³⁺ cations should be in neighboring positions near the O vacancy, which is probably due to the higher structural flexibility of the nanoparticles. Our calculations showed that in the most stable structures Y³⁺ cations replace Ce⁴⁺ cations from the (100) facet or the layers beneath it, including the only subsurface cation site.

4.3 Models of larger Ce₄₀O₈₀ nanoparticle

In the most stable structures O vacancies are created in the small (100) facets, where O centers are low-(two-)coordinated. It is not favorable to remove four-coordinated O centers from subsurface regions. In the latter case the energy of the structure is high, or a surface O center moves to the subsurface vacancy during geometry optimization and fills it. A Y³⁺ cation is stable when it replaces a Ce⁴⁺ cation in the same coordination environment; it is not necessary to remove an O center from around the Y³⁺ cation. The presence of Y³⁺ dopants can increase the formation energy of an O vacancy to an insignificant extent.

When four Y³⁺ cations are introduced into the nanoparticle it is not favorable that both O vacancies should be created in one (100) facet. In the most stable structures oxygen vacancies are located in different facets of the nanoparticle. Combinations of two stable positions of the Y³⁺ cations lead to a stable structure. There is no clear preference regarding whether Y³⁺ cations should be close to each other or well apart. This is also the case for the distance from Y³⁺ ions to O vacancies.

4.4 General conclusions

Our calculations demonstrate that the smaller is an yttrium-doped ceria system the fewer rules and trends can be elucidated regarding the energetic and structural features of the most stable locations of the dopant ions with respect to each other and to the created charge-compensating O vacancies. In the slab models one can observe some trends in the preferred locations of the Y³⁺ cations and the O vacancy. It is most stable when the Y³⁺ cations stay close to each other and the created O vacancy is in their immediate vicinity, both preferably in the subsurface region. In contrast, the smaller Ce₂₁O₄₂ nanoparticle, which has various exposed low-coordinated corner and edge sites and a very flexible structure, does not exhibit a clear tendency either for the locations of Y³⁺ cations or for the position of the O vacancy. A similar statement can be made for the results of computational modeling we obtained for the doped larger Ce₄₀O₈₀ nanoparticles. However, one trend was found in the latter case regarding the location of the O vacancies, which are preferentially formed on the small (100) facets, where less strongly bound two-coordinated O centers are available for removal.

In addition, we found that in the structures with four Y³⁺ cations and two O vacancies the most stable structures can be obtained by combining the positions in the most stable structures with two Y³⁺ cations and one O vacancy.

Hence, the main effect of yttrium doping on the ceria slab model, which is representative of large ceria particles, is the lowering of the formation energy of an oxygen vacancy for surface oxygen centers by ca. 0.3 eV, and it thus becomes almost the same as the energy required for the formation of a subsurface vacancy. In general, the doping of ceria with yttrium has a negligible effect on the reducibility of ceria surfaces and nanoparticles. One reason may be the very similar ionic radii of Y³⁺ and Ce⁴⁺ ions, which differ by only 0.03 and 0.05 Å for 6- and 8-coordinated ions, respectively.⁴³

On the other hand, the lack of preference regarding the positions of yttrium cations and oxygen vacancies in the doped ceria nanoparticles, in combination with their flexibility, would facilitate local rearrangements between structures of similar stability. This may also be related to the migration of oxygen in these systems, which can enable the dedicated use of CeO₂–Y₂O₃ nanocomposites in catalytic applications.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

Support by the Bulgarian National Science Fund (project DN-09/5), the Horizon2020 program of the European Commission (project Materials Networking, grant agreement 692146), the European Regional Development Fund and the Operational Program “Science and Education for Smart Growth” under contract UNITe no. BG05M2OP001-1.001-0004-C01 (2018–2023) is gratefully acknowledged. HAA is grateful to the Operational Program “Science and Education for Smart Growth” project BG05M2OP001-2.009-0028. KMN acknowledges support by the Spanish MINECO-FEDER grant CTQ2015-64618-R and by the grant 2017SGR13 from the Generalitat de Catalunya. The authors thank the Red Española de Supercomputación for the computer resources and the technical support.

References

1. A. Trovarelli, Catal. Rev.: Sci. Eng., 1996, 38, 439–520.
2. Catalysis by Ceria and Related Materials, ed. A. Trovarelli, Imperial College Press, London, 2002.
3. Catalysis by Ceria and Related Materials, ed. A. Trovarelli and P. Fonasiero, 2nd edn, Imperial College Press, London, 2013.
4. J. Paier, C. Penschke and J. Sauer, Chem. Rev., 2013, 113, 3949–3985.
5. D. R. Mullins, Surf. Sci. Rep., 2015, 70, 42–85.
6. T. Montini, M. Melchionna, M. Monai and P. Fornasiero, Chem. Rev., 2016, 116, 5987–6041.
7. P. Fornasiero and A. Trovarelli, Catal. Today, 2015, 253, 1–218.
8. E. W. McFarland and H. Metiu, Chem. Rev., 2013, 113, 4391–4427.
9. S. M. Collins, S. Fernandez-Garcia, J. J. Calvino and P. A. Midgley, Sci. Rep., 2017, 7, 5406.
10. I. Atribak, A. Bueno-López and A. García-García, J. Mol. Catal. A: Chem., 2009, 300, 103–110.
11. Y. She, L. Li, Y. Zhan, X. Lin, Q. Zheng and K. Wei, J. Rare Earths, 2009, 27, 411–417.
12. L. Ilieva, P. Petrova, G. Pantaleo, R. Zanella, L. F. Liotta, V. Georgiev, S. Boghosian, Z. Kaszkur, J. W. Sobczak, W. Lisowski, A. M. Venezia and T. Tabakova, Appl. Catal., B, 2016, 188, 154–168.
13. E. O. Jardim, S. Rico-Francés, F. Coloma, E. V. Ramos-Fernández, J. Silvestre-Albero and A. Sepúlveda-Escribano, Appl. Catal., A, 2014, 487, 119–129.
14. P. P. Dholabhai, J. B. Adams, P. Crozier and R. Sharma, J. Chem. Phys., 2010, 132, 094104.
15. P. P. Dholabhai, J. B. Adams, P. Crozier and R. Sharma, Phys. Chem. Chem. Phys., 2010, 12, 7904–7910.
16. M. Burbano, S. T. Norberg, S. Hull, S. G. Eriksson, D. Marrocchelli, P. A. Madden and G. W. Watson, Chem. Mater., 2012, 24, 222–229.
17. S. Grieshammer, B. O. H. Grope, J. Koettgen and M. Martin, Phys. Chem. Chem. Phys., 2014, 16, 9974–9986.
18. Z. Hu and H. Metiu, J. Phys. Chem. C, 2011, 115, 17898–17909.
19. M. Nolan, J. Phys. Chem. C, 2011, 115, 6671–6681.
20. J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh and C. Fiolhais, Phys. Rev. B: Condens. Matter Mater. Phys., 1992, 46, 6671–6687 (ibid, 1993, 48, 4978).
21. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558–561.
22. Version VASP.5.3; http://cms.mpi.univie.ac.at/vasp/.
23. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
24. V. I. Anisimov, F. Aryasetiawan and A. I. Lichtenstein, J. Phys.: Condens. Matter, 1997, 9, 767–808.
25. S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys and A. P. Sutton, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 57, 1505–1509.
26. C. Loschen, J. Carrasco, K. M. Neyman and F. Illas, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 75, 035115.
27. E. A. Kümmerle and G. Heger, J. Solid State Chem., 1999, 147, 485–500.
28. A. Migani, G. N. Vayssilov, S. T. Bromley, F. Illas and K. M. Neyman, Chem. Commun., 2010, 46, 5936–5938.
29. A. Migani, G. N. Vayssilov, S. T. Bromley, F. Illas and K. M. Neyman, J. Mater. Chem., 2010, 20, 10535–10546.
30. G. N. Vayssilov, Y. Lykhach, A. Migani, T. Staudt, G. P. Petrova, N. Tsud, T. Skála, A. Bruix, F. Illas, K. C. Prince, V. Matolín, K. M. Neyman and J. Libuda, Nat. Mater., 2011, 10, 310–315.
31. H. A. Aleksandrov, K. M. Neyman and G. N. Vayssilov, Phys. Chem. Chem. Phys., 2015, 17, 14551–14560.
32. H. A. Aleksandrov, K. M. Neyman, K. I. Hadjiivanov and G. N. Vayssilov, Phys. Chem. Chem. Phys., 2016, 18, 22108–22121.
33. I. Z. Koleva, H. A. Aleksandrov and G. N. Vayssilov, Catal. Sci. Technol., 2017, 7, 734–742.
34. M. Y. Mihaylov, E. Z. Ivanova, H. A. Aleksandrov, P. St. Petkov, G. N. Vayssilov and K. I. Hadjiivanov, Mol. Catal., 2018, 451, 114–124.
35. A. Figueroba, A. Bruix, G. Kovács and K. M. Neyman, Phys. Chem. Chem. Phys., 2017, 19, 21729–21738.
36. Y. Lykhach, A. Bruix, S. Fabris, V. Potin, I. Matolínová, V. Matolín, J. Libuda and K. M. Neyman, Catal. Sci. Technol., 2017, 7, 4315–4345.
37. A. Bruix, Y. Lykhach, I. Matolínová, A. Neitzel, T. Skála, N. Tsud, M. Vorokhta, V. Stetsovych, K. Ševčíková, J. Mysliveček, R. Fiala, M. Václavů, K. C. Prince, S. Bruyère, V. Potin, F. Illas, V. Matolín, J. Libuda and K. M. Neyman, Angew. Chem., Int. Ed., 2014, 53, 10525–10530.
38. A. Figueroba, G. Kovács, A. Bruix and K. M. Neyman, Catal. Sci. Technol., 2016, 6, 6806–6813.
39. M. V. Ganduglia-Pirovano, J. L. F. Da Silva and J. Sauer, Phys. Rev. Lett., 2009, 102, 026101.
40. S. M. Kozlov and K. M. Neyman, Phys. Chem. Chem. Phys., 2014, 16, 7823–7829.
41. T. Tabakova, L. Ilieva, I. Ivanov, M. Manzoli, R. Zanella, P. Petrova and Z. Kaszkur, J. Rare Earths, 2018.
42. G. N. Vayssilov and N. Rösch, J. Catal., 1999, 186, 423–432.
43. WebElements, https://www.webelements.com, accessed September 2018.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c8ra07014h

---