Computational Modeling Study of Bulk and Surface of Yttria-Stabilized

Jul 9, 2009 - Computational Modeling Study of Bulk and Surface of Yttria-Stabilized Cubic Zirconia ... E-mail: [email protected]. ... In the Mott−Li...
0 downloads 8 Views 1021KB Size
Download PDF
3576 Chem. Mater. 2009, 21, 3576–3585 DOI:10.1021/cm900417g

Computational Modeling Study of Bulk and Surface of Yttria-Stabilized Cubic Zirconia Xin Xia,* Richard Oldman, and Richard Catlow Department of Chemistry, University College London, 20 Gordon Street, London WC1H 0AJ, United Kingdom Received February 12, 2009. Revised Manuscript Received June 20, 2009

Interatomic potential simulation techniques have been applied to both the bulk and surfaces of yttrium-stabilized cubic zirconia (YSZ). In the Mott-Littleton bulk calculations representing the infinitely dilute system, the yttrium dopants tend to exist as a pair with two yttriums close to each other and preferentially occupying the next-nearest neighbor (NNN) sites to the oxygen vacancy. However, the energy of the YSZ system as well as the configuration of defect cluster depends on the dopant concentration level. The calculated lattice energy for supercell models linearly increases with yttria content. At around 10% mol Y2O3, the simulations indicate the existence of two stable cubiclike phases differencing in term of the detached arrangement of oxygen ions. Surface-energy calculations confirm the dominance of the (111) surface in c-ZrO2 and YSZ. The yttrium solution energy at the surface has been estimated as a function of the Y dopant-vacancy cluster depth. The calculations demonstrate that for yttrium segregation to the top layers of the (111) surface. However, there is no evidence for strong segregation to the (110) surface. 1. Introduction Zirconia (ZrO2) is an important technological material that has been exploited both as a solid-state oxygen ion conductor and electrode in several key industrial applications including solid-oxide fuel cells (SOFC) and catalytic sensors. Importantly, the addition of trivalent cation dopants e.g., yttrium, substantially improves the mechanical, electrical, and chemical properties of ZrO2. The high conductivity of the doped material promotes the use of yttria-stabilized zirconia (YSZ) as a solid electrolyte; moreover, YSZ is proving to be a promising material for catalytic partial oxidation of methane (CPOM) to synthesis gas.1-3 Pure ZrO2 displays three polymorphs at atmospheric pressure: cubic (c-), tetragonal (t-), and monoclinic (m-).4-8 At ambient pressure, the ground-state ZrO2 has a monoclinic structure up to a temperature of ∼1450 K4 when the m- to t-phase transition occurs accompanied by a decrease in unit-cell volume because the 7-fold coordination of Zr changes to 8-fold.7 The cubic fluorite phase of zirconia is only thermodynamically stable at tempera*Corresponding author. E-mail: [email protected].

(1) Zhu, J.; Rahuman, M. S. M. M.; van Ommen, J. G.; Lefferts, L. Appl. Catal., A 2004, 259, 95. (2) Zhu, J.; van Ommen, J. G.; Lefferts, L. J. Catal. 2004, 225, 388. (3) Zhu, J.; van Ommen, J. G.; Knoester, A.; Lefferts, L. J. Catal. 2005, 230, 291. (4) Teufer, G. Acta Crystallogr. 1962, 15, 1187. (5) Yoshimura, M. Am. Ceram. Soc. Bull. 1998, 67, 1950. (6) Stevens, R. Zirconia and Zirconia Ceramics; Magnesium Elecktron Ltd: Manchester, ::U.K., 1986 (7) Wilson, M.; Schonberger, U.; Finnis, M. W. Phys. Rev. B. 1996-I, 54, 9147 (8) Ruff, O.; Ebert, F. Z. Anorg. Allg. Chem. 1929, 180, 19.

pubs.acs.org/cm

tures higher than 2650 K.8 In simple terms, Zr4+ appears to be too large for an efficiently packed rutile structure and too small to form the fluorite structure. The tetragonal phase can be consider as a distorted cubic structure with displacement of oxygen ions along the Æ001æ direction (c axis).9 Phase stabilization of ZrO2 is important in high-temperature applications where cycling occurs through the m-t transition temperature during material preparation or use. For example, in catalytic steam reforming of methane, the substantial change in unit-cell volume leads to mechanical failure of catalyst pellets. Stabilization of t- or c-ZrO2 down to room temperature can be achieved in a number of ways but most importantly by doping with Y3+.10,11 Doping to create oxygen vacancies according to the defect reaction below, expressed in the Kroger-Vink notation, reduces the average coordination number and stabilizes the cubic phase 0 33 Y2 O3 þ 2ZrX Zr0 þ 1=2O2 f 2YZr þ VO þ 2ZrO2 ð1Þ

In this context, we can also note that the D53 sesquioxide structure of Y2O3 has a similar oxygen anion sublattice to fluorite ZrO2.12 The oxygen vacancies in YSZ, in addition to stabilizing ZrO2 in its cubic form, are thought to play an important role in high temperature catalytic processes (9) Ho, S. M. Master. Sci. Eng. 1982, 54, 23. (10) Yashima, M.; Ohtake, K.; Kakihana, M.; Arashi, H.; Yoshimura, M. J. J. Am. Ceram. Soc. 1996, 57, 17. (11) Yashima, M.; Takanhashi, H.; Ohtake, K.; Hirose, T.; Kakihana, M.; Arashi, H.; Yoshimura, M. J. J. Am. Ceram. Soc. 1996, 57, 289. (12) Wells, A. F. Structural Inorganic Chemistry, 4th ed.; Claredon Press: Oxford, U.K., 1975

Published on Web 07/09/2009

r 2009 American Chemical Society

Article

such as direct partial oxidation of methane or steam reforming. In particular, Zhu et al. have studied catalytic partial oxidation of methane over pure ZrO2 (monoclinic) and YSZ (cubic/tetragonal).13 Their experiments have shown that YSZ has superior activity compared to pure ZrO2. From a combination of 18O isotope exchange and pulsed CH4 and O2 studies they have concluded that CH4 is selectively oxidized at surface lattice oxygen vacancies, where oxygen is rapidly replenished by lattice diffusion because oxygen mobility is enhanced by oxygen vacancies, leading to higher activity. The results also indicate that the oxygen vacancies in YSZ are important in the reduction of carbon deposition, which otherwise limits catalytic lifetime. The importance of oxygen vacancies has also been discussed in relation to steam reforming over Ni supported on a variety of perovskites, e.g., SrTiO3, by Urasak et al.14 and CO2 reforming of CH4 over Ni on Sm- and Gd-doped CeO2 by Huang et al.15 There is extensive experimental data on bulk Y-doped ZrO2. The phase diagram for YSZ has been well-described by Yashima et al.16 in terms of two metastable tetragonal phase, t0 and t00 , in addition to the three equilibrium phases m, t, and c discussed above. The formation of m-ZrO2 is confined to a small region of the phase diagram below 3 mol % Y2O3. At concentrations less than 10 mol %, the system contains the two metastable tetragonal phases. Both of these phases belong to the P42/nmc space group and do not transform directly to the equilibrium phases. The t0 form (c/a > 1) exists up to 8 mol % Y2O3 and is frequently described as partially stabilized zirconia (PSZ) because it phase separates as a function of temperature and time to give a mixture of mand fully stabilized c-ZrO2 (FSZ).17 The t00 form, described as cubiclike (c/a ≈ 1) but detected by laser Raman spectroscopy as well as electron and neutron diffraction,11 is observed between 8 and 10 mol %. This phase is often described as fully stabilized, but its stability remains to be fully tested under aggressive conditions. Beyond 10 mol %, Y2O3 zirconia is fully stabilized in its cubic form down to room temperature. According to Pascual and Dur an, for yttria content higher than 40 mol %, the system crystallizes as Y4Zr3O12. Mixtures of the cubic Y-doped ZrO2 solid solution and segregated microdomains of Y4Zr3O12 have been proposed for yttria content between 15 and 40 mol %.18 The surface structure of Y doped ZrO2 is less well understood than the bulk, although a number of experimental studies have indicated a degree of Y surface segregation, for example, from X-ray photoelectron spec(13) Zhu, J.; van Ommen, J. G.; Henny, J.; Bouwmeester, M.; Lefferts, L. J. Catal. 2005, 233, 434. (14) Urasaki Appl. Catal., A 2005, 286, 23. (15) Huang, T. J.; Lin, H. J.; Yu, T. C. Catal. Lett. 2005, 105, 239. (16) Yashima, M.; Morimoto, K.; Ishizawa, N.; Yoshimura, M. J. J. Am. Ceram. Soc. 1993, 76, 1745. (17) von Schusterius, C.; Padurow, N. N. Berlin Deutsch. Keram. Gesell (DKG) 1953, 30, 235. (18) Pascual, C.; Duran, P. J. Am. Ceram. Soc. 1983, 66, 22. (19) Bernasik, A.; Kowalski, K.; Sadowski, A. J. Phys. Chem. Solids 2002, 63, 233. (20) Morterra, C.; Cerrato, G.; Ferroni, L.; Montanara, L. Mater. Chem. Phys. 1994, 37, 243.

Chem. Mater., Vol. 21, No. 15, 2009

3577

troscopy (XPS)19,20, Auger electron spectroscopy (AES),21 and low energy ion spectroscopy (LEIS).22,23 Morterra et al. from XPS measurements have proposed that Y segregation is favored under conditions of marginal stabilization, i.e.. ∼2-3 mol % Y2O3.20 A detailed surface model resulting from LEIS measurement for cZrO2 has suggested that there is an yttrium-rich monolayer on top of a thicker (4 nm) undoped subsurface layer. The authors also proposed that segregation to the exposed surface is more extensive than to individual grain boundaries.22,23 Not only is there uncertainty about the dispersion of Y, but also about the location of compensating oxygen vacancies. There have also been several computational studies of the surface structure of zirconia. For pure cubic ZrO2, interatomic potential.24 and quantum mechanical (QM) Hartree-Fock calculations25 have shown that the most stable surface for c-ZrO2 is (111)c, followed by (011)c. Similar results are obtained by density functional theory (DFT) calculations26 on the t-phase where the most stable (101)t surface is closely equivalent to (111)c. Christensen and Carter have also confirmed that the stoichiometric oxygen (O) terminated surface is considerably more stable than the Zr or O-O terminations independent of the surface orientation.26 For Y-doped t- ZrO2, Eichler et al. have demonstrated that the (101)t surface remains the most stable and that there is an energetic driving force for Y segregation. In contrast, for the less-stable (001)t face, equivalent to (110)c, Y is not surface segregated.27 These conclusions are supported by other DFT calculations of Ballabio et al.28 Almost all the above studies have concerned the dispersion of Y. However, for a better understanding of catalyst behavior, it is more important to understand the dispersion of the compensating oxygen vacancies in relation to defining the nature of the active sites in high temperature oxidation catalysis. To this end, we have carried out atomistic calculations to study Y-doped cZrO2 in more detail, to elucidate the formation, location, and orientation of the compensating defect clusters as a function of Y concentration. We have chosen to employ the interatomic potential-based methods, as they allow us to exploit a wide range of options in terms of concentration and dispersion of the defects. Besides the surface studies, we have performed bulk calculation as a function of Y concentration with particular reference to the phase behavior of 8-10 mol % Y-doped ZrO2, as this region is (21) Winnubst, A. J. A.; Kroot, P. J. M.; Burggraaf, A. J. J. Phys. Chem. Solids 1983, 44, 955. (22) Ridder, M. de; Vervoort, A. G.; van Welzenis, R. G.; Brongersma, H. H. Solid State Ionics 2003, 156, 255. (23) Ridder, M. de; van Welzenis, R. G.; Brongersma, H. H.; Kreissig, U. Solid State Ionics 2003, 158, 67. (24) Balducci, G.; Kaspar, J.; Fornasiero, P.; Graziani, M.; Islam, M. S. J. Phys. Chem. B 1998, 102, 557. (25) Gernnard, S.; Cora, F.; Catlow, C. R. A. J. Phys. Chem. B 1999, 103, 10158. (26) Christensen, A.; Carter, E. A. Phys. Rev. B 1998, 58, 8050. (27) Eichler, A; Kresse, G. Phys. Rev. B 2004, 69, 045402. (28) Ballabio, G.; Bernasconi, M.; Pietrucci, F.; Serra, S. Phys. Rev. B 2004, 70, 075417.

3578

Chem. Mater., Vol. 21, No. 15, 2009

Xia et al.

Table 1. Potential Parameters Employed in the Atomistic Calculations32 (i) Short Range interaction

A (eV)

F (A˚)

C (eV A˚6)

cut-off (A˚)

Zr4+ 3 3 3 O2O2- 3 3 3 O2Y3+ 3 3 3 O2-

985.87 22764.00 1345.10

0.3760 0.1490 0.3491

0.0 27.88 0.0

10 12 10

(ii) Shell Model species

k/eV A˚ -2

Y/e -2.077 1.35

O2Zr4+

27.290 169.617

of considerable technological interest because of its superior mechanical properties in relation to catalytic pellet performance.3 2. Computational Methods The bulk and surfaces of cubic zirconia systems were modeled using standard lattice simulation techniques. The bulk simulations that used the GULP code (General Utility Lattice Program) are formulated within the framework of the Born model of ionic solids in which the dominant longrange interactions are Coulombic. Buckingham potentials eq 2 describe the combination of the short-range repulsion between neighboring electron clouds with Van de Waals attraction.29 Eshort range

    -r C - 6 ¼ Aexp F r

ð2Þ

Ionic polarization is taken into account using the shell model,30 which represents each ion as a massless shell with charge Y connected to an inner core through a harmonic potential of the form 1 Ecore -shell ¼ k2 d 2 2

ð3Þ

where k is the spring constant and d is the relative displacement of the core and shell. The total charge of the shell charge Y plus core system equals the integral charge of the ionic species. The relevant parameters employed in eqs 2 and 3 are taken from earlier work by Dwivedi and Cormack31 and tabulated in Table 1. The defect structures and energies have been modeled using both the supercell and the Mott-Littleton (ML) approaches. The supercell calculations, which apply 3D periodic boundary conditions, are used for systems with a higher concentration of defects, whereas in the ML approach, the isolated cluster is embedded in the crystal with a surrounding region of ions relaxed to zero force. Unlike the periodic techniques it can be applied to the charged system as well as the systems that are charged neutral with compensating. In both approaches, a cubic fluorite structure was assumed for the initial undoped structure, becausee the high symmetry of the cubic relative to the monoclinic phase makes the system less computationally demanding. In addition, the most interesting experimental results have been obtained with materials with cubic structure.3 The supercell calculations investigated Y3+ cation doping and the creation of (29) Catlow, C. R. A.; Kotomin, E. Computational Materials Science; IOS Press: Amsterdam, 2001. (30) Dick, B. G.; Overhauser, A. W. Phys. Rev. 1958, 112, 90. (31) Dwivedi, A.; Cormack, A. N. Philos. Mag., A 1990, 61, 1.

oxygen anion vacancies in bulk ZrO2 by locating the compensating yttrium defects either in nearest-neighbor (NN) sites or next-nearest-neighbor (NNN) sites to the oxygen vacancy. The simulations were carried out as a function of Y concentration by varying the supercell size from 111 to 422. Lattice energies were calculated using the standard procedure available in the GULP program and employing the same potentials as used in the rest of the study; the values reported are obtained after full minimization of cell dimensions and unit-cell coordinates. Phonon dispersion curves have also been examined for the relaxed structures obtained from the GULP calculations because experimental studies have suggested there are displacive phase transitions of ZrO2. To deal with the imaginary modes indicated by negative phonon frequencies, the Rational Function Optimiser (RFO) has been employed.33 Born-Haber cycles, derived from eq 1 and incorporating the lattice energies of ZrO2 and Y2O3, were employed to calculate the solution energy of Y2O3 in ZrO2. Calculations representing the dilute case employing the Mott-Littleton approach are based on explicitly relaxing ions within an inner spherical region I surrounding the defect with polarization of the outer region II being calculated via a quasi continuum approximation.34 The sizes of region I (10 A˚) and region II (20 A˚) were determined from the convergence of the calculated energy, i.e., the lattice energy changes less than 0.01 eV on expansion of either. The binding energy of defect clusters has been estimated from eq 4 Ebind ¼ Ecluster -f

X

Eisolated defect g

ð4Þ

component

where Ecluster is the energy of a charge neutral cluster of compensating defects and Eisolated defect is the energy of an individual defect, i.e., Y3+ or O2- vacancy. The surface studies of ZrO2 and YSZ systems were performed within a 2-dimensional slab model applying periodic boundary conditions using the MARVIN code.35 The slab model, which consists of one or more blocks, is split into two regions (I and II): region I includes the surface layer in which all ions are relaxed explicitly to zero net force, and region II includes those atoms more distant form the surface, which are kept fixed at their bulk equilibrium positions. The oxygen (O)-terminated surfaces have been chosen for the calculations, because they were identified to be more stable than Zr or O-O termination from the first principle studies by Eichler et al.27 The (100) surface was (32) Khan, M. S.; Islam, M. S.; Bates, D. R. J. Mater. Chem. 1998, 8(10), 2299. (33) Banerjee, A.; Adams, N.; Simons, J.; Shepard, R. J. Phys. Chem. 1985, 89, 52. (34) Mott, N.; Littleton, M. Trans. Faraday Soc. 1938, 34, 485. (35) Gay, D. H.; Rohl, A. L. J. Chem. Soc, Fraday Trans. 1995, 91(5), 925.

Article

Chem. Mater., Vol. 21, No. 15, 2009

3579

Table 2. Calculated and Observed Structure Parameters of Cubic ZrO2 and Tetragonal ZrO2 structure parameters Zr 3 3 3 OI Zr 3 3 3 OII OI 3 3 3 OII Zr 3 3 3 Zr a c c/a dz R = β = γ (deg) lattice energy per ZrO2 (eV) a

r A˚ (expt)

cubic ˚ r A (calcd)

2.204 2.204a 2.563a 3.589b 5.070c

2.198 2.198 2.538 3.589 5.076

1.000

1.000

90 -109.76

90 -109.76

ΔA˚

Δ%

-0.006 -0.006 -0.025

-0.27 -0.27 -0.97

0.006

0.11

r A˚ (expt)

tetragonal ˚ r A (calcd)

2.065 2.463 2.635 3.683 3.640d 5.270d 1.448d 0.065e 90 -109.74

2.049 2.416 2.663 3.639 3.588 5.217 1.454 0.060 90 -109.78

ΔA˚

Δ%

-0.016 -0.047 0.028 -0.044 -0.052 -0.053 0.006 0.005

-0.77 -1.91 1.06 -1.19 -1.43 -1.01

From ref 31. b From ref 39. c From ref 40. d From ref 41. e From ref 42.

Figure 1. Schematic diagrams for the NN and NNN sites of yttrium to oxygen vacancy (Zr, blue-green; O, red; O vacancy, white; alternative sites for Y, yellow).

reconstructed to quench the dipole perpendicular to the surface plane in the manner proposed by Tasker36 before energy minimization is performed. To achieve insensitivity to region size defined in the terms of a surface energy change of less than 0.01 J/m2 on expansion of either, we determined the sizes for region I and II required for simulation of the (111) and (110) surfaces to be 6 and 12 unit cells; for the reconstructed (100) surface, 6 and 9, and for the (310) surface, 6 and 20 unit cells, respectively. All the structures were then relaxed using the Broyden-FletcherGoldfarb-Shanno algorithm (BFGS).37 The surface energy is defined as the excess energy at the surface of a material compared to the bulk according to γ ¼

Esurf - Ebulk A

ð5Þ

When modeling the pure surface, this excess energy γ is evaluated from the energy difference for equal numbers of ions in the surface and bulk blocks per unit cell. A is the surface area. However, in the case of the Y-doped surface, oxygen vacancies are created. Because the values of the bulk energy automatically used in calculations employing the MARVIN code are the energies of the perfect lattice in region II of the model, a correction term is required to obtain the surface energy. For this purpose, the difference in bulk energy for doped and pure cZrO2 generated in GULP is applied to correct the bulk energy (36) Tasker, P. W.; Colbourn, E. A.; Mackrodt, W. C. J. Am. Ceram. Soc. 1985, 68, 74. (37) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T. Flannery, B. P., Numerical Recipes, 2nd ed.; Cambridge University Press: Cambridge, U.K., 1992.

from MARVIN for the presence of defects. The details have been discussed in our previous studies,38 and this correction method has been adopted throughout the present work.

3. Results and Discussion 3.1. Bulk Structures. The validity of the atomistic potential model was first assessed by simulating the three low-pressure phases of zirconia. The cubic and tetragonal structures of ZrO2 were well-reproduced by the GULP code. The lattice parameters obtained (Table 2) show excellent agreement with the experimental results and first principle quantum mechanical calculations.31,39-42 However, our attempt to model the monoclinic phase using the GULP code predicts an energetically more favorable orthorhombic structure, which has also been suggested by the classical potential calculations of Stefanovich, Shluger, and Catlow.43 The structure of yttria, which has a sesquioxide cubic lattice, has also been modeled. The lattice parameters obtained are within 1% of the experimental data. The lattice energy per Y2O3 unit for the equilibrated structure (38) Khan, S.; Oldman, R. J.; Cora, F.; Catlow, C. R. A.; French, S. A.; Axon, S. A. Phys. Chem. Chem. Phys. 2006, 44, 5207. (39) Ruhle, M.; Heuher, H. Adv. Ceram. 1984, 12, 14. (40) Catlow, C. R. A. J. Chem. Soc., Faraday Trans. 1990, 86, 1167. (41) JCPDS pattern no. 42-1164. (42) Orlando, R.; Pisani, C. Roetti, C. Phys. Rev. B. 1992-II, 45, 592. (43) Stefanovich, E. V.; Shluger, A. L.; Catlow, C. R. A. Phys. Rev. B. 1994, 49, 11560.

3580

Chem. Mater., Vol. 21, No. 15, 2009

Xia et al.

Table 3. Binding Energies and Solution Energies Derived from Mott-Littleton Calculationsa isolated defect defect configuration d1 (A˚) d2 (A˚) d3 (A˚) d3rel (A˚) EDEF (eV) Ebind per Y (eV)

Vo

14.78

Y

34.81

Y-Vo NN

Y-Vo-Y cluster NNN

2.195

4.204

49.67 0.07 0.18b -0.10c

49.26 -0.33 -0.26b -0.44c

Es per Y (eV)

NN

NNN

NN-NNN

2.195 2.195 3.585 3.608 83.62

4.204 4.204 3.585 3.608 82.98

4.204 4.204 8.016 7.980 83.59

2.195 4.204 3.585 3.543 83.45

2.195 4.204 5.07 5.118 84.06

2.195 4.204 6.209 6.276 84.03

-0.26

-0.58

-0.27

-0.34

-0.04

-0.05

a d1 and d2 are the distances between yttriums and the center of oxygen vacancy (Y-V); d3 and d3rel are the distance of two yttriums (Y-Y) before and after relaxation in the various configurations of the defect cluster. b From ref 32. c From ref 47

Figure 2. Lattice energy as a function of yttrium concentration.

is calculated to be -135.39 eV and is employed in the subsequent calculation of solution energies in the YSZ system. 3.2. Y-Doped ZrO2. For the Y-doped zirconia system, according to eq 1, two Y3+ ions substitute for two Zr4+ ions to create one oxygen vacancy. This pair of yttrium dopants together with the charge compensating oxygen vacancy is considered in these calculations as a bonded defect cluster (Y-V-Y). Different configurations of the Y-V-Y defect cluster have been studied by locating the cation defects either in nearest-neighbor (NN) sites or next-nearest-neighbor (NNN) sites (2.2 A˚ and 4.2 A˚, respectively) to the oxygen vacancy. These simple pair clusters showing NN and NNN sites configurations are illustrated in Figure 1. As presented in Figure 1, various cluster configurations (three series) have been considered in the Mott-Littleton calculations for the infinitely dilute system. The solution energy derived from Born-Haber cycle and the binding

energy obtained by comparison with isolated defects have been summarized in Table 3. In particular, the NNN site yttrium-vacancy (Y-V) cluster which has a negative binding energy of -0.33 eV is favored by 0.40 eV with respect to the NN site (Y-V) cluster. The preference for yttrium to occupy the next-nearest-neighbor (NNN) site to the oxygen vacancy generally agrees with early reported experimental observations44-46 as well as the previous interatomic study by Khan et al.32 and ab initio calculations by Stapper et al.47 However, for the charge neutral Y-V-Y cluster, the energically favorable configuration can be indentified here as two yttrium next to each other and both occuping the next-nearest neighbor (44) Li, P.; Chen, W.; Penner-Hahn Phys. Rev. B. 1993, 48, 10063. J. Am. Ceram. Soc. 1994, 77, 118. (45) Catlow, C. R. A.; Chadwick, A. V.; Greaves, G. N.; Monroney, L. M. J. Am. Ceram. Soc. 1986, 69(3), 272. (46) Komyoji, D.; Yoshiasa, A.; Moriga, T.; Emura, S.; Kanamaru, F.; Koto, K. Solid State Ionics 1992, 50, 291. (47) Stapper, G. Bernasconi, M. Phys. Rev. B. 1999-II, 59, 797.

Article

Chem. Mater., Vol. 21, No. 15, 2009

3581

Figure 3. Schematic diagrams for the calculated YSZ cubic-like structures. (Zr, blue-green; Y, yellow; O, red). Table 4. Comparisons of Undoped ZrO2 and Y-Doped ZrO2 Cubiclike Structures parameter a (A˚) b (A˚) c (A˚) R (deg) β (deg) γ (deg) lattice energy per Zr(1-X)YXO(2-X/2) (eV)

2  2  2 Zr32O64

structure I Zr26Y6O61

Δ%

structure II Zr24Y8O60

10.14 10.14 10.14 90 90 90 -109.76

10.25 10.27 10.26 90.00 90.04 90.00 -101.94

1.12 1.24 1.18

10.33 10.33 10.33 90.00 90.00 89.69 -99.21

0.04

Δ% 1.90 1.87 1.85 -0.15

Table 5. Calculated Surface Energies of Zirconia by Interatomic Potential (IP) Approach and Previous Hartee-Fork (HF) Study ZrO2 IP (J/m2 )

ZrO2 HF (J/m2) (previous work)48

Miller index

region I: II (A˚)

unrelaxed

relaxed

(111) (110) (100) (310)

6:12 6:12 6:9 6:20

1.448 3.638 10.062 13.142

1.210 2.147 2.704 2.955

sites to the oxygen vacancy. Such a configuration gives a solution energy of -0.58 eV per yttrium in contrast to the considering more negative value of -3.1 eV suggested by Khan et al. The contradiction with the reported solution energies may be due to different values being taken for the lattice energies of the solute; in our case we employ the energy for the cubic structure which has been fully relaxed employing the reported interatomic potentials. We consider our value to be more realistic because the figure from Khan’s work would represent a strongly exothermic process. In addition to the Mott-Littleton approach, the YSZ systems with varying Y concentrations were investigated by introducing one defect cluster to the cell and extending the unit cell size from a 111 (Y2Zr2O3∼50% Y) to a 4 2 2 (Y2Zr62O127∼3% Y) supercell. Studies of phonon dispersion curves show negative frequencies for the YSZ systems especially those with YO1.5 content less than 5% mol, which corresponds experimentally to the monoclinic region in the Y2O3-ZrO2 phase diagram.16 Figure 2 plots lattice energy as a function of YO1.5 concentration. Three series of data are presented by choosing the value for the lowest energy configuration in each series. According to Figure 2, there is a clear trend that the lattice energy linearly increases with yttrium doping level. Nevertheless,

unrelaxed 1.514 3.035

relaxed 1.485 2.406

the lattice energies for YSZ systems with the same yttrium concentration are slightly influenced by the configuration of the defect cluster. Which site, NN or NNN, is more favorable for yttrium varies in each case, the energy differences being within (0.1 eV per unit cell. To investigate the experimentally claimed cubic-like structure of PSZ (Y2O3 concentration g8 mol %), 222 supercells doped with 3 and 4 Y-V-Y defect clusters have been employed equivalent to doping level of 18.75% and 25% mol YO1.5. After relaxation, two cubiclike YSZ structures were formed with the unit cell expanded 1%2% compared to pure ZrO2 as shown in Figure 3. The stability of these cubic-like structures has been confirmed by examination of the phonon dispersion curves, no negative frequency observed. Table 4 compares the two yttrium doped structures with pure cubic zirconia. Structure I splits the 6 Y dopants into 4 and 2 and locates these two groups of dopants alternately with a single ZrO2 slab between them. For structure II, with a 25% YO1.5 (Figure 3-II) content, i.e., fully stabilized material, the oxygen vacancies are located in four corner sites of the supercell and yttrium occupies the NNN sites with respect to the oxygen vacancies. 3.3. Surface Structure. A. Pure Zirconia. To examine the surface structure of cubic ZrO2, we studied the

3582

Chem. Mater., Vol. 21, No. 15, 2009

Xia et al.

Figure 4. Plots of Y solution energy and surface energy versus Y concentration for the (111) and (110) surfaces of c-ZrO2. NN defect clusters are located in the top layer of the surface and orientated in the X direction.

Figure 5. Schematic representation of the defect cluster moving down from the top of the c-ZrO2 (111) surface in a sequence of 1, 2, 3, 4, etc., sites. (Zr, blue-green; Y, yellow; O, red; O vacancy, white).

three low index surfaces (111), (110), and (100) together with the (310) surface because according to Balducci et al.,24 the (310) surface is suggested to have a comparable surface energy to that of the (100) surface after relaxation. The simulations were carried out by the MARVIN code using a 11 surface repeat. Oxygen (O)-terminated faces have been chosen for the calculations, since they were identified as being more stable than Zr or O-O terminations from the first-principles studies of Eichler et al.27 The dipolar type III (100) surface was reconstructed by removing the dipole moment as in ref 36. The calculated surface energies are presented in Table 5, which clearly predicts that the (111) and (110) surfaces are the most stable for cubic ZrO2 with the (111) surface being dominant because the energy is approximately 1 J/m2 lower. The reconstructed (100) surface after relaxation is more energetically favored than before but still less stable than the (111) or (110) surfaces. For the (310) surface, following a large displacement of atoms in the top layers (about 4 A˚ depth), the surface energy of the optimized structure becomes comparable to the (100) surface. Our study

precisely reproduces the results reported by Islam et al. and the order of surface stability for c-ZrO2, (111) > (110) > (100) > (310), is in agreement with many previous theoretical studies.20,24,27,48 B. Y-Doped ZrO2. For the YSZ system, the surface energies and yttrium solution energies were calculated as a function of the defect cluster depth and orientation and yttrium concentration at the surface. Two Y atoms and one oxygen vacancy are introduced as a defect cluster (Y-V-Y) throughout the surface calculations as well as in the bulk. This study focused on the two most energetically favorable surfaces of c-ZrO2, i.e., the (111) and (110). The solution and surface energies for different yttrium concentrations (CY) in the outmost surface layer of YSZ system are plotted in Figure 4. The variation of yttrium concentration (CY) was achieved by extending the surface repeat unit from (11) to (23). Energetically favorable yttrium solution energies were obtained at a Y surface content level of less than 25% for the (111) and (110) surfaces, which is accords with the earlier atomistic study of the Ca/ZrO2 system carried out by Kenway et al., they suggested that Ca2þ cations also have favorable solution energies for the (111)c and (110)c surfaces with all compositions, and segregate at the surface.49 Compared to the (111) and (110) surfaces of undoped ZrO2 (Table 5), the surface energies were reduced slightly over the same yttrium content region, which indicates that both these two surfaces of c-ZrO2 (48) Gernnard, S.; Cora, F.; Catlow, C. R. A. J. Phys. Chem. 1999, 103, 10158. (49) Kenway, P. R.; Oliver, P. M.; Parker, S. C.; Sayle, D. C.; Sayle, T. X. T.; Titiloye, J. O. Mol. Simul. 1992, 9, 83.

Article

Chem. Mater., Vol. 21, No. 15, 2009

3583

Figure 6. Defect solution energy versus cluster depth for the 1  1 and 12 repeat units of the (111) surface. The trend line presents the solution energy of the most favorable configuration at each depth.

The results in Figure 4 show that yttrium solution at the surface becomes decreasingly favorable with increasing Y concentration; and indeed, both the defect solution and the surface energies have a strong dependence on the yttrium content in the surface layer of YSZ. The tendency for yttrium segregation at the (111) and (110) surfaces of c-ZrO2 was investigated directly by calculating the segregation energy (Eseg) and trap energy (Etrap). Eseg is defined as the energy difference between the defect cluster in the bulk crystal (Edef bulk), or at a consider50 able depth below the surface, and at the surface (Edef surf) Figure 7. Schematic representation of the defect cluster moving from the top of the c-ZrO2 (110) surface in sequence of 1, 2, etc., sites. (Zr, bluegreen; Y, yellow; O, red; O vacancy, white).

def def Eseg ¼ Esurf - Ebulk

can be stabilized by yttrium doping. However, the relative stabilities of the (111) and (110) surfaces do not change.

(50) Redfern, S. E.; Stanek, C. R.; Grimes, R. W.; Rawlings, R. D. Philos. Mag. Lett. 2005, 85, 445.

ð6Þ

3584

Chem. Mater., Vol. 21, No. 15, 2009

Xia et al.

Figure 8. Plots of defect solution energy as a function of defect cluster depth for the 1  1 and 1  2 repeat units of the (110) surfaces. Table 6. Summary of Segregation and Trap Energies for the Y-Doped (111) Surface of c-ZrO2a surface unit

defect orientation (axis)

defect configuration

ESEG per Y (eV)

ETAPper Y (eV)

segregation depth (A˚)

11 11 1 2 12 11 11 12 1 2

X X X X Z Z Z Z

NN NN-NNN NN NNN NN NNN NN NNN

0.7 -0.5 -0.3 -0.7 -0.6 -0.6 -0.5 -0.7

-1.1

∼5 3-4 ∼5 ∼5 3-4 ∼5 3-4 3-4

a

-0.2 -0.1

Segregation depth indicated the depth of the energetically favorable sites for yttrium to the top of the surface.

Similarly, ET is defined as the difference between the energy of the defect cluster at the surface (Edef surf) and in the lowest energy site (Edef trap) def def Etrap ¼ Etrap - Esurf

ð7Þ

A negative value for Eseg indicates the surface sites are more favorable than the bulk and segregation can occur. The trap energy Etrap can be zero, which means no trap sites are found, or negative, i.e., there are more stable sites near the surface but not at the surface.

The segregation energy was evaluated by moving the Y-V-Y clusters from the top layer of the surfaces toward the bottom of region I. Both 11 and 12 surface repeat units were examined. We studied two series of cluster sets with different defect orientations, i.e., Y-VY along the X and Z directions, were studied with nearest neighbor (NN) or next-nearest neighbor (NNN) configurations as displayed in Figures 5 and 7. The variation of yttrium solution energy versus the depth of these defect clusters from the surface is presented in Figures 6 and 8 for the (111) and (110) surfaces respectively.

Article

The most significant results for the (111) surface are listed in Tables 6. In general, for both X and Z orientations of the Y-V-Y cluster, the solution energy no longer depends strongly on cluster orientation at distances of 5 A˚ or more from the surface (Figure 6.). However, as the cluster approaches the surface, a strong configurational dependence is apparent. The surfaces with the defect cluster in the Z orientation (Figure 5-II) possess negative segregation energies up to -0.7 eV, and hence the defect cluster tends to locate along the Æ100æ direction on the (111)c surface. An enrichment of yttrium toward the (111) surface would occur with formation of a higher concentration layer extending over 5 A˚ deep (depth equivalent to four oxygen layers). The values of Eseg in Table 6 are mostly negative, which is consistent with the results of the surface energy calculations. One notable result in Table 6 concerns the high yttrium concentration level of 50%, i.e., for the 11 surface repeat unit with the defect cluster in the X orientation where a positive segregation energy in the outmost layer accompanied by a large energy trap in the subsurface layer are indicated by the calculations for the NN sites configuration. In such a case, the yttrium tend to disperse to the NNN sites, so the defect clusters are orientated with the vacancy at the top surface or subsurface (i.e., Z direction). These results suggest that saturation of yttrium has been achieved at the outmost layer with an Y/Zr ratio 1:1 (12.5% mol Y2O3) and thus the defect clusters prefer to locate in the second and third layers of the subsurface. Therefore, from these calculations we can conclude that yttrium segregates at the (111) surface up to a depth of the top four layers. A series of calculations for the less stable (110) surface was performed in a similar manner (Figure 7.) to those for the (111) surface. Figure 8 presents the solution energy as a function of the defect depth using 11 and 12 surface repeat units. In these models, the NNN defect sites adopt an identical configuration to the NN sites after relaxation. For the X orientation, only the 1  2 repeat can be considered. In this case, with substitution of yttrium at half of the Zr sites at the (110) surface layer, the segregation energy oscillates, and no significant trap energy is observed. Such a variation in ES with defect depth is ascribed to the cluster dipole interacting with the local dipole along the surface. A slightly lower yttrium solution energy is estimated for the outmost layer of the (110) surface for Y-V-Y cluster in the X orientation. However, there is no strong indication (0 > ESEG > -0.3 eV) that the yttrium dopants will segregate to the (110) surface. Equally, Y solution energies as a function of depth for the Z orientation do not suggest surface segregation. These results from atomistic calculations are consistent with the previous first-principles study, which also indicates little driving force for yttrium segregation at the (001) surface in t-ZrO2,27 which is equivalent to the (110) surfaces of c-ZrO2. Overall, from the above results, it is expected that an enrichment of yttrium toward the (111) surface would

Chem. Mater., Vol. 21, No. 15, 2009

3585

occur with the formation of a higher concentration layer extending over 4-5 A˚ (i.e., four layers from the top) and no clear driving force for Y segregation at the (110) surface. This finding is consistent both with previous calculations of t-ZrO227,50 and c-ZrO228 as well as several experimental observations3,19-23. In particular, the production of an upper limit to the level of Y at the surface is in accordance with the experimental study by Zhu et al., which suggested the composition of the outmost surface YSZ is independent of Y concentration in the bulk and reaches a limit of 12 ( 2 mol % Y2O3 because of segregation of Y2O3.3 The accompanying surface segregation of the charge compensating of oxygen anion vacancies will have a substantial effect on surface chemistry. 4. Summary and Conclusions The simulation studies reported in this paper have successfully modeled the bulk cubic phase ZrO2 and the (111), (110), (100), and (310) surfaces of c-ZrO2. For bulk supercell calculation in the range 9-12 mol % yttria, two stable cubiclike phases were indicated differing in their detailed oxygen arrangement in line with experimental observations. The order of surface stability for c-ZrO2: (111)c > (110)c > (100)c > (310)c was confirmed by the calculations, i.e. the (111)c surface is the dominant surface of undoped c-ZrO2 in good agreement with previous experimental and theoretical studies. Doping with yttrium stabilizes the cubic phase ZrO2 and the energy of the YSZ lattice is found to linearly increase with yttrium concentration in the bulk. According to the Mott-Littleton calculations, yttrium preferentially occupies the next nearest neighbor (NNN) sites to the vacancy rather than nearest neighbor (NN) sites; in addition, the two yttriums tend to be located together. However, the supercell calculations suggest that the dispersion of yttrium with respect to the oxygen vacancy depends on yttria content level and the NN and NNN sites are almost equally favorable. The energetically favorable (111) and (110) surfaces of pure cubic zirconia could be further stabilized by the yttrium doping. However, doping does not change the relative stabilities of the (111)c and (110)c surfaces. Yttrium tends to segregate at the top layers (up to 5 A˚) of the dominant (111)c surface. In contrast, no clear tendency for yttrium segregation of yttrium has been observed on the (110)c surface of YSZ. The present work has advanced our understanding of the effect of yttrium doping and the creation of oxygen vacancies in ZrO2 surfaces. An investigation of the electronic structure of ZrO2 and YSZ surfaces using quantum mechanical method will be reported subsequently. Acknowledgment. We are grateful for the financial support from University College London and the computer facilities funded by UCL or EPSRC.