Gadolinium-Vacancy Clusters in the (111) Surface of Gadolinium-Doped Ceria: A Density Functional Theory Study

: Solid-oxide fuel cells are promising devices for sustainable power generation. Electrolyte materials play an important role in connecting the anode and cathode, and they in ﬂ uence the performance of the device. In this context, gadolinium-doped ceria (GDC) has proven to be an e ﬃ cient electrolyte material, although the presence of dopant clusters can lower its e ﬃ ciency. After usage, dopant clusters start appearing at dislocations, translocations, grain boundaries, or surfaces. Hence, the study of dopant clustering at the atomic level near these regions becomes of vital importance, as it allows us to understand the reasons for the occurrence of this phenomenon and its impact on the oxygen conduction. In this context, the present paper studies the impact of dopant clustering near the (111) GDC surface. We have studied two di ﬀ erent gadolinium concentrations in the material, of approximately 7% and 14%, which are close to the optimum concentration of 10%. Our results indicate that surface relaxation is a key factor in determining the preference of defect clusters to be found in the surface. We have also calculated the relative abundance of di ﬀ erent defect clusters at di ﬀ erent temperatures, including the con ﬁ gurational entropy term. It was revealed that working temperatures (650 − 1100 K) show the relative abundance of di ﬀ erent cluster structures, displaying that, at high concentrations, preferred dopant clusters resemble the structure of Gd 2 O 3 , showing the formation of gadolinia domains. Finally, we show that oxygen di ﬀ usion will be a ﬀ ected by the formation of these domains. After evaluating the oxygen mobility, we conclude that oxygen vacancies will be trapped by the


INTRODUCTION
The need for sustainable energy sources has become a priority to minimize the consequences of CO 2 production and to meet increasing energy demands. 1,2 Among the devices and materials designed to perform efficiently and under environmentally acceptable conditions, fuel cells represent an important development. In particular, solid-oxide fuel cells (SOFC) are promising devices for clean and efficient energy conversion. 3 Considering both oxygen and hydrogen as fuels, the general reactions occurring at the anode and the cathode are Anode Oxide anions generated in the cathode migrate through the electrolyte to finally react with hydrogen in the anode. Usual electrolyte materials for SOFC are ionic solids with the fluorite structure at working conditions, e.g., zirconia (ZrO 2 ) or ceria (CeO 2 ). These materials are normally doped with trivalent cations, usually lanthanides and in some cases transition metals, generating oxygen vacancies through charge compensation. 4−7 Owing to these vacancies, O 2− anions migrate by means of a hopping mechanism that has been studied extensively from a theoretical point of view. 8−11 Considering M 2 O 3 dissolved in CeO 2 , the oxygen vacancy formation reaction expressed in terms of Kroger-Vink notation 12 is where M is a trivalent metal cation. A more detailed explanation of the Kroger-Vink notation can be found in the Supporting Information.
The mobility of oxide anions is directly related to the ionic conductivity in the material, which mainly depends on two factors: electrostatic interaction between dopants and vacancies and elastic distortion generated by the dopants. 13 Several studies, either theoretical or experimental, have suggested that a 10% concentration of gadolinium oxide (Gd 2 O 3 ) is one of the best options for doping the cerium dioxide and obtaining high ionic conductivities. 14,15 In addition, gadolinium-doped ceria (GDC) allows operation at intermediate temperatures between 775 and 1075 K, 16 in comparison with other electrolyte materials like yttria-stabilized zirconia (YSZ), which operates between 873 and 1273 K. 17 Simulations of doped ionic solids often assume perfect crystals 18−20 or focus on the study of specific grain boundaries 21,22 or surfaces. 23−26 More realistic models have used molecular dynamic simulations to include dislocations and lattice misfit, although their computational cost is still very high. 27 Other types of defects to be considered originate from the segregation of dopants and associated vacancies toward the surface or the grain boundaries of the materials. 28 For example, in the case of GDC, the Gd concentration is markedly higher at surfaces than in the bulk material, 29 which has considerable impact on the ionic conductivity.
To improve our understanding of the dopant distribution in surfaces for gadolinium-doped ceria, we have used density functional theory (DFT) calculations to carry out a systematic study into Gd clustering in the dominant (111) surface, its abundance at typical SOFC working temperatures, and its impact on the mobility of O 2− near the surface.

COMPUTATIONAL METHODS AND SURFACE MODEL
We have performed ab initio calculations within the density functional theory by means of the Vienna Ab-initio Simulation Package (VASP). 30−33 Assuming the Projector-Augmented Wave (PAW) approximation, 34 we have used the Perdew−Burke−Ernzerhof (PBE) exchange correlation functional, 35 setting the kinetic energy cutoff at 550 eV. The reciprocal space was described by a 3 × 3 × 1 Monkhorst−Pack scheme for the slab model. 36 We have considered the following explicit electrons for each atomic species: 4f 1 5s 2 5p 6 5d 1 6s 2 for Ce, 4f 7 5s 2 5p 6 5d 1 6s 2 for Gd, and 2s 2 2p 4 for O. The experimental magnetization of GDC is very low, 37,38 which is why we have assumed that the total spin for the system is null, despite performing spin polarized calculations. Finally, all optimizations were performed to self-consistency, with convergence parameters of 1 × 10 −5 eV and 1 × 10 −2 eV·Å −1 for the threshold of the electronic and ionic steps, respectively. The optimization of the transition state geometries (TS) was done with the improved dimer method of Heyden et al., 39 and numerical frequencies were also calculated to confirm the presence of a single imaginary frequency for TS and none for ground state geometries. It is well-known that pure DFT fails in its description of the localized state of the 4f 1 electron in Ce 3+ . This is a consequence of the electron self-interaction error derived in the Coulomb interaction as treated in the Kohn−Sham formalism. 40 Hence, we have introduced an on-site Coulombic interaction for the 4f Ce electrons, in order to enhance their localization. The value of the U parameter was set at 5 eV, 41 under the Liechtenstein algorithm. 42 We also considered nonspherical contributions from the gradient corrections to the PAW spheres. This setup entails an expansion of the lattice parameter for CeO 2 bulk up to 5.501 Å, only 1.7% larger than the experimental one, 5.411 Å. 43 However, our goal was to model GDC, with an experimental lattice parameter of 5.423 Å, 18,44 leading to a mismatch of just 1.4%.
Our surface model was obtained by cutting the bulk material using the dipole method implemented in the METADISE code. 45,46 There are two possible terminations without perpendicular dipoles: an oxygen-terminated and cerium-terminated surface. After several tests, the most stable slab was shown to be the oxygen-terminated surface, encompassing 12 and 18 atomic layers (Ce 16 O 32 ), as shown in Figure  1. The area of the surface is 52.48 Å 2 , corresponding to a 2 × 2 (111) unit cell, with a surface energy (γ R ) of 0.72 J·m −2 , in agreement with previous results reported by several groups. 41,47 The six bottom atomic layers are fixed to model the bulk material, whereas the six top atomic layers are allowed to relax explicitly, representing the surface. These layers are labeled O a to O d and Ce a and Ce b .
Finally, sampling of all nonequivalent dopant configurations was achieved using an adaptation of the SOD code 48 to allow for double substitutions. Taking advantage of the system's symmetry, the code permits a reduction of the number of total configurations to only those that are nonequivalent by symmetry, determining at the same time their configurational degeneracy.
Ce substitutions and V O were only considered in Ce a , Ce b , and O a to O d atomic layers (See Figure 1). Since the resulting slabs were not symmetric, dipole corrections were included during optimization of the different structures using the VASP code.

RESULTS AND DISCUSSION
3.1. 6.67% Gadolinium Concentration. To model a situation where the concentration of Gd is 6.67%, i.e., slightly lower than the usual 10%, 44,49 we have substituted two cerium atoms per gadolinium and introduced an oxygen vacancy in the model depicted in Figure 1a, leading to Ce 14 Gd 2 O 31 . It is worth highlighting, however, that considering only the region where we place Gd (top two metal layers), the local concentration becomes higher, i.e., 25%, which is simulating a situation where Gd has migrated toward the surface. We first considered a structure without V O and with at least one Gd atom placed in the Ce a layer, where we identified three nonequivalent configurations, labeled as n 1 , n 2 , and n 3 (see Figure S1 for further information). Starting from these configurations, we introduced an oxygen vacancy and optimized all nonequivalent geometries for the vacancy: 8 structures for n 1 , 12 for n 2 , and 4 out of 8 for n 3 (we discuss below the reason for not optimizing all 8 geometries of n 3 ). The nomenclature that we will use from now on is n x -g y , where x refers to the Gd configuration and y to the different geometries of the V O for each x ( Figure S2).
3.1.1. Energetics, Geometry, and Electronic Structure. After optimization, we evaluated the relative energies for all the structures (collected in Tables 1 and S1). The most stable structures (n 3 -g 3 , n 2 -g 5 , n 2 -g 4 , n 3 -g 4 , n 1 -g 1 , and n 2 -g 6 , all depicted in Figure 2) are all within a range of 0.25 eV, suggesting that Gd can be placed fairly randomly in the surface, whereas the vacancy is always found in the O b atomic layer. The relative energy increases as the vacancy is located in O a , O c , and O d , respectively, reaching values up to 2 eV for n 1 -g 7 and n 1 -g 8 . The zero point energy (ZPE) for all configurations is found to be between 1.20 and 1.23 eV (Table 1), indicating that vibrational contributions are comparable for the different structures, and

Chemistry of Materials
Article the energy ordering remains unaltered once ZPE is considered. In view of the constant effect of the vibrational contributions on all configurations, we have assumed from now on that vibrations do not play an important role in our comparative study.
In a representation of the V O as a function of the relative energies, in Figure 3, it is easy to distinguish that, the deeper into the bulk from layer O b the vacancy is sited, the higher are the relative energies. This trend was first observed when analyzing all geometries for n 1 and n 2 , which therefore avoided the need to optimize all structures for n 3 , apart from those which had V O in the O a or O b layers. In fact, these results are in full agreement with previous studies of Ganduglia-Pirovano et al., Murigida and Ganduglia-Pirovano, Nolan et al., and Fabris et al., whose work has shown that, for CeO 2 (111) surfaces, oxygen vacancies are found in the oxygen sublayer. 50−53 We define the segregation energy (E seg ) as the difference between the slabs with a given defect cluster placed in the bulk and in the surface (E seg = E bulk def − E surf def ). We have used the 18 atomic layer model, depicted in Figure 1b), with the 12 topmost layers relaxed. It is worth noting that, when we mention the surface, we are referring to the six atomic layers highlighted in red in Figure 1, whereas the bulk refers to the six atomic layers highlighted in blue in Figure 1b. Using this model, we evaluated E seg for the five most stable structures, leading to positive energies in all cases, which unequivocally shows that defects are preferentially found in the surface rather than in the bulk ( Table 1). The same trend is observed at grain boundaries, as Dholabhai et al. have shown recently, with segregation energies that are slightly higher than the ones we obtain here for the surfaces. 54 This finding validates our results, as segregation is expected to be stronger at the grain boundaries.
To understand why defects are preferentially found in the surface, we have studied the distortion that occurs in the system, both in the surface and in the bulk. We have determined the relaxation energies of both the bulk and the surface (E relax,bulk and E relax,surf ), defined as the difference in energy between the optimized (E opt ) and unrelaxed (E un ) system, i.e., E relax = E opt − E un (Table 1). We have also evaluated the average distances from metal and oxygen atoms to V O (considering that it was placed in the position of the former oxygen atom), all collected in Table 2. In almost all cases studied, E relax is more negative when defects are placed in the surface, which indicates that surface relaxation is an important

Chemistry of Materials
Article factor since it enhances the accommodation of the defect clusters. For n 1 -g 1 , E rel in the surface is not as low as in the other cases, but it still shows a clear trend in terms of segregation energy. What becomes clear is that there is a thermodynamic trend for the defect clusters to be found in the surface for all systems, as indicated by the segregation energies, and the relaxation energy represents an important incentive for that segregation to occur. Compared with the ideal CeO 2 structure, Ce (as well as Gd) relaxes away from the V O (increasing the distance by ∼0.12 Å in all cases, regardless of the ion), whereas oxygen atoms move closer to V O (by between 0.20 and 0.35 Å). Bulk distortions involving neighboring oxygen atoms are only 0.20 Å. In fact, this phenomenon has been previously observed in similar systems. 9 These distortions are a direct consequence of the resulting positive charge of the oxygen vacancy, which attracts the negative charges around the point defect. The main difference between placing the defect cluster in the bulk or in the surface is that, in the latter case, the oxygens in the top oxygen layer, O a , have lower coordination numbers and this layer therefore distorts more easily. This allows a better compensation of the positive charge of the oxygen vacancy, which is reflected in the relaxation energies, since a major distortion leads to more stable final geometries.
Finally, we have analyzed the electronic structure of n 3 -g 3 (the most energetically favorable configuration). As depicted in Figure 4a, the oxygen states are below the Fermi level (0 eV). The gap between oxygen and Ce(d) states is 5.23 eV, very close to previous computational studies of ceria surfaces, 52,53 but lower than the experimental gap for CeO 2 , which is about 6 eV. 55 This disagreement between experimental and DFT+U results is, however, well documented and not relevant to the discussion of our results. 40 Density of states are present between 4.5 and 5.3 eV as a consequence of small contributions from all states, but no specific band can be assigned. The band gap between O(p) and Ce(f) is 2.20 eV. Finally, Gd(f) are found between Ce(f) and Ce(d) states, indicating no overlap between them.
3.1.2. Vacancy Formation Energy and Configurational Contribution. We have extrapolated the relative abundance of each defect structure at different temperatures. From a thermodynamic point of view, the cost of saturating an oxygen vacancy is determined via the vacancy formation energy, ruled by the stoichiometric equation (eq 1): The free energy (ΔF i v ) of this process is associated with eq 2: Surf 2 (2) where F Surf−V and F Surf are, respectively, the free energies of the doped system and that with the oxygen vacancy. Considering a constant volume, the free energy is approximated to the electronic energy plus the entropic and the temperature terms. Both entropy and enthalpy can be approximated to the product of the electronic, vibrational, and configurational contributions.
Since the electronic ground state and the vibrational contributions remain practically constant for all the different systems, we have considered their contributions as negligible and only considered the configurational entropy.
The SOD code provides us with the configuration degeneracy (Ω i ) of each nonequivalent structure. Hence, we have calculated the configurational entropy as defined by Grau-Crespo et al.: S i conf = k B ln(Ω i ), 48 and with this definition, we have deduced the vacancy formation energy (ΔF i v ) as As the energy of the oxygen molecule is overestimated by GGA, ΔF i v is therefore underestimated. 56 Several corrective terms can be found in the literature, all of them conditioned by the functional, the pseudopotential, or the consideration of Hubbard parameters in the calculation. 57 In our case, the aim of calculating these energies is not to perform a direct but a comparative analysis. Therefore, the overbinding correction becomes irrelevant.
Assuming that all different configurations are in thermodynamic equilibrium, we have assigned an occurrence probability for each configuration i, which can be understood as the configurational molar fraction (χ i ): where k B is the Boltzmann constant (8.6173 × 10 −5 eV·K −1 ). Using this equation, we extrapolated the statistical weight of each configuration between 400 and 1400 K (represented in Figure 5a). We have only highlighted those structures that show χ i higher than 0.05. At SOFC working temperatures, between 650 and 1100 K, the most stable geometry (n 3 -g 3 ) is actually not the most  abundant one, which is n 2 -g 5 . This effect is observed because configurational entropy is considered in the molar fraction analysis. It is also interesting that, at temperatures higher than 1000 K, n 3 -g 4 and n 2 -g 6 become more abundant than n 3 -g 3 .

Chemistry of Materials
Despite these variations, the molar fraction does not reveal a unique preferred site. In fact, we would suggest, from a thermodynamic point of view, that different geometries are present at a given temperature, with the proviso that Gd 3+ is found within either the first or the second metal layer, with the vacancy located in the second oxygen layer. 3.1.3. Oxygen Migration in the Surface. As mentioned in the Introduction, the activation energy (E act ) for oxygen diffusion depends on two different contributions, that can be assessed as the association energy (E ass ) and the migration energy (E mig ), where both are dependent on neighboring cations (Ce 4+ or Gd 3+ ). E ass is associated with the electrostatic interaction between the dopants and the oxygen vacancies, whereas E mig can be described as the sum of the system distortions during oxygen migration from one position to another. We should consider both when considering diffusion, but since E ass can be overcome easily at working temperatures, 58,59 we will assume that E act = E mig .
Bearing in mind how the surface influences the vacancy position and the dopant segregation, we have determined how dopant segregation can also affect the oxygen migration near the surface. Hence, we have evaluated separately all the inequivalent transition states (TS) for oxygen migration between O a and O b in the n 2 and n 3 configurations. We have not considered the n 1 geometries because of their low molar fraction at working temperatures. A further detailed explanation of the transition state geometry can be found in the Supporting Information. It is worth mentioning that, similar to our analysis of the segregation energy, here we consider O a and O d to be part of the surface region, as the accommodation of the oxygen vacancy involves oxygen atoms from both atomic layers.
All these inequivalent transition states are represented in Figure 6 for n 2 (a) and n 3 (b). Interestingly, activation energies for both sets of geometries, n 2 and n 3 , show the same trend. When an oxygen atom migrates from the second to the first oxygen layer, E act is between 0.10 and 0.27 eV for n 2 and between 0.08 and 0.23 eV for n 3 . However, migration of oxygen from the first to the second oxygen layer requires higher energies, between 0.36 and 0.59 eV for n 2 and between 0.41 and 0.49 eV for n 3 . Previously reported activation energies for bulk GDC are found to be between 0.5 and 1.0 eV. 58−61 Since activation energies are lower near the surface compared to the bulk, the vacancy diffusion to the surface is favorable. However, to drive an oxygen vacancy from the second oxygen layer to the surface (or equally, to drive an oxygen anion from the first to the second oxygen layer) is around two times higher than the

Chemistry of Materials
Article reverse process, indicating that vacancies may be trapped in the second oxygen layer, as experimental results have also suggested. This effect has also been observed near grain boundaries, where activation energies are also lower than in the bulk. 62, 63 3.2. 14.29% Gadolinium Concentration. Doubling the number of Gd centers, in a slab model, Ce 12 Gd 4 O 30 leads to a dopant concentration of 14.29% in our model (50% in surface). In this context, the number of nonequivalent configurations increases up to 476, although we have only calculated a subset of 41 structures, where we have used the following selection criterion: (i) We evaluated the relative energies of the nonreduced systems (Ce 12 Gd 4 O 32 , labeled from s 1 to s 9 ) from which we considered the six most stable structures (relative energies are listed in Table S3); (ii) we then preferentially placed the two V O in the first and the second oxygen layers, since previous results showed that they are most likely to be found in these locations; and (iii) we also tested some structures with vacancies placed deeper in the bulk for comparison.
After optimization, we found that the most stable structure is s 5 -g 5 (Figure 7b), which has both vacancies in the oxygen sublayer ( Figure S5). From that structure and in a range of 0.2 eV, we observed that structures mainly show s 3 and s 5 gadolinium distribution (See Table 3). Note that s 5 -g 3 and s 5g 2 lead to s 5 -g 7 by migration of an oxygen from the O b layer to the O a during optimization (see Figure S5).
The slab relaxation helps to accommodate vacancies in the oxygen sublayer. However, under high Gd concentration, there is a different restructuring of the surface due to the higher concentration of defects. We have compared the (111) surface of Gd 2 O 3 (crystallizing in its Ia3 group), with the six top layers involved in the Gd-Vacancy clusters. The formation of these domains has been reported previously from neutron diffraction data, 64 and comparing the resulting structure on the surface with the Gd 2 O 3 (111) surface, we have found exactly the same motif, depicted in Figure 7. This motif is observed in all structures found below 0.20 eV, showing evidence of the clustering effect at higher Gd 3+ concentrations.
The electronic structure for s 5 -g 5 (Figure 4b) shows the differences as a result of higher concentrations, where now the electronic gap is 2.10 eV between Ce(f) and O(p) states, Gd(f) orbitals are found between Ce(f) and Ce(d), and the Ce(d)− O(p) electronic gap is between 5.10 and 5.20 eV. Between 4 and 5 eV, we observe the presence of non-negligible DOS, similar to the n 3 -g 3 structure when the dopant concentration was 6.67%.
Finally, we evaluated χ i for several Ce 12 Gd 4 O 30 configurations between 400 and 1400 K, as represented in Figure 5b. As a consequence of oxygen migration occurring in s 5 -g 3 and s 5g 2 , the configurational degeneracy for s 5 -g 7 increases from 24 to 72, which finally results in this structure becoming the most abundant above 400 K, rather than s 5 -g 5 . Considering both s 5 -g 5 and s 5 -g 7 , together they represent almost 80% of the structures that could be present at fuel cell working temperatures, with the only difference between the two structures being the location of the oxygen vacancies.

CONCLUSIONS
The gadolinium-doped ceria (GDC) (111) surface has been studied using DFT+U techniques. We have analyzed at different Gd concentrations practically all nonequivalent distributions of Gd and V O at dopant concentration of 6.67%, concluding that there is no clearly preferred location for the Gd. Vacancies are preferentially sited in the oxygen sublayer at the nearest or the next nearest neighbor site from Gd, which is accommodated by surface relaxation. The segregation energies indicate that defect clusters are thermodynamically stabilized at the surface. When considering configurational entropy, the relative abundance for all the systems studied reveals that, under SOFC working temperatures, we can expect a mixture of different Gd−V O configurations. We have also observed that the mobility of oxygen vacancies is enhanced via lower activation energies, when they move from the surface to the second oxygen layer, indicating that oxygen diffusion is driven preferentially in one direction and explaining why, in the presence of Gd clusters, vacancies remain trapped near the interface. We did not observe any significant change in the geometry or the electronic structure when increasing the concentration of gadolinium dopant. At higher gadolinium concentrations, when we only considered the relative energies, we still observed a mixture of different dopant clusters. However, when we included the configurational entropy and the temperature, a single dopant cluster emerged as the most

Chemistry of Materials
Article abundant structure in the material. This finding is in line with experimental works that indicate the formation of dopant domains near the surface when the dopant concentration is high.
In conclusion, configurational contributions clearly play an important role in doped systems. The segregation of dopants and vacancies to the surface, identified in this work, may help in the design of new SOFC materials and processes.

* S Supporting Information
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.chemmater.5b02861.
A short description of the Kroger-Vink notation, oxygen migration pathway, and all the different structures optimized in this work ( Figures S3−S7), as well as a detailed list of all the relative and oxygen vacancy formation energies and configurational degeneracy for all the systems studied (Tables S1−S4) (PDF)

Author Contributions
The manuscript contains contributions from all authors. All authors have given approval to the final version of the manuscript.