Adsorption of Water on Yttria-Stabilized Zirconia

: Water adsorption on the cubic (111) surface of yttria-stabilized zirconia (YSZ) was investigated using density functional theory calculations. Relaxation of atomic positions away from ideal cubic lattice sites, in particular of the oxygen anion sublattice, is observed on including both phase-stabilizing dopants and water adsorption. A large slab model has been used to explore the e ﬀ ects of extended relaxation throughout the anionic sublattice and the role of vacancy − vacancy interactions on water adsorption. Dissociative adsorption of water to ﬁ ll a surface vacancy site, accompanied by concerted oxygen movement in the vacancy cluster region of the slab, leads to a very strong adsorption of − 2.20 eV, blocking surface sites for oxygen activation. We show that the use of larger slab models leads to a more detailed representation of the YSZ surface system.


INTRODUCTION
Yttria-stabilized zirconia (YSZ) is an oxide material with diverse and desirable properties for a wide range of applications. It is resistant to corrosion and radiation and has high mechanical strength and chemical stability, and the presence of oxygen vacancies causes it to become a fast-ion conductor at elevated temperatures. Subsequently, YSZ has found use in thermal barrier coatings, 1 passivating media in nuclear pressure tubes, 2,3 gas sensors, 4 solid oxide fuel cell (SOFC) electrolytes, 5 and heterogeneous catalysis. 6 Of particular interest to this study is the use of YSZ in steam reforming and catalytic partial oxidation of methane (CPOM) to synthesis gas (syngas). CPOM has received increased interest as a production route for syngas due to the large abundance of natural gas and the favorable H 2 /CO ratio of CPOM compared to traditional steam reforming. 7 YSZ is commonly used as a support for metal particles for methane reforming and CPOM. However, it has been shown to be active for these reactions in its own right. Zhu et al. 8 demonstrated this behavior in a two-bed reactor setup in which CPOM occurred on YSZ starting at 550°C with complete oxygen conversion and 42% methane conversion at 950°C. Carbon dioxide and steam reforming reactions were found to take place at higher temperatures.
Aliovalent doping of zirconia with yttria, shown below in Kroger−Vink notation, serves to stabilize the cubic fluoritestructured phase while introducing charge-compensating oxygen vacancies within the material The presence of oxygen vacancies, and the fast diffusion of oxygen ions they facilitate, is thought to be crucial to the catalytic activity, as CPOM on YSZ has been shown to proceed via a Mars van Krevelen mechanism with surface vacancies providing sites for oxygen activation. 9 As one of the products of CPOM and a reactant in methane steam reforming, water is guaranteed to be present at the YSZ surface during these reactions. The water−YSZ interaction is therefore integral to the development of the atomic-scale understanding of CPOM and steam reforming while also applicable to hydrogen fuel SOFC interfaces. This fundamental water−surface interaction has been subjected to a number of theoretical studies on pure zirconia for monoclinic 10−12 and tetragonal 12−14 systems. However, study of the doped YSZ surface has so far been limited.
Theoretical surface studies of YSZ compared with those on pure ZrO 2 must deal with the additional complexity introduced into the model when doping the ZrO 2 system. Numerous experimental 15−19 (EXAFS, NMR, and neutron scattering) and computational 20−24 investigations into the nature of the yttrium−vacancy−yttrium (Y−V−Y) defect have been performed, in which the general consensus is that zirconium cations are situated in nearest neighbor (NN) sites to the oxygen vacancies, with yttrium preferentially located at next nearest neighbor (NNN) sites to the vacancies. However, the difference in energies of configurations having both yttrium ions NNN to the vacancy or having the yttrium located in a mixed NN/NNN configuration is not large. This observation is commonly explained as zirconium occupying NN sites so that it can adopt the lower 7-fold coordination normally seen in the low-temperature monoclinic zirconia polymorph or alternatively that the larger yttrium cation is better able to accommodate 8-fold coordination and so remains on NNN sites. However, the computational studies mostly focus on the introduction of a single Y−V−Y defect into a large bulk supercell or the ordered Zr 3 Y 4 O 12 δ-phase and are often not representative of real dopant concentrations used in YSZ supports. Xia et al. 23 confirmed the stability of a cubic-like phase at a 10 mol % yttria dopant concentration using interatomic potentials; meanwhile Pietrucci et al. 25 studied realistic dopant concentrations using ab initio molecular dynamics (AIMD). This AIMD study of bulk YSZ systems shows that many possible vacancy configurations are locally unstable but that multiple-vacancy concerted jumps may occur to stabilize a locally unstable arrangement. They conclude that vacancy−vacancy interactions are not fully understood but that existing models need to be improved by including these interactions. 25 Further complexity is introduced to the surface system when yttrium segregation effects are considered. Low-energy ionscattering characterization of YSZ catalysts by Zhu et al. 26 showed that the yttrium surface concentration remained at roughly 12 ± 2 mol % Y 2 O 3 , independent of the bulk yttria concentration. This level of surface yttrium saturation was predicted computationally by Xia et al. 23 and is also seen in calculations by Stanek et al. 27 However, a combined theoretical and experimental study conducted by Lahiri et al. 28 suggested yttrium surface enrichment between 30 and 45% by fitting calculated photoemission intensity ratios to high-resolution Xray photoemission spectroscopy data and a 49% surface yttrium concentration as calculated by ReaxFF simulations. The catalysts used by Zhu et al. 26 were zirconia-based powders stabilized with 5−9 wt % Y 2 O 3 and calcined in air, whereas Lahiri et al. 28 used an epi-polished single-crystal YSZ sample with no specified dopant concentration. It is probable that different preparation routes for the catalysts are responsible for the different surface yttrium concentrations observed.
The majority of computational YSZ surface studies to date have used a fixed slab (2 × 2) expansion vacuum-slab model (36 atoms, 3 trilayers). These models achieve the target dopant concentration of 9 mol % Y 2 O 3 through the introduction of just one Y−V−Y defect cluster due to the small number of atoms in the surface slab. However, these models only include the very specific vacancy−vacancy interaction case of single defect clusters repeated under periodic boundary conditions, which were employed in these calculations. Shorter-range nonperiodic vacancy−vacancy interactions are not included. Xia et al. 29 investigated oxygen adsorption on one of these fixed slab models; Cooper et al. 30 then used the same model to study methane activation; while Shishkin 31 et al. and later Gorski 32 et al. studied hydrogen adsorption on a YSZ + O surface. Only recently have larger symmetric slabs been used to study the YSZ surface, 33 although so far they have just been used to investigate the structure of a single Y−V−Y defect cluster.
In this work, we identify the modes of water adsorption on the YSZ surface at industrially relevant dopant concentrations by using large symmetric slab models: one containing a number of distinct surface regions and the other consisting of periodically repeated single defect clusters. The defect cluster region of the multiregion slab has allowed concerted vacancy− vacancy interactions upon adsorption to be modeled in this system for the first time. We show that larger slab models are needed for adsorption calculations on YSZ in order to describe the effects of extended lattice relaxation upon adsorption and to explore short-range vacancy−vacancy interactions.

COMPUTATIONAL DETAILS
Density functional theory calculations were carried out in a plane wave basis set with a 500 eV plane wave kinetic energy cut off, using the Perdew−Burke−Ernzerhof (PBE) exchangecorrelation functional 34 as implemented in the VASP code. 35,36 The core electrons are kept frozen and interact with the valence electrons in accordance with the projector augmented wave (PAW) method. 37,38 The bulk zirconia unit cell was geometry optimized using a 4 × 4 × 4 Monkhorst−Pack k-point mesh 39 giving a lattice parameter of 5.13 Å. The clean (111) surface of The Journal of Physical Chemistry C Article cubic ZrO 2 was cleaved from the optimized bulk unit cell, generating a (4 × 4) supercell expansion (288 atoms, 6 trilayers) with a slab thickness of 16.36 Å and a vacuum separation of 15 Å between the periodically repeated slabs. Two slab systems with different arrangements of defect clusters have been studied for adsorption calculations. Four defect clusters were introduced symmetrically on either surface of each slab system. The enthalpy of mixing Y 2 O 3 into the cubic ZrO 2 surface slab, for the 9 mol % dopant concentration used in these models, was calculated to be −0.08 eV per formula unit, in agreement with previous calculations. 40 Adsorption calculations on the symmetric YSZ slab were carried out with a 1 × 1 × 1 k-point mesh and were deemed converged when the forces on each ion were <0.02 eV/Å. Adsorbates were placed symmetrically on either surface of the slab, and all atoms in the simulation were relaxed. The adsorption energy (E ads ) was calculated as where E s+a is the energy of the surface−adsorbate system; E s is the energy of the bare surface; and E a is the energy of the gasphase adsorbate molecule. The factor of 2 for E a accounts for two adsorbate molecules in the simulation cell, and the resulting value of E ads is the adsorption energy per molecule. Activation barriers were calculated via transition state searches using nudged elastic band (NEB) calculations, 41 and Bader charge analysis was performed to evaluate atomic charges. 42 3. RESULTS AND DISCUSSION 3.1. Surface Models. As noted, in this work we use larger symmetric slab models to understand better the effects of surface adsorptions on YSZ. Two separate slab models with different defect cluster arrangements have been studied. The surface models used are shown in Figure 1. The multiregion slab model allows for a variety of surface topologies to be described as it contains: (i) an isolated subsurface vacancy, (ii) subsurface regions free from vacancies and (iii) a vacancy "cluster" region to probe the effect of vacancy−vacancy interactions; whereas the periodic vacancy model maintains the same dopant concentration as the multiregion model, but only contains one type of defect topology. In this case four isolated subsurface vacancies, with oxygen vacancies NNN to yttrium ions, are repeated periodically across the surface. It is important to note that for both slab systems described here only stoichiometric oxygen vacancies, resulting from charge compensation due to the aliovalent doping of Y 3+ onto Zr 4+ sites, are considered. No further surface reduction is treated within the scope of this paper.
The surface investigated in this study is the (111) terrace of cubic YSZ, widely studied in previous YSZ adsorption calculations, 32,31,29 as it has been demonstrated experimentally 43 and computationally 13,44 to be the most stable surface. The c-(111) surface was cleaved from geometry-optimized bulk ZrO 2 with a lattice parameter of 5.127 Å, which is in good agreement with previous density functional theory calculated lattice parameters of c-ZrO 2 . 44,14 Additionally, density of states simulations calculate a band gap of ∼3 eV, which agrees with previous calculated band gaps using the same exchangecorrelation functional. 45 While the lattice parameter is slightly larger and the band gap smaller than the experimental values of 5.09 Å 46 and 4−6 eV, 47,48 these respective over-and underestimations are well-known consequences of the use of GGA-DFT. The yttria concentration of the slab is 9 mol % Y 2 O 3 , and in order to account for yttrium surface segregation, four Y−V−Y defect clusters were introduced into the top two layers of the surface only, resulting in a surface dopant concentration of 14.3 mol % Y 2 O 3 . In accordance with previous findings on the nature of the Y−V−Y defect cluster, the oxygen vacancies are situated in NNN positions to yttrium cations, with the exception of the lower-most vacancy in the "vacancy cluster" region of the multiregion slab, which is NN to one yttrium and NNN to another.
In order to create the starting surfaces used for adsorption calculations, three optimization runs were performed to obtain "prerelaxed" slabs. This procedure was used to generate both multiregion and periodic vacancy prerelaxed slabs. The prerelaxation process is described in detail below for the multiregion slab; similar findings are obtained for the periodic vacancy slab. First, yttrium cations and oxygen vacancies were introduced symmetrically on either side of the slab model of the ideal c-ZrO 2 (111) surface. The surface with the vacancies introduced was geometry optimized, keeping the supercell size and shape constrained, while allowing the ionic positions to move. A large relaxation of the ions, particularly the oxygen anion sublattice, was observed, as shown in Table 1. This phenomenon has been previously reported when introducing vacancies into c-ZrO 2 models. 33 Second, a water molecule was then adsorbed symmetrically on either surface of the slab model; there is very little relaxation of the cationic positions (<0.1 Å for all atomic coordinates), yet a further large relaxation of the oxygen ion sublattice was observed (0.17 and 0.16 Å in the x and y directions). Finally, the water molecules were removed, and the slab geometry was reoptimized. The resulting structure had very little relaxation of the atomic positions (<0.055 Å for all atomic coordinates), and so this "prerelaxed" surface was then used for adsorption calculations from which the adsorption energies and geometries reported in the next section were derived. These three geometry optimization steps show that on introduction of oxygen vacancies local relaxation of atomic positions around the oxygen vacancy occurs, while atomic positions far from the The Journal of Physical Chemistry C Article vacancies remain on the ideal lattice sites. Water adsorption induces further relaxation of the oxygen sublattice, and as the sublattice remains distorted upon removal of water this result suggests that the slab generated after the first run is slightly metastable and that the adsorption of water pushes the atomic positions toward the global minimum. This prerelaxation process highlights the complexity of the energy landscape of YSZ, as the stabilized system is best described as a metastable tetragonal state in which the oxygen ions are displaced away from the ideal lattice sites of the cubic-fluorite phase. 49 The changes in ionic positions upon subsequent geometry optimizations can be seen in Figure 2. A single electronic SCF calculation for the initial atomic positions is initially performed, and a large stabilization occurs when geometry optimization of the introduced vacancies is carried out. The "prerelaxation" process of adsorbing and then removing water further stabilizes the surface slab, by 0.79 eV, due to further relaxation of the oxygen ion sublattice. It is worth noting that while this stabilization is much smaller than that of the initial geometry optimization, and relatively small when compared to the absolute total energy of the system, it is of the same order of magnitude calculated for molecular adsorptions on the surface. As discussed earlier, the current models employed for molecular adsorption calculations on the YSZ system use threelayer vacuum slabs with the bottom most layer held fixed in bulk lattice positions. A dopant concentration of 9 mol % is built into these 36 atom slabs through the introduction of a single Y−V−Y defect cluster repeated under periodic boundary conditions. However, this structure contains only the specific case of periodic vacancy−vacancy interactions, without addressing the possibility of short-range concerted vacancy interactions. An early experimental paper by Goff et al. 50 studied the defect structure of YSZ with X-ray and neutron diffraction techniques. They found vacancies to align along the ⟨111⟩ direction at yttrium concentrations below 15 mol %, and these vacancy pairs pack together at higher yttrium concentrations to resemble the Zr 3 Y 4 O 12 δ-YSZ phase. Bogicevic et al. 40 investigated the nature and strength of defect interactions in YSZ computationally and in agreement with Goff's experimental work and conclude that vacancies will order as third nearest neighbors (3NN) along the ⟨111⟩ direction, in both dilute solid solutions and the ordered 40 mol % Y 2 O 3 δ-YSZ phase, as a result of balancing electrostatic and elastic effects. Furthermore, Pietrucci et al. 25 investigated the effect of vacancy−vacancy interactions on vacancy diffusion (oxygen migration) starting with vacancies as 3NN. They found that both single vacancy and concerted multiple vacancy jumps can occur in YSZ and that, while this many-body problem is not well understood, existing models can be improved by considering vacancy−vacancy interactions. In this work, the periodic vacancy−vacancy interaction case used in earlier computational YSZ adsorption studies is modeled in the periodic vacancy slab, while the vacancy cluster region included in the multiregion surface model contains the experimentally observed 3NN oxygen vacancy separation and serves to probe the vacancy−vacancy interactions computationally as suggested in the previous studies.
The absolute energy of the periodic vacancy system is lower than that of the multiregion surface, with an energetic stabilization of −0.06 eV per formula unit of the surface, indicating that the thermodynamic ground state has a more ordered periodic state. However, the vacancies are not expected to remain in this perfect periodic arrangement under catalytic operation at high temperature. The multiregion slab, with its extended oxygen anion relaxation, therefore serves to provide insight into the nature of the surface under operating conditions.
3.2. Water Adsorption. We now explore in detail the adsorption of water both as a molecule and in a dissociative mode. Molecular water adsorption is calculated for both the periodic vacancy slab and the multiregion model, whereas full dissociative adsorption into separated surface hydroxyls was only found to be possible in the vacancy cluster region of the multiregion slab.
3.2.1. Molecular Adsorption. Water adsorption has been probed for a variety of different surface sites on the "prerelaxed" 9 mol % YSZ vacuum slab system. Geometry optimization The Journal of Physical Chemistry C Article simulations were performed placing a water molecule initially ∼1.5 Å above the surface. Water was found to adsorb by molecular chemisorption outside the defect cluster region; a typical structure is shown in Figure 3, in which the water oxygen coordinates to a surface cation and a slight elongation of one water O−H bond occurs as water hydrogen-bonds to a surface oxygen ion. This adsorption mode is the same as molecular water adsorption found on pure ZrO 2 , which has been well described previously. 14 The adsorption energies and structural properties of the molecularly adsorbed water states on the multiregion slab are given in Table 2 and are discussed in detail below.
The water molecules in configurations 1, 5, and 6 were started near the isolated subsurface vacancy; in configurations 2 and 3 they were above the cluster region and in system 4 began above a surface region with no defects. The water adsorption energies have a range of 0.5 eV, and there appears to be no single parameter on which the strength of the interaction depends. The site of initial adsorption does not provide a clear picture of the energy, as 1 and 6 are the second weakest and strongest adsorptions, yet both start near the isolated vacancy. The extent of anionic sublattice relaxation also does not fully explain the variations in the adsorption energy. The largest stabilization on adsorption occurs for system 6, in which the largest sublattice relaxation occurs. However, the second largest degree of relaxation is seen in the configuration with the lowest adsorption energy. The lack of a single parameter controlling the adsorption energy indicates that the factors controlling the strength of the interaction are complex and depend on the local structure of the adsorption site, although the relative similarities in energies indicate a flat potential energy surface on which a variety of molecularly adsorbed water species will be found.
The water adsorption modes on the periodic vacancy slab observed after structural relaxation were all molecular chemisorption configurations similar to those described above for the multiregion slab. Very little relaxation of the anionic sublattice is observed, and the water adsorption energies are ∼−0.5 and −0.8 eV when the water oxygen is coordinated to surface zirconium or yttrium, respectively. The lack of extended relaxation in the oxygen sublattice upon water adsorption and the lower adsorption energies observed on the periodic vacancy slab suggest that the energy of the interaction is dependent on the specific state of the surface due to extended oxygen sublattice relaxation, but as seen before this relationship is not straightforward.
3.2.2. Dissociative Adsorption. Dissociative chemisorption of water to two surface hydroxyls, in the defect cluster region, was found to occur with a much larger adsorption energy of −2.20 eV. The starting and final structures for the dissociative adsorption mode are shown in Figure 4. A concerted movement of oxygen ions (or vacancies) is observed on structural relaxation in which the water hydroxyl fills a surface vacancy and an interlayer oxygen ion relaxes into a subsurface vacancy site. A nudged elastic band calculation of this water dissociation shows the process to be barrierless.
Bader analysis has been performed to evaluate the charges on the ions in the water adsorption systems. For reference the   The Journal of Physical Chemistry C Article average values for ionic charges of the bare surface, gas-phase water, and the respective bulk systems are given alongside the Bader charges of the adsorbed water systems in Table 3. The average ionic charges of the atoms in the slab adsorption systems are indicated by (a) in the table, while specific individual atomic charges for oxygen ions originating from the water molecule or the YSZ surface are denoted (m*) and (s*), respectively. The positive charge for zirconium ions in the surface slab models decreases slightly compared with bulk ZrO 2 , particularly in the case of undercoordinated surface ions, and the lower values compared to ideal formal charges indicate appreciable covalent bonding character. Molecular adsorption of water leads to increased negative and positive charge on the water oxygen and hydrogen species, respectively, as electron density moves from hydrogen toward the oxygen coordinated to a surface cation. The negative charge on the surface oxygen ion H-bonded to the molecular water also increases. On dissociation to hydroxyl species, the negative charge on the water oxygen decreases, while the negative charge on the surface oxygen species, which accepts a proton from molecular water, increases. The incorporation of water oxygen into a surface vacancy as a surface hydroxyl is indicated by this reduction in negative charge. While the negative charges on the water (m*) and surface (s*) oxygen species in the dissociated water system are larger than the average oxygen ion charge of −1.19 e − in the bare surface system, they are both representative of surface oxygen species as oxygen Bader charges of up to −1.28 e − are observed in both the bulk and top layer of the bare surface system. This mode of dissociative water adsorption could not be replicated on the periodic vacancy model as no surface vacancy is available for the water hydroxyl; our model contains the most stable defect cluster arrangement of a subsurface oxygen vacancy NNN to two Y cations. When attempting to create a surface vacancy periodic defect cluster model, the surface vacancies become filled with subsurface oxygen ions during the prerelaxation process. We propose the defect cluster region of the multiregion surface, with concerted vacancy−vacancy movement, is therefore necessary for this dissociative water interaction to occur.
To date, the water−YSZ surface interaction has only been studied via investigations of hydrogen dissociation to surface hydroxyls, on a YSZ+O model and subsequent desorption of water from the surface to the gas phase. In the YSZ+O systems, the fixed slab model does not contain a charge-compensating oxygen vacancy, as it is proposed that under SOFC operation oxide ions would have diffused from the bulk to fill vacancy sites at the surface. Gorski 32 and Shishkin 31 both calculate exothermic hydrogen dissociation to surface hydroxyl species on the YSZ+O surface and find that systems of two surface hydroxyls and molecularly adsorbed water are energetically favorable compared to desorption of water to the gas phase, leaving behind a stoichiometric YSZ surface. However, the activation barriers to these processes are very different. The adsorption energies and activation energies from both these studies are given in Table 4. The adsorption energies for molecular water adsorption and dissociative chemisorption in this study are larger than have been found previously, owing to the thicker vacuum slab allowing for the modeling of extended lattice relaxation effects. However, the barrierless dissociation of adsorbed water to two surface hydroxyl species is in agreement with the work carried out by Shiskin et al. The much larger energy of the chemisorbed hydroxyls is due to the concerted oxygen movement possible in the defect cluster region of the slab, resulting from vacancy−vacancy interactions, which have not beforehand been incorporated into models of the YSZ surface system.

SUMMARY AND CONCLUSIONS
We have developed a new slab model for YSZ surface adsorption calculations, which contains distinct surface regions enabling the simulation of different surface environments within the same model. The isolated vacancy and defect-absent regions probe the common, previously studied, surface regions of YSZ, while the inclusion of a defect cluster region attempts to include experimentally and computationally observed vacancy−vacancy interactions. The much larger slab allows the inclusion of extended anionic sublattice relaxation, previously unattainable in the models used to study the YSZ surface, while the symmetric adsorption of molecules on either side of the slab negates any cross-slab dipoles introduced on adsorption. Molecular water adsorption on the surface has been well characterized, and the dissociation of surface water to surface hydroxyl species is demonstrated to be a barrierless process. The strong adsorption energies observed in our calculations for both molecular and dissociated water explain the experimentally observed strong adsorption of water on YSZ catalysts, characterized by a delayed product signal for water compared to other methane partial oxidation products, during CPOM reactions. 9 The defect cluster region of the surface model has allowed concerted oxygen (vacancy) movements, on The Journal of Physical Chemistry C Article surface adsorption, to be identified for the first time on this system. The very strong binding of a water hydroxyl into a stoichiometric surface vacancy in the defect cluster region is important when considering the catalytic activity of YSZ for CPOM and methane steam reforming. Formaldehyde has been proposed as a major intermediate for CPOM. 51 Lattice oxygen on the YSZ surface, replaced by activated molecular oxygen in a Mars van Krevelen scheme, is considered to abstract hydrogen from methane to form formaldehyde as the major intermediate. This surface formaldehyde species can then either be further oxidized by lattice oxygen to surface formate and subsequently total oxidation products or be converted directly to the partial oxidation products forming syngas. As surface vacancies are proposed as active sites for the activation of molecular oxygen, 9 strong binding of water derived hydroxyls at these sites could compete with molecular oxygen adsorption. Blocking oxygen activation could then inhibit further oxidation of formaldehyde to formate, which could possibly explain the increased selectivity to syngas observed when water is added to the CPOM reactant mixture. 52 The vacancy cluster region, where water dissociation occurs, also introduces coordinatively unsaturated zirconium ions into the surface. A related effect has been seen whereby water will compete with N 2 O for undercoordinated zirconium surface sites. 53 Further calculations will be reported to determine the surface behavior of formaldehyde and the impact that surface hydroxyls will have on formaldehyde oxidation, as well as the behavior of oxygendeficient reduced YSZ surfaces.