Protonated Forms of Monoclinic Zirconia: A Theoretical Study - The

8014

J. Phys. Chem. C 2010, 114, 8014–8025

Protonated Forms of Monoclinic Zirconia: A Theoretical Study Yves A. Mantz* and Randall S. Gemmen U.S. Department of Energy, National Energy Technology Laboratory, 3610 Collins Ferry Road and PO Box 880, Morgantown, West Virginia 26507 ReceiVed: December 2, 2008; ReVised Manuscript ReceiVed: February 22, 2010

In various materials applications of zirconia, protonated forms of monoclinic zirconia may be formed, motivating their study within the framework of density-functional theory. Using the HCTH/120 exchange-correlation functional, the equations of state of yttria and of the three low-pressure zirconia polymorphs are computed, to verify our approach. Next, the favored charge state of a hydrogen atom in monoclinic zirconia is shown to be positive for all Fermi-level energies in the band gap, by the computation of defect formation energies. This result is consistent with a single previous theoretical prediction at midgap as well as muonium spectroscopy experiments. For the formally positively (+1e) charged system of a proton in monoclinic zirconia (with a homogeneous neutralizing background charge density implicitly included), modeled using up to a 3 × 3 × 3 arrangement of unit cells, different stable and metastable structures are identified. They are similar to those structures previously proposed for the neutral system of hydrogen-doped monoclinic zirconia, at a similar level of theory. As predicted using the HCTH/120 functional, the lowest energy structure of the proton bonded to one of the two available oxygen atom types, O1, is favored by 0.39 eV compared to that of the proton bonded to O2. The rate of proton transfer between O1 ions is slower than that for hydrogen-doped monoclinic zirconia, whose transition-state structures may be lowered in energy by the extra electron. 1. Introduction Zirconia (zirconium dioxide, ZrO2) is a material that, due to its strength and thermal stability, is useful for a variety of applications. In some of these applications, protonated forms of monoclinic zirconia can form, motivating their study both experimentally and theoretically. One example is in solid oxide fuel cells during the process of fuel oxidation. These devices, typically operating between 770-1270 K, are composed of an anode [generally a Ni-YSZ cermet, where YSZ is yttria (Y2O3)-stabilized zirconia], an electrolyte (generally YSZ), and a cathode (such as Sr-doped lanthanum manganate, LaMnO3), along with other components. For the case of molecular hydrogen as fuel, the overall process of fuel oxidation is 2H2,gas + OYSZ f H2Ogas + 2eNi

(1)

This electrochemical conversion occurs by different competing pathways or reaction schemes,1-8 such as hydrogen spillover or interstitial hydrogen transfer, where each pathway includes a charge-transfer step. Electrochemical impedance spectroscopy is used to characterize these steps and help formulate hypotheses as to the most important pathway(s).9,10 To test these hypotheses, kinetic modeling can be applied in order to predict impedance spectra for a particular (set of) pathway(s).1,11-13 A goal is to construct kinetic models that are able to represent all pathways as accurately as possible. Thus, it is essential to understand the thermodynamics/kinetics of the reactions comprising these pathways, including reactions or processes that might affect the concentrations of species involved in the charge-transfer steps.1 One such pathway, previously applied to interpret impedance spectra,9,10 is interstitial hydrogen transfer. It consists of the * E-mail: [email protected].

10.1021/jp810601j

dissolution of adsorbed hydrogen atoms into the Ni anode, Hi,Ni, followed by their oxidation and transference across the planar + , anode/electrolyte boundary, forming protons in YSZ, Hi,YSZ according to + Hi,Ni f Hi,YSZ + e-

(2)

These protons may then react with oxygen ions in YSZ to form + water.1,9,10 The concentration of Hi,YSZ , and therefore the rate of the charge-transfer step, reaction 2,1 will be affected if the diffusion of the proton in the electrolyte is slow. Because the surface region of the YSZ electrolyte appears dominated by a monoclinic zirconia phase ≈6-nm thick,14 an understanding of hydrogen dissolution and kinetics of the proton in monoclinic zirconia is needed for a better understanding of this mechanistic pathway. Other examples of applications in which protonated monoclinic zirconia can form are quite diverse. The cladding material of light water nuclear reactors can form a layer of monoclinic zirconia. To minimize the hydrogen embrittlement of the cladding material, an understanding of hydrogen solubility and diffusion in monoclinic zirconia is required.15 Lastly, monoclinic zirconia is a candidate as a high-k dielectric replacing silica (silicon dioxide, SiO2) as the gate dielectric in field effect transistors. Because hydrogen is a ubiquitous impurity affecting the performance of this and other prospective materials for this application, the behavior of hydrogen in oxides must be understood.16 The development of a theory of hydrogen in oxides may play a key role for the improvement of microelectronics devices.17-23 Our knowledge of hydrogen in monoclinic zirconia is appreciable. In contrast to a proton in YSZ,10,24,25 a proton in monoclinic zirconia is difficult to probe experimentally due to a lack of oxygen vacancies permitting the incorporation of

This article not subject to U.S. Copyright. Published 2010 by the American Chemical Society Published on Web 04/07/2010

Protonated Forms of Monoclinic Zirconia

J. Phys. Chem. C, Vol. 114, No. 17, 2010 8015

protons (in the form of hydroxide ions filling the vacant sites) via the dissociative adsorption of water from the gas phase.26 Despite this fact, this system can be studied by various means. If monoclinic zirconia is placed in an oxygen- and water-rich atmosphere, interstitial oxygen ions, O2i,m-ZrO2, interstitial protons, + , and electron holes, h+, will form:15 Hi,m-ZrO 2

1 2h Oi,m-ZrO + 2h+ O 2 2 2,gas

(3)

1 + H2Ogas + 2h+ h 2Hi,m-ZrO + O2,gas 2 2

(4)

By taking advantage of the equilibrium in reactions 3 and 4, hydrogen can be shown to possess a limited solubility in monoclinic zirconia.15 In a different experimental study, using the technique of muonium spectroscopy, essentially a form of ion implantation, with positive muons (Mu+’s) in place of protons, combined with a form of magnetic resonance known as muon spin rotation, a shallow Mu+ state is observed in monoclinic zirconia.22,23 Thus, it can be concluded that a hydrogen atom will relinquish its electron in pure (undoped) monoclinic zirconia. This conclusion is substantiated by theoretical studies. By performing static computations as well as dynamical simulations at the level of density-functional theory (DFT), stable structures of a hydrogen atom in monoclinic zirconia can be identified, resulting in hydrogen relinquishing its electron and reacting with an oxygen ion to form a hydroxide ion, and a diffusion mechanism proposed involving the hopping of a proton from an oxygen ion of a certain type to another of the same type.16 Furthermore, a theory of hydrogen in oxides has been developed to explain the creation of hydrogen in different charge states in a variety of oxide materials.19-21 A central point is that the hydrogen atom will donate its electron if the conduction band edge lies below some critical energy above the valence band, i.e., if the band gap is less than a critical value. This band gap criterion is also a key aspect of a second theory more recently proposed.22,23 To enhance our understanding of proton transfer in monoclinic zirconia and further test these theories of hydrogen in oxides, in this work the thermodynamics of protonated monoclinic zirconia is characterized in detail, and the dynamics of the proton is studied at different temperatures at the level of DFT. The organization is as follows: In section 2, a description of the methodology is given. The choice of methodology is extensively validated by performing computations of yttria and different zirconia polymorphs that are summarized in section 3. The remainder of section 3 is devoted to static computations of protonated monoclinic zirconia. In section 4, dynamical simulations of monoclinic zirconia-based systems are described. Conclusions are drawn in section 5. 2. Computational Approach Three types of computationsswave function optimizations, geometry optimizations, and Born-Oppenheimer molecular dynamics (BOMD) simulationssare performed in this work. Technical details of the computations are described below. The CPMD v3.11 software package is used throughout.27,28 The CPMD computations are performed at the level of DFT, within the generalized gradient approximation (GGA), using the HCTH/120 exchange-correlation functional,29 hereafter, HCTH. The wave function is expanded in a plane-wave basis set with a kinetic energy cutoff of 80 Ry. The core electrons and nuclei

are described using Troullier-Martins norm-conserving pseudopotentials.30 The electronic configuration previously chosen to generate the hydrogen pseudopotential is 1s0.85 (mathematically well-defined), the oxygen pseudopotential is [He] 2s22p4, and the zirconium pseudopotential is [Kr] 4d2 (+2 ion).28 Note that the semicore 4s and 4p as well as the 4d electrons are treated as valence states in the generation of the zirconium pseudopotential, which is important to ensure the transferability of this pseudopotential. The wave functions of a variety of systems are optimized using direct methods, which can be more efficient than iterative diagonalization methods when using a large (plane-wave) basis set.27 Computations are performed at the level of spin-restricted DFT, using the default algorithm, direct inversion in the iterative subspace (DIIS), with 10 vectors.27 The convergence criterion for the maximum component of the wave function gradient, defined as the partial derivative of the Kohn-Sham (potential) energy with respect to a wave function coefficient,31 is set equal to either 5 × 10-7 or 10-6 Ha, yielding well converged Kohn-Sham (potential) energies. In addition, geometry optimizations are performed. For the perfect (periodic) crystalline models studied, the limited-memory Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm32 is employed along with adaptive convergence criteria. The maximum component of the wave function gradient for the initial geometry is set equal to 10-4 Ha, the ratio of the maximum component of the wave function gradient to the maximum component of the force on any ion (reached so far) is set equal to 0.02 Bohr, and this criterion is switched off once the maximum component of the wave function gradient is 10-5 Ha. Thus, fewer DIIS steps are performed during earlier steps of the geometry optimization. For the defective (periodic) model systems of protonated monoclinic zirconia, with several shallow minima of potential energy, adaptive convergence criteria are not used. Instead, standard values for the convergence criteria are selected, namely 10-6 Ha for the maximum component of the wave function gradient and 5 × 10-4 Ha/Bohr for the maximum component of the force on any ion. In addition, DFT-based BOMD simulations are performed within the microcanonical (NVE) ensemble. During a BOMD simulation, the wave function is optimized at every time step.27 To optimize the wave function, an initial guess is extrapolated from the three previous wave functions, resulting in a factor of 2-3 enhanced efficiency. During selected optimizations, the irreducible Brillouin zone is sampled using a Monkhorst-Pack33 k-point sampling scheme. However, the BOMD simulations are performed strictly at the Γ point, which is a reasonable approximation (vide infra). In Appendix A, a model system, molecular hydrogen, is examined. The minimum-energy geometry of this species is predicted using the default method of DIIS, with five vectors, combined with a quasi-Newton method (using BFGS). Also, BOMD simulations are performed, during which different values of the time step are sampled, the wave function convergence criterion is fixed at 10-6 Ha, and wave function extrapolation is not performed. Additional computational details for other systems that are studied in the Supporting Information are described there. The dependence of certain results on the choice of functional and k-point grid density is assessed using the Vienna ab initio simulation package (VASP).34-37 This program evaluates the total energy of periodically repeated geometries based on DFT and the pseudopotential approximation. Computations are performed at the GGA-PW91 level,38,39 with a plane-wave cutoff

8016

J. Phys. Chem. C, Vol. 114, No. 17, 2010

Mantz and Gemmen

Figure 1. Minimum volume (in black) and conventional (in gray) unit cells of the three ambient or low-pressure zirconia polymorphs. For monoclinic zirconia, zirconium atoms are in gray, and the two types of oxygen atoms O1 and O2 are in red and green, respectively.

of 400 eV, while the electron-ion interaction is described using the projector augmented wave (PAW) pseudopotential method.40,41 3. Static Properties of Zirconia-Based Systems Zirconia can exist in three ambient or low-pressure phases (Figure 1). The monoclinic (“baddeleyite”) phase is of space group C2h5 or P21/c (no. 14).42 The unit cell contains 12 atoms and three atom types, Zr, O1, and O2, following the established labeling convention. Thus, it is completely specified by four lattice parameters as well as nine Wyckoff atomic coordinates. This phase is stable from 0 to ≈1400 K.43 The tetragonal phase of zirconia is of space group D4h15 or P42/nmc (no. 137)42 and is stable from ≈1400 to ≈2570 K.43 The cubic (“fluorite”) phase is of space group Oh5 or Fm3m (no. 225)42 and is stable from ≈2570 K to the melting temperature of 2983 K.43 In this section, a validation of the methodology is given. Next, the favored charge state of a hydrogen atom in monoclinic zirconia is predicted. Selected stable and metastable structures are described and further characterized by performing a Bader atomic charge and a harmonic vibrational analysis, the latter providing an estimate of entropic effects at finite temperature. Several tests are performed to establish that the DFT-GGAHCTH methodology, previously applied with success to study liquid water, including hydrogen bond breaking and formation events as well as self-diffusion,29,44 is appropriate for characterizing the behavior of a proton in monoclinic zirconia. The lattice constants of monoclinic zirconia predicted by scanning the potential energy surface, i.e., by incrementing a, b, c, and β independently with fractional atomic coordinates held fixed and performing wave function optimizations, are a ) 5.190 Å, b ) 5.240 Å, c ) 5.370 Å, and β ) 99.0°, overestimating by only 1% experimental results extrapolated to 0 K, a ) 5.137 Å, b ) 5.199 Å, c ) 5.330 Å, and β ) 99.0°.43,45 Furthermore, the Murnaghan equations of state46 of cubic, tetragonal, and monoclinic zirconia are predicted (Figure 2). Results for monoclinic zirconia are in reasonable agreement with those previously obtained using a GGA functional.47,48 The energy differences between monoclinic-tetragonal and monoclinic-cubic phases predicted here, 140 and 320 meV/ZrO2, respectively (Figure 2), as well as elsewhere47,48 are larger by a factor of 2-3 than those obtained within the local density approximation,49-52 which can reproduce quite accurately the measured enthalpy difference at the monoclinic-tetragonal phase transition temperature, 63 meV/ZrO2, as well as the sum of enthalpy differences at the monoclinic-tetragonal and tetragonal-cubic phase boundaries, 120 meV/ZrO2.53 This appears to be a deficiency of GGA functionals, as previously suggested.47 Nonetheless, monoclinic zirconia is, correctly, predicted to be

Figure 2. Computed cohesive energy versus volume data fitted by Murnaghan equations of state for cubic (c-ZrO2), tetragonal (t-ZrO2), and monoclinic (m-ZrO2) zirconia.

most stable at 0 K, followed by tetragonal and then the cubic phase. Furthermore, a GGA functional is expected to be more accurate for describing protonated monoclinic zirconia and any bonds involving the proton,29,44 the main focus of this work. Thus, the selected methodology can be applied with confidence to characterize monoclinic zirconia-based systems. A detailed description of the validation procedure, results obtained, and comparison to previous theoretical and experimental results is provided in the Supporting Information. To study a hydrogen atom in monoclinic zirconia in various charge states, a realistic model system of monoclinic zirconia is constructed. First, the lattice constants of a single unit cell are set equal to values of a ) 5.1828 Å, b ) 5.2117 Å, c ) 5.3731 Å, and β ) 98.835° that are within 1% of those extrapolated to 0 K and are the same as the experimentally determined values at 1100 K.43,54 A geometry optimization is then performed, preserving the crystal symmetry, and utilizing a sufficiently dense 3 × 3 × 3 Monkhorst-Pack k-point mesh producing 10 nonsymmetric k points in the irreducible Brillouin zone. Next, a 2 × 2 × 2 unit cell arrangement is constructed from the optimized cell, and a second geometry optimization is performed. Each atomic coordinate from this geometry optimization, sampling only at the Γ point, is converged to within 5 × 10-4 Å of results obtained for the same model system but sampling over either two or 10 nonsymmetric k points in the irreducible zone. Furthermore, the predicted electronic density of states (Figure 3, dashed black curve) is in excellent agreement with previous theoretical results.47,55 The broad features at -16.5 and -3 eV are primarily due to mixing of O 2s and O 2p atomic orbitals, respectively; that at >3 eV is due to mixing of Zr 4d orbitals. A careful procedure is used to determine the initial placement of hydrogen in monoclinic zirconia, represented by the 96-atom


Figure 3. Electronic density of states of monoclinic zirconia, the configuration 1A (Figure 4), and a “transition state” structure taken from simulation #2 (section 4). For 1A, the two O-H bonding states are indicated by arrows. All curves are shifted by the same energy amount. For the first two structures, the top of the valence band is the zero of energy.

Figure 4. Minimum energy structure of a proton bonded to an oxygen atom in a 2 × 2 × 2 arrangement of monoclinic unit cells, denoted configuration 1A in the text. The proton (in blue) is bonded to an oxygen ion (in red) that would be of type O1 in pure monoclinic zirconia. This oxygen ion is located near the center of a plane formed by three Zr ions (and is connected to them by thick black lines). The interactions between the proton and its first two oxygen atom neighbors H · · · O1st and H · · · O2nd are shown (thin dashed black lines).

supercell described above. First, computations of a hydrogen atom in monoclinic zirconia are performed in which the net charge of the periodic system is taken to be q ) -1e, 0e, or +1e. In the first and third model systems, a homogeneous jellium background of neutralizing charge is added to avoid divergence of the (Coulombic) energy.56,57 Throughout this work, periodic models of a hydrogen atom in monoclinic zirconia placed in formal charge states of -1e, 0e, or +1e will be denoted H + m-ZrO2 (q ) -1, 0, or +1, respectively). Similar energetic minima are formed only by the hydrogen atom interacting with the different kinds of oxygen atoms, while the interaction between hydrogen and zirconium atoms is repulsive, for all three charge states q. This finding is consistent with a previous theoretical study of the system H + m-ZrO2 (q ) 0), in which a hydrogen atom interstitial not bonded to an oxygen ion was metastable by 0.95 eV with respect to a hydrogen atom donating its electron into the conduction band and making a dative bond to O1, forming (OH- + e-).16 Thus, the hydrogen atom is placed manually along the vectors formed by either O1 or O2 and other oxygen atoms within 3.83 Å of these ions, sampling thoroughly the volume around both of them. The initial O1-H or O2-H bond length is chosen to be 1.00 Å. In all, 13 initial configurations of a hydrogen atom bonded to O1, and 13 of it bonded to O2, are constructed, and each is optimized for q ) -1, 0, and +1. The most stable configuration of H + m-ZrO2 (q ) +1) identified is designated 1A (Figure 4). It is quite similar geometrically to the most stable configurations of H + m-ZrO2 (q ) 0) and H + m-ZrO2 (q ) -1). In particular, the optimized

J. Phys. Chem. C, Vol. 114, No. 17, 2010 8017 O-H distances are 0.99 Å (q ) +1) versus 0.98 Å (q ) 0 or -1), while the shortest hydrogen bond is 1.91 Å (q ) +1) versus 1.97-1.98 Å (q ) 0 or -1). The fact that all three geometries, obtained at the Γ point, are essentially the same can be explained by orbital arguments. As is true for a variety of protonated wide-gap oxides, the O-H bond is the strongest bond in the system, such that the O-H antibonding σ* state lies above the conduction band minimum.19-21 Therefore, the added electron or two in the systems H + m-ZrO2 (q ) 0 or -1, respectively) will relax to the conduction band minimum and leave unperturbed the local electronic structure around the O-H moiety. Indeed, two new peaks are observed in the electronic density of states at -19.4 and -7.0 eV (Figure 3, solid black curve, indicated by arrows) corresponding to H 1s-O 2s and H 1s-O 2p bonding states, respectively, as confirmed by a projection of these states onto pseudo atomic orbitals. The situation is different in cubic zirconia, with presumably a similar O-H bond strength19,20 but a slightly larger band gap than that of monoclinic zirconia (although this latter point is debatable),49,58-60 resulting in a population of the O-H σ* state and an O-H bond elongation from 1.02 Å (q ) +1) to ≈1.07 Å (q ) 0 or -1).21 Note that the density-of-states plots obtained for H + m-ZrO2 (q ) +1, 0, or -1) are essentially identical, implying that the effect of the neutralizing jellium background is small. To compare the relative stabilities of the lowest energy H + m-ZrO2 (q ) -1, 0, and +1) configurations, defect formation energies ∆EHf q’s are determined. This quantity can be defined as the total energy of the system containing a defect minus the total energy of the perfect system:17,61-64 f

∆EHq ) Etot(Hq) - Etot(0) - µH + qµe

(5)

where Etot(Hq) is the total energy of the defective system in charge state q, Etot(0) is the total energy of the neutral perfect system, µH is the hydrogen atom’s chemical potential, which is dependent on environmental conditions, and µe is the chemical potential of electrons in the reservoir, i.e., the neutral perfect system, with which electrons are exchanged. An example of a defective system is the configuration 1A with q ) +1. The perfect system is neutral monoclinic zirconia. Under hydrogenrich conditions, µH is taken as half the total energy of a H2 molecule:17,63,64

1 µH ) Etot(H2) 2

(H-rich)

(6)

Under water- and oxygen-rich conditions, it is defined with respect to the total energy of an isolated H2O and O2 molecule:63,64

µH )

1 1 tot E (H2O) - Etot(O2) 2 2

[

]

(O-rich)

(7)

Lastly, in a semiconductor or insulator, µe can assume values ranging from the energy of the valence band maximum Evbm to the energy of the conduction band minimum Ecbm.65 These limiting values of µe can be obtained as Evbm ) [Etot(0) Etot(+1)] and Ecbm ) [Etot(-1) - Etot(0)],63,64 where Etot(q) denotes the total energy of a perfect lattice supercell after removing q electrons.66

8018


Mantz and Gemmen

Figure 5. Defect formation energies, ∆EHf q, as a function of the Fermi energy µe and the hydrogen chemical potential µH, i.e., under either hydrogen- or oxygen-rich conditions, computed at the Γ point and at the HCTH level of theory.

Figure 6. Defect formation energies, ∆EHf q, as a function of the Fermi energy µe and the hydrogen chemical potential µH, i.e., under either hydrogen- or oxygen-rich conditions, computed over 16-32 k points at the PW91 level of theory.

The formation of a proton in monoclinic zirconia is computed to be favorable for all values of the chemical potential of electrons, i.e., the Fermi energy, in the theoretical band gap (Figure 5). An upper limit to the band gap energy Eg can be obtained as the energy difference between the lowest occupied and highest occupied one-electron states at the Γ point and is 3.88 eV. This value is less than the difference in total energies at the Γ point of the same 96-atom supercell model, Eg ) Ecbm - Evbm ) [Etot(-1) - Etot(0)] -[Etot(0) - Etot(+1)] ) 4.29 eV, which is in error due to orbital relaxations. Still, our band gap is smaller than very accurate theoretical values of 5.46 eV49 and 6.04 eV60 as well as the gap from experiment, 4.2-7.1 eV.58,59 This underestimation is typical of DFT and is due to electronic correlation effects.67 Note that no defect levels are induced by hydrogen within the band gap (Figure 3, solid black curve). Thus, if a band gap correction is applied, the defect energy levels that possess conduction band character will follow this upward shift, while those that possess valence band character will be left unchanged.61,63 That is, both ∆EHf 0 and ∆EHf -1 will be shifted to higher energies,61,63 and it can be concluded that the proton is the lowest energy state for all µe’s in the experimental band gap, as well. Lastly, note that the difference between eqs 6 and 7 is predicted to be 1.71 eV, in reasonable agreement with the experimental formation energy of an H2O molecule extrapolated to 0 K, 1.24 eV.64 The ∆EHf q’s obtained are based on energies computed at the Γ point, resulting in equal ∆EfHq’s at µe ) Eg, and would change minimally if additional k points were included. To illustrate this fact and obtain results using a more popular GGA functional, computations are performed using the VASP software package and the PW91 functional. Previously, this computational protocol was used to model successfully four different polytypes of zirconia.47 Similar results for cubic, tetragonal, and monoclinic zirconia are obtained using our computational setup with different pseudopotentials as in that study,47 serving as a validation. Next, the most stable structures of the proton in a 2 × 2 × 2 arrangement of monoclinic unit cells for different charge states of the system (q ) -1, 0, and +1), obtained at the HCTH level of theory, as well as those structures that are less stable by 1.5-2.0 eV, are reoptimized at the Γ point. They are minimally changed. These structures are further refined by sampling over a 2 × 2 × 2 Monkhorst-Pack k-point mesh/4 special k points in the irreducible Brillouin zone, and then over a 4 × 4 × 4 grid/32 special k points. The identity of the most stable configurations is preserved. The OH and nearest-neighbor hydrogen bond lengths of 0.99 Å and 1.90-1.93 Å, respectively, are similar to the values predicted at the Γ point (under either functional description). Energies for pure monoclinic zirconia in different charge states are converged using a 4 × 4 × 4 grid/

16 special k points. A band gap of 3.59 eV is determined as the energy difference between one-electron states, versus 3.70 eV computed as the difference in total energies. The ∆EHf q’s obtained at the Γ point (not shown) or sampling over 16-32 k points (Figure 6) are rather similar to those in Figure 5, because the exchange-correlation functional descriptions are similar, and because the symmetry properties of the electronic density are relatively unperturbed by the hydrogen dopant, causing the added k points to have little effect. The difference between eqs 6 and 7 is 1.56 eV at the PW91 level of theory. Other tests are performed using the HCTH functional. Finitesize effects are confirmed to be small by performing a test computation using a larger model system, consistent with the results of previous work.64,68 This finding is expected, given that the strength of the O-H bond is the main driving force for stabilization of the configurations.61 Including spin polarization, Etot(H0), Etot(+1), and Etot(-1) are lowered by only a few hundredths of an eV. Thus, a proton is predicted to be the favored charge state of a hydrogen atom in monoclinic zirconia for all values of µe. This is the reason that new peaks are not introduced in the electronic density of states by the addition of electrons to the system H + m-ZrO2 (q ) +1). That is, the (O-H)- moiety is quite stable, and the added electrons are placed into the bottom of the conduction band and not the O-H σ* state, thereby preventing the formation of a free hydrogen atom in the crystal and a 1s orbital state. This type of behavior was observed in other wide-gap oxides by Robertson and co-workers.19,21 The fact that a proton is predicted to be the favored charge state of a hydrogen atom in monoclinic zirconia for all values of µe is consistent with previous theoretical and experimental results.16,19,20,22,23 Interstitial hydrogen exhibits two basic types of behavior in semiconductors:17,21 It can form a deep-gap state and exist in all three charge states depending on the chemical potential of electrons, or it can act exclusively as a donor, an example of which is reported here (Figures 5 and 6). In the system H + m-ZrO2 (q ) 0), hydrogen is predicted to relinquish its electron,16 indicating that the positive charge state is favored at a single value of µe corresponding to the midgap, Eg/2. Using the experimental technique muonium spectroscopy, a shallow Mu+ state is observed in monoclinic zirconia, along with metastable Mu0 states.22,23 Given our computed defect formation energies of hydrogen in monoclinic zirconia as well as those of hydrogen in cubic zirconia,21 the chemical potential at which positive and negative charge states are equal in energy would be ≈5 eV above the valence band maximum after correcting for the band gap error, as anticipated by Peacock and Robertson.19,20 Thus, the present computations are providing a proof of hydrogen acting exclusively as a donor in monoclinic



TABLE 1: Stable and Metastable Configurations of the Proton Bonded to the Oxygen Atom O1 (or O2) in a 2 × 2 × 2 Arrangement of Monoclinic Unit Cells, Denoted 1A and 1B (or 2A, 2B, and 2C)a 1A 1B 2A 2B 2C

∆E

OX-H

H · · · O1st

H · · · O2nd

∠OX-H · · · O1st

∠OX-H · · · O2nd

νs

νb1

νb2

0.00 0.02 0.39 0.63 1.22

0.979 0.969 0.990 1.011 0.994

1.91 2.19 1.86 1.67 1.57

2.33 2.27 2.05 2.34 2.40

127.6 121.7 137.5 154.6 157.9

113.5 106.8 118.9 106.5 102.3

3505 3684 3234

968 894 1055

795 787 899

The energy ∆E of a structure relative to that of the configuration 1A is given in eV. The OX-H covalent bond, where X ) 1 or 2, and the distance of the proton from its two nearest oxygen atom neighbors H · · · O1st and H · · · O2nd are given in Å. The angles ∠OX-H · · · O1st and ∠OX-H · · · O2nd are given in degrees. The stretching and bending mode frequencies of the proton νs, νb1, and νb2 are given in cm-1. a

Figure 7. Clusters of atoms extracted from the periodic structures 1A, 1B, 2A, 2B, and 2C. In 1A and 1B, the proton (in blue) is bonded to an oxygen ion (in red) that would be of type O1 in pure monoclinic zirconia. This oxygen ion is located near the center of a plane formed by three Zr ions (and is connected to them by thick black lines). In 2A-2C, the proton is bonded to an oxygen ion (in green) that would be of type O2 in monoclinic zirconia. This oxygen ion is initially located near the center of a tetrahedron formed by four Zr ions (and is connected to them by thick black lines). The interactions between the proton and its first two oxygen atom neighbors H · · · O1st and H · · · O2nd are shown (thin dashed black lines).

zirconia, and the agreement with previous results can be taken as further validation of the models/methodology employed in this work. Having established that a proton is always favored, different stable and metastable configurations of H + m-ZrO2 (q ) +1) are identified (Table 1). The lowest energy configurations, labeled 1A and 1B, are formed by the proton bonded to the oxygen atom O1. In monoclinic zirconia, the 3-coordinated oxygen atom O1 is located approximately in a plane formed by three Zr ions (Figure 1).16 To minimize the repulsion of the proton from these cations, the O1-H bond formed is roughly perpendicular to this plane. Thus, there are two possible ways of proton attachment (Figure 7), as previously described for the system H + m-ZrO2 (q ) 0).16 The fact that similar configurations are obtained for these different systems is physically reasonable, as the extra electron present in the other system is occupying a delocalized state, the conduction band minimum. This comparison is an explicit verification for a system in the solid state that the homogeneous neutralizing background charge density is not influencing the predicted structures, as similarly demonstrated for aqueous-phase systems.69,70 Because the local environment above and below the plane of Zr ions is slightly different, configurations 1A and 1B are not isoenergetic. In configuration 1A, a fairly typical three-center bond is formed, in which the proton is bonded to three atoms, forming one covalent bond and two hydrogen bonds that are unsymmetrical.71 However, the proton is not within the plane formed by the oxygen ions, as indicated by the sum of the three ∠O-H-O angles (329° < 360°). In configuration 1B, a relatively rare fourcenter bond is formed, in which the proton is bonded to three acceptor groups, and all three hydrogen-bond angles are greater than 90°.71 Note that the distance of the proton from its third

nearest neighbor is only 2.34 Å, whereas the hydrogen-bond angle is 101.4°. In both 1A and 1B, O1 is forced slightly out of the plane formed by the Zr ions, as indicated by the sum of the three ∠Zr-O1-Zr angles (356°). The three metastable configurations of the proton bonded to O2, labeled 2A, 2B, and 2C, are significantly less stable than either 1A or 1B. In monoclinic zirconia, the 4-coordinated oxygen atom O2 is located within a tetrahedron formed by four Zr ions (Figure 1).16 However, one of the ∠Zr-O2-Zr angles is significantly larger than the tetrahedral angle (131 versus 109.5°). Thus, a metastable configuration 2A, 0.39 eV higher in energy than 1A, can be formed by the displacement of O2, making it 3-coordinated (sum of the three ∠Zr-O2-Zr angles is 359°; Figure 7). In configurations 2B and 2C, the energy gained from O2-H bond formation is offset by the Coulombic repulsion of the proton from the Zr ions (Figure 7). Configurations of the proton initially oriented toward the fourth face of the tetrahedron are observed to relax to either 2A or 2B. Thus, the oxygen atom O2 and the nearest neighbor cations are significantly displaced in the configurations 2A, 2B, and 2C, leading to a greater distortion of the ordered lattice. This can be quantified by computing the root-mean-square displacement of the heavy atoms as a function of their distance from the proton (Figure 8). To test the convergence of results, optimal geometries of the proton in a 3 × 3 × 3 arrangement of monoclinic unit cells are predicted. Models are constructed by adding unit cells to stable configurations of the proton centered in a 2 × 2 × 2 arrangement of monoclinic unit cells. Results obtained are essentially identical to those reported (Table 1), e.g., the lowest energy structure of the proton bonded to O1 is favored by 0.39 eV with respect to that of the proton bonded to O2. Starting from the optimized

8020


Mantz and Gemmen

Figure 8. Root-mean-square displacement, rmsd, of the heavy atoms (from their ordered lattice positions) versus their distance from the proton for every structure identified (Table 1).

TABLE 2: Bader Charges in Units of e of the Atoms in the 96-Atom Supercell Model of Monoclinic Zirconia and in the Structures 1A (Figure 4), 1B, and 2Aa m-ZrO2 1A 1B 2A a

Zr

O1

O2

H

O(H)

2.595 2.602 2.602 2.603

-1.270 -1.265 -1.265 -1.266

-1.325 -1.322 -1.322 -1.322

1.000 1.000 1.000

-1.693 -1.695 -1.720

O(H) is the oxygen atom bonded to the hydrogen atom H.

structures in Table 1, geometry optimizations are performed sampling over four nonsymmetric k points in the irreducible Brillouin zone. Again, results obtained are very similar to those reported. Partial charges of the atoms in the 96-atom supercell model of monoclinic zirconia, as well as in the highly relevant structure 1A, are computed by performing a Bader analysis (Table 2),72 decomposing the charge density using the algorithm developed by Henkelman and co-workers.73-75 For both of these systems, the charge density is computed using a grid of 2003 points. Because the bonding in monoclinic zirconia is largely ionic in character,55 with most of its semicore and valence charge density localized near the nuclei, the addition of a core electronic density is not required to obtain accurate charges.73 Dividing the partial charge of the zirconium atom in monoclinic zirconia (Table 2) by its formal charge, the percent ionic bonding character of monoclinic zirconia is predicted to be 65%, compared to a consensus value of about 80%.55 In the structures 1A, 1B, and 2A, the Bader charges of Zr, O1, and O2 are computed by averaging the values obtained for atoms that are at least 4.5 Å removed from the proton (Table 2). The charges are similar to those in monoclinic zirconia. In these structures, the proton is predicted to have a charge of 1.000e, while the charge of the oxygen atom that is bonded to it is always much more negative than that of the other oxygen atoms. Interestingly, the oxygen ions that are located closer to the proton have relatively larger Bader charges. Specifically, in 1A, the first and second nearest neighbor oxygen ions, O1st and O2nd, have Bader charges that are respectively 0.021e and 0.007e more negative than the appropriate bulk values in Table 2. In 1B, the three nearest neighboring oxygen ions have charges that are 0.010e, 0.009e,

and 0.006e more negative, respectively. Because a correspondence between hydrogen bond length and Bader charge magnitude is predicted, the structures 1A and 1B can be labeled “three-center” bond and “four-center” bond configurations, following the geometric definitions given previously from ref 71. A Mulliken and Bader charge analysis for all three lowpressure zirconia polymorphs, considering convergence with respect to the number of charge density grid points and k points, is provided in the Supporting Information. The stretching and bending motion of the proton in the structure 1A is characterized by performing a vibrational analysis. Harmonic frequencies are calculated by finite differences of first derivatives, using the default step length of 0.02 Bohr. All of the atoms are included in the analysis. A high frequency mode is identified as the stretching of the O-H bond, with a frequency of νs ) 3505 cm-1. The next highest frequency modes are identified as corresponding to bending motions of the proton, with frequencies of νb1 ) 968 cm-1 and νb2 ) 795 cm-1. However, additional low frequency modes are coupled to these latter modes. Because of this coupling, the librational motion of the proton cannot, strictly speaking, be used as done elsewhere to estimate which structures would be favored by entropy.76 The values computed for the structures 1B and 2A are slightly different (Table 1). This is due to the fact that the hydrogen bond length H · · · O1st decreases in the order 1B > 1A > 2A, and (slightly) shorter hydrogen bonds are accompanied by a (slight) red-shifting of the O-H stretch frequency and blueshifting of the O-H wag frequencies.64,67 It is expected that entropy is not strongly favoring any of these structures. Because kT ) 0.09 eV at 1000 K, the binding partner O1 is concluded to be more relevant than O2 at fuel cell operating temperatures. In the following, the effect of entropy on the lowest energy configuration 1A is explored in greater detail. 4. Dynamical Properties of Zirconia-Based Systems In this section, BOMD simulations of monoclinic zirconia and of the configuration 1A are described. First, the structure of monoclinic zirconia is predicted at 1000 K. Next, the dynamics of the proton is characterized at ≈700 K and ≈1000 K. At 690 K, proton transfer is not observed, and the hydrogen bonding and vibrational motion of the proton can be quantified. At 990 K, diffusive jumps of the proton between O1 atoms are observed, allowing a proton transfer barrier height to be estimated. A standard protocol is used to simulate monoclinic zirconia at 1000 K. The 2 × 2 × 2 unit cell model with dimensions of 10.3656 Å × 10.4234 Å × 10.7462 Å and β ) 98.835° is equilibrated for 2 ps by rescaling the temperature, with a tolerance of ( 150 K. Starting from the final configuration obtained, a production run is performed for 8.26 ps, during which the average temperature, with root-mean-square fluctuations, is 1000 ( 60 K. The size of the simulation time step is reduced from 2.056 to 1.028 fs after 0.82 ps of data collection and from 1.028 to 0.411 fs after 2.88 ps. The fluctuations in the conserved (total energy) quantity are decreased by these changes (Figure 9), in this case to a value on the order of 0.1 K/atom. The wave function convergence criterion is kept fixed, 5 × 10-7 Ha, yielding a constant-in-magnitude drift of the total energy equal to a few K/(ps-atom) from 0-7.4 ps (Figure 9), ensuring sufficiently accurate sampling during a 4-5 ps simulation. The role played by these simulation parameters during BOMD simulations of both isolated and periodic systems is examined further in Appendix A and the Supporting Information.



Figure 9. Conserved (total energy) quantity Etot(t) at time t relative to Etot(0) for the 96-atom supercell model of monoclinic zirconia at 1000 K. Points at which the time step is reduced are indicated by arrows.

Figure 11. Bond lengths formed between the proton and nearby oxygen atoms during simulation #1. Each oxygen atom is identified by a unique line color. In monoclinic zirconia (without the proton), the atoms represented by black, blue, and green lines would be O1 atoms, while those represented by the red and yellow lines would be O2 atoms.

Figure 12. Hydrogen-oxygen radial distribution functions during selected time intervals of simulation #1. Figure 10. Radial distribution functions for the 96-atom supercell model of monoclinic zirconia from the first half (represented by solid curves) and second half (represented by dashed curves) of a BOMD simulation at 1000 K. Shown below the top panel are the predicted distances between pairs of atoms in this model system (represented by times, plus and circled dot symbols), with degeneracies of either one or two.

The structure of monoclinic zirconia at 1000 K is characterized by computing the oxygen-oxygen, gO-O(r), zirconiumoxygen, gZr-O(r), and zirconium-zirconium, gZr-Zr(r), radial distribution functions. The computed functions are well converged; those obtained from both halves of the production run are nearly identical (Figure 10, solid and dashed curves). The peaks and troughs are determined by the distribution of atomic distances (Figure 10, symbols). Thus, they would be expected to move inward slightly if the temperature were decreased and the dimensions of the supercell reduced. Consistent with the interpretation of neutron diffraction data,42,43,54 the Zr-O coordination number nZr-O(2.97 Å) is seven, obtained by integration of gZr-O(r) from 0 Å to its first minimum r ) 2.97 Å using

nX-Y(r) ) 4πFY

∫0r r2gX-Y(r) dr

(8)

where FY is the atomic density of atom type Y. Also, nZr-Zr(3.73 Å) ) 6.96 ≈ 7, and nO-O(3.39 Å) ) 9.53 ≈ 9.5. The latter value of 9.5 is the average number of O1-O and O2-O distances that are less than 3.4 Å (Figure 10, blue crosses). Next, two BOMD simulations are performed starting from the lowest energy structure 1A (Figure 4). The simulation time step is 0.484 fs, and the wave function convergence criterion is 10-6 Ha. Simulations #1 and #2 are comprised of an equilibration period followed by a data collection period of length 4.76 and 5.11 ps, during which the average temperature, with root-

mean-square fluctuations, is 690 ( 40 and 990 ( 60 K, respectively. The magnitude of the energy fluctuations and drift are similar to those illustrated in Figure 9, after 2.88 ps. Proton transfer is not observed during simulation #1, of the proton in monoclinic zirconia at 690 K. Rather, the proton is bonded to a single oxygen atom for the duration of the simulation (represented by the black curve in Figure 11). The motion of the proton is generally confined within an irregular tetrahedral volume that is formed by the proton’s favored hydrogen-bonding partner at the apex (represented by the blue curve in Figure 11), and three neighboring oxygen atoms at the triangular base (Figure 11, green, red, and yellow curves). The structure of the oxygen atoms around the proton can be evaluated by computing the hydrogen-oxygen radial distribution function, gH-O(r) (Figure 12, thick black line). Because large displacements of the proton are not observed, reasonably wellconverged results are obtained, as judged by an explicit comparison of the functions computed from both halves of the simulation (Figure 12, thin solid and dashed black lines). The peak at 0.98 Å is intramolecular; the calculated coordination number nH-O(1.29 Å) ) 1.00. The first intermolecular peak at 2.00 Å, providing the most probable hydrogen bond distance, is assigned to nH-O(2.21 Å) - nH-O(1.29 Å) ) 1.25 nearest neighbors, where r ) 2.21 Å is the first minimum of gH-O(r). Assuming that this distance is the only criterion for hydrogen bonding, the proton is interacting with its favored hydrogen bonding partner 87% of the time (Figure 11, blue curve) and with three remaining oxygen atoms 31, 6, and 1% of the time, respectively. The second intermolecular peak of gH-O(r) at 2.73 Å is assigned to nH-O(3.37 Å) - nH-O(2.21 Å) ) 7.11 second nearest neighbors, where r ) 3.37 Å is the second minimum of gH-O(r). During simulation #1, the distribution of the hydrogen-bond angle, R ) ∠O-H · · · O, is sharply peaked at 123° (Figure 13, thick black line). The hydrogen bonds formed are far from linear, in contrast to those in liquid water,78,79 but these can form in

8022


Mantz and Gemmen

Figure 13. Normalized distributions of the angle R obtained during selected time intervals of simulation #1 by considering oxygen atoms within 2.21 Å of the proton.

Figure 14. Evolution of the O-H bond formed and the conserved (total energy) quantity Etot(t) at time t relative to Etot(0) during simulation #1.

Figure 15. Power spectrum computed from the velocity-velocity autocorrelation function of hydrogen extracted from simulation #1.

the solid state:71 the packing of the OH- group in monoclinic zirconia is determined by its shape and a variety of forces, of which hydrogen bonding is only one. Results are similar to those described in a density-functional study of protonated BaZrO3 perovskite oxide.77 The stretching frequency of the proton can also be determined. This is accomplished by computing the fast Fourier transform32 of the O-H bond length during one of several time intervals in which the hydrogen bonds are rather long, and the O-H stretching motion can be described as approximately harmonic (Figure 14).76 The resultant spectrum exhibits a single sharp peak at 3501 cm-1. This result is quite consistent with that obtained from a harmonic vibrational analysis, 3505 cm-1, in which the hydrogen bond length is a comparable 1.91 Å (section 3). Alternatively, the power spectrum can be computed (Figure 15)32,80 using

P(ν) ∝ |H(ν)| 2, H(ν) )

∫-∞∞ e2πiνt〈vH(0)vH(t)〉 dt

(9)

where 〈vH(0)vH(t)〉 is the velocity-velocity autocorrelation function of hydrogen.81 The values of νs, νb1, and νb2 obtained are 3485, 960, and 770 cm-1, which can be compared to 3505, 968, and 795 cm-1 obtained from the harmonic vibrational analysis (section 3). The agreement implies weakly anharmonic vibrational motion, on average, during the simulation. Interestingly, the most rapid fluctuations in the conserved (total energy) quantity have a period that is equal to one-half that of the O-H stretch (Figure 14). Additionally, the restoring force (not shown) acting on the hydrogen atom is observed to

Figure 16. Bonds formed between the proton and selected oxygen atoms during simulation #2. Each oxygen atom is identified by a unique line color consistent with those in Figure 11. In monoclinic zirconia (without the proton), all of these atoms would be O1 atoms.

be maximal at times equal to one-quarter and three-quarters that of the O-H stretching period, as expected for harmonic motion. In a simulation, the degree of energy conservation and the choice of time step is dictated by how rapidly the potential energy changes with distance.81,82 Accordingly, this correlation of energy fluctuations with the force magnitude is physically reasonable. It is further investigated in Appendix A. A second BOMD simulation, performed at 990 K, is more eventful than that at 690 K. Frequent hops of the proton between O1 atoms are observed, initiated by the lattice vibrations, occurring within 0.1 ps, and separated by ≈2 ps (Figure 16). The barrier Eo to proton transfer between O1 atoms can be roughly estimated using the transition-state theory expression83

kTST )

kBT Q‡ kBT exp(-Eo /kBT) ≈ [1 h Q h exp(hυs /kBT)] exp(-Eo /kBT) (10)

where kB is Boltzmann’s constant, h is Planck’s constant, T is the absolute temperature, and Q and Q‡ are the partition functions for the reactant and transition state, respectively. In simplifying eq 10, the approximation is made that the proton stretching frequency νs is the only affected mode during a transition.77,84 Substituting kTST ) 0.5 ps-1, T ) 990 K, and νs ) 3500 cm-1 into eq 10, Eo ) 0.32 eV is obtained. During a dynamical simulation of the system H + m-ZrO2 (q ) 0),16 the same qualitative mechanism of proton transfer between O1 ions is observed as just described for the system H + m-ZrO2 (q ) +1). However, the average delay between jumps is only 0.5 ps at 600 K, versus 2.0 ps at 990 K, and the estimated proton activation energy is ∼0.1 eV versus 0.32 eV. This discrepancy is remarkable, given that the details of the simulations are similar, i.e., both are DFT-GGA computations with a plane-wave basis set and pseudopotentials. It is suggested to be due to the effect of the extra electron present in the other simulated system. As established in section 3, the O-H σ* level is located above the conduction band minimum in stable structures of H + m-ZrO2 (q ) +1, 0, and -1). During a proton transfer event, the O-H bond is stretched and weakened, lowering the O-H σ* level toward the conduction band minimum. If the conduction band minimum is occupied by an extra electron, the O-H bond will be weakened further, decreasing the barrier to proton transfer. This effect is illustrated by examining the electronic density of states of a “transition state” structure extracted from simulation #2 at 0.65 ps, in which the shortest O-H bond distance is 1.19 Å. Compared to those of the lowest energy configuration 1A (Figure 3, black curve), the bands of this structure are shifted to higher energy, denoting

Protonated Forms of Monoclinic Zirconia a weakening of all bonding states, especially those corresponding to O 2s-H 1s and O 2p-H 1s states (Figure 3, gray curve) and, by implication, a conduction band minimum with an enhanced fraction of H 1s character relative to that of 1A. In the extreme case of a hydrogen interstitial, the 1s electron is seen to occupy a localized band gap level located 1.2 eV above the valence band maximum.16 This explanation is not meant to suggest that n-type doping will enhance the proton conductivity of monoclinic zirconia or other wide-gap oxides. In general, proton diffusion is rapid in materials such as the perovskite ABO3 oxides with open, flexible oxygen sublattices that can provide momentarily short O-O distances to permit the breaking of the O-H bond and that are not hindering the rotational diffusion/reorientation of the proton.85-87 Less positively charged dopants with small ionic radii that are substituted for the B-site cations can permit the formation of strong hydrogen bonds and, by enhancing the negative charge of oxygen ions surrounding the dopant, may also increase the O-H bond strength.77,88-90 Thus, the role of dopants is complex and an active area of research. 5. Conclusions In this work, hydrogen in monoclinic zirconia is studied within the framework of DFT, using the HCTH/120 exchangecorrelation functional. To validate the methodology, the Murnaghan equations of state of yttria and the three low-pressure zirconia polymorphs are predicted, for the first time using the HCTH/120 functional. Results are supporting the conclusion from a previous theoretical study that GGA functionals can predict the correct energetic ordering of phases, but that the energy differences between phases are overestimated. The formation energies of a hydrogen defect in monoclinic zirconia placed in different charge states are computed. A proton is predicted to be the favored charge state of a hydrogen atom in monoclinic zirconia for all values of the chemical potential of electrons, under different environmental conditions (i.e., under either hydrogen- or oxygen-rich conditions). This finding is consistent with a previous theoretical prediction at a single value of the chemical potential of electrons as well as muonium spectroscopy experiments and is supporting the theory of hydrogen in wide-gap oxides developed by Peacock and Robertson based on their study of hydrogen in various oxides (not including monoclinic zirconia).19-21 In this theory, the chemical potential of electrons at which positive and negative charge states are equal in energy is predicted to lie at a relatively constant energy above the top of the valence band. The geometries of protonated monoclinic zirconia and binding energies of the proton (as formation energies) are predicted. Different stable and metastable structures of the proton in monoclinic zirconia bonded to both types of oxygen atoms O1 and O2 are described. Two of these minima are associated with O1 and are nearly isoenergetic but involve the formation of either three-center or four-center bonds. The most stable of these structures is favored by at least 0.39 eV compared to configurations of the proton bonded to O2, which is explained by Coulombic repulsion between the proton and zirconium ions. Our findings are unchanged when a larger model system is considered and, when compared to those described previously,16 suggest that our results are not affected by the homogeneous neutralizing background charge density implicitly included, a reasonable finding for a system in the solid state. The most relevant structures are characterized in detail by performing a Bader charge and vibrational analysis. In both monoclinic zirconia and protonated monoclinic zirconia, a

J. Phys. Chem. C, Vol. 114, No. 17, 2010 8023 percent ionic bonding character of 65% between zirconium and oxygen atoms is computed. For the most stable structure identified, the vibrational frequency of the proton is computed from a harmonic (phonon) mode analysis, yielding 3505, 968, and 795 cm-1. In structures with stronger hydrogen bonds between the proton and neighboring oxygen atoms, a red-shifting of the O-H stretch frequency and blue-shifting of the O-H wag frequencies is observed. Lastly, the dynamics of the proton in monoclinic zirconia is explored by performing BOMD simulations at two different temperatures. At 690 K, the motion of the proton is still largely harmonic, based on its computed vibrational frequencies of 3485, 960, and 770 cm-1. At 990 K, a hopping mechanism of proton transfer between O1 ions is observed. However, the rate is slower than previously predicted at 600 K when an extra electron is added to the system.16 It is suggested that, when the O-H bond is stretched/weakened, the conduction band minimum is seen to acquire an appreciable fraction of O-H σ* character. Thus, dissociation of the O-H bond is facilitated by this extra electron, the energy of the transition state is lowered, and a lower transfer barrier is obtained. Acknowledgment. We would like to thank Dr. Jacob Gavartin for providing a reprint of ref 16. This research was supported in part by the National Science Foundation through TeraGrid resources provided by Pittsburgh Supercomputing Center. The Bader analyses were performed using the program bader (v0.26a). The density-of-states plots were generated using the program dosplot. Selected figures were prepared using the VMD Molecular Graphics Viewer.91 Appendix A. Energy Conservation for Different Time Steps To understand the origin of the fluctuations in the conserved (total energy) quantity during simulations of a proton in monoclinic zirconia, BOMD simulations are performed of a single hydrogen molecule in the gas phase at 50 K. The molecule’s translational and rotational degrees of freedom are constrained by fixing five of the six atomic coordinates. To select an initial bond length, a geometry optimization is performed. The result, req ) 0.747 Å (dashed line bisecting Figure 17, top panel), is in accord with spectroscopy, 0.741 Å.92 During a simulation of 0.12 ps, the H-H vibration is largely, but not entirely, harmonic. For harmonic vibrational motion at a given temperature, the maximum bond length minus req will equal req minus the minimum bond length. As can be seen from the evolution of the H-H bond length (Figure 17, top panel, solid line), the H-H bond is slightly more amenable to elongation than to contraction. Furthermore, the magnitude of the force acting on the movable hydrogen atom is larger when the molecule is fully compressed versus extended, e.g., 0.395 eV/Å at 25 fs versus 0.376 eV/Å at 30 fs (Figure 17, bottom panel, solid line). Both of these facts imply weakly anharmonic motion, even at 50 K. It was confirmed that the force magnitude at a given time is precisely equal to the value predicted upon taking the numeric derivative32 of the Kohn-Sham (potential) energy at that particular distance. The evolution of these quantities, the bond length and the force magnitude, is the same for different time steps ranging from 0.024-0.484 fs. For each of the simulations performed with different time steps, the error in the conserved (total energy) quantity is composed of a systematic drift that is linear in time, as well as periodic fluctuations (Figure 17, middle panel, solid lines). This has been previously observed for other systems,27,93 and both

8024


Figure 17. (Top) Bond length evolution during a BOMD simulation of molecular hydrogen initially at 50 K and with constrained degrees of freedom; (middle) the conserved (total energy) quantity Etot(t) at time t relative to Etot(0) for different values of the time step, where the dashed lines are the best linear fits to the curves; (bottom) atomic force magnitude.

errors are dependent on the choice of time step. The magnitude of the energy drift, which is most severe for the smallest time step chosen, can be reduced by a more stringent wave function convergence criterion.27,93 Meanwhile, the increasing amplitude of the energetic fluctuations as a function of the time step can be easily explained. If the potential energy is changing rapidly with distance, and the atomic force magnitude is relatively large, then the error in the conserved (total energy) quantity will be large and more sensitive to the value of the time step. If a simulation is performed at 50 K during which the molecule is allowed to translate and rotate, the net force acting on either hydrogen atom is always the same and, interestingly, observed to be greater when the molecule is fully extended versus compressed, because of the centripetal force. However, the period of the energy oscillations is still identical to that of the force magnitude or half that of the H-H bond length evolution, as illustrated for the fixed molecule (Figure 17) and as similarly observed for protonated monoclinic zirconia (Figure 14). Supporting Information Available: The predicted equations of state of yttria and the three low-pressure zirconia polymorphs, with detailed comparison to previous results; predicted Bader and Mulliken atomic charges for these zirconia polymorphs; minimum energy structures of the proton in monoclinic zirconia; an examination of energy conservation during BOMD simulations of a simple periodic system, eight silicon atoms; references to studies of commonly used integrators. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Bessler, W. G.; Warnatz, J.; Goodwin, D. G. Solid State Ionics 2007, 177, 3371–3383. (2) Horita, T.; Kishimoto, H.; Yamaji, K.; Xiong, Y.; Sakai, N.; Brito, M. E.; Yokokawa, H. Solid State Ionics 2006, 177, 1941–1948. (3) Sun, C.; Stimming, U. J. Power Sources 2007, 171, 247–260. (4) Anderson, A. B.; Vayner, E. Solid State Ionics 2006, 177, 1355– 1359. (5) Shishkin, M.; Zeigler, T. J. Phys. Chem. C 2008, 112, 19662–19669. (6) Shishkin, M.; Zeigler, T. J. Phys. Chem. C 2009, 113, 21667–21678.

Mantz and Gemmen (7) Hansen, K. V.; Norrman, K.; Mogensen, M. J. Electrochem. Soc. 2004, 151, A1436–A1444. (8) Mogensen, M.; Jensen, K. V.; J gensen, M. J.; Primdahl, S. Solid State Ionics 2002, 150, 123–129. (9) Primdahl, S.; Mogensen, M. J. Electrochem. Soc. 1997, 144, 3409– 3419. (10) Holtappels, P.; Vinke, I. C.; de Haart, L. G. J.; Stimming, U. J. Electrochem. Soc. 1999, 146, 2976–2982. (11) Bieberle, A.; Gauckler, L. J. Solid State Ionics 2002, 146, 23–41. (12) Vogler, M.; Bieberle-Hu¨tter, A.; Gauckler, L.; Warnatz, J.; Bessler, W. G. J. Electrochem. Soc. 2009, 156, B663–B672. (13) Goodwin, D. G.; Zhu, H.; Colclasure, A. M.; Kee, R. J. J. Electrochem. Soc. 2009, 156, B1004–B1021. (14) de Ridder, M.; van Welzenis, R. G.; Brongersma, H. H.; Kreissig, U. Solid State Ionics 2003, 158, 67–77. (15) Yamanaka, S.; Nishizaki, T.; Uno, M.; Katsura, M. J. Alloys Comp. 1999, 293-295, 38–41. (16) Shluger, A. L.; Foster, A. S.; Gavartin, J. L.; Sushko, P. V. Models of defects in wide-gap oxides: Perspectives and challenges. In Nano and Giga Challenges in Microelectronics; Greer, J., Korkin, A., Labanowski, J., Eds.; Elsevier: Amsterdam, 2003. (17) Van de Walle, C. G.; Neugebauer, J. Nature 2003, 423, 626–628. (18) Kılıç, C.; Zunger, A. Appl. Phys. Lett. 2002, 81, 73–75. (19) Peacock, P. W.; Robertson, J. Appl. Phys. Lett. 2003, 83, 2025– 2027. (20) Robertson, J.; Peacock, P. W. Thin Solid Films 2003, 445, 155– 160. (21) Xiong, K.; Robertson, J.; Clark, S. J. J. Appl. Phys. 2007, 102, 083710. (22) Cox, S. F. J.; Gavartin, J. L.; Gil, J. M.; Vilaõ, R. C.; Lord, J. S.; Davis, E. A. Physica B 2006, 376-377, 385–388. (23) Cox, S. F. J.; Gavartin, J. L.; Lord, J. S.; Cottrell, S. P.; Gil, J. M.; Alberto, H. V.; Piroto Duarte, J.; Vilaõ, R. C.; Ayres de Campos, N.; Keeble, D. J.; Davis, E. A.; Charlton, M.; van der Werf, D. P. J. Phys.: Condens. Matter 2006, 18, 1079–1119. (24) Wagner, C. Ber. Bunsenges. Phys. Chem. 1968, 72, 778–781. (25) Bay, L.; Horita, T.; Sakai, N.; Ishikawa, M.; Yamaji, K.; Yokokawa, H. Solid State Ionics 1998, 113, 363–367. (26) Norby, T.; Widerøe, M.; Glo¨ckner, R.; Larring, Y. Dalton Trans. 2004, 3012–3018. (27) Marx, D.; Hutter, J. Ab Initio Molecular Dynamics: Theory and Implementation. In Modern Methods and Algorithms of Quantum Chemistry; Grotendorst, J., Ed.; John von Neumann Institute for Computing: Jülich, 2000. (28) CPMD, Copyright IBM Corp. 1990-2003, Copyright MPI fu¨r Festko¨rperforschung Stuttgart 1997-2001. (29) Boese, A. D.; Doltsinos, N. L.; Handy, N. C.; Sprik, M. J. Chem. Phys. 2000, 112, 1670–1678. (30) Troullier, N.; Martins, J. L. Phys. ReV. B 1991, 43, 1993–2006. (31) Hutter, J.; Lu¨thi, H. P.; Parrinello, M. Comput. Mater. Sci. 1994, 2, 244–248. (32) Press, W. H.; Teukolsky, A. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes, 2nd ed.; Cambridge University Press: Cambridge, 1994. (33) Monkhorst, H. J.; Pack, J. D. Phys. ReV. B 1976, 13, 5188. (34) Kresse, G.; Hafner, J. Phys. ReV. B 1993, 47, 558. (35) Kresse, G.; Hafner, J. Phys. ReV. B 1994, 49, 14251. (36) Kresse, G.; Furthmu¨ller, J. Comput. Mater. Sci. 1996, 6, 15. (37) Kresse, G.; Furthmu¨ller, J. Phys. ReV. B 1996, 54, 11169. (38) Perdew, J.; Chevary, J.; Vosko, S.; Jackson, K.; Pederson, M.; Singh, D.; Fiolhais, C. Phys. ReV. B 1992, 46, 6671. (39) Perdew, J.; Chevary, J.; Vosko, S.; Jackson, K.; Pederson, M.; Singh, D.; Fiolhais, C. Phys. ReV. B 1993, 48, 4978. (40) Blo¨chl, P. E. Phys. ReV. B 1994, 50, 17953. (41) Kresse, G.; Joubert, D. Phys. ReV. B 1999, 59, 1758. (42) Howard, C. J.; Hill, R. J.; Reichert, B. E. Acta Crystallogr. 1988, B44, 116–120. (43) Leger, J. M.; Tomaszewski, P. E.; Atouf, A.; Pereira, A. S. Phys. ReV. B 1993, 47, 14075–14083. (44) VandeVondele, J.; Mohamed, F.; Krack, M.; Hutter, J.; Sprik, M.; Parrinello, M. J. Chem. Phys. 2005, 122, 014515. (45) Nevitt, M. V.; Chan, S. K.; Liu, J. Z.; Grimsditch, M. H.; Fang, Y. Physica B 1988, 150, 230. (46) Murnaghan, F. D. Proc. Natl. Acad. Sci. U.S.A. 1944, 30, 244– 247. (47) Jomard, G.; Petit, T.; Pasturel, A.; Magaud, L.; Kresse, G.; Hafner, J. Phys. ReV. B 1999, 59, 4044. (48) Terki, R.; Bertrand, G.; Aourag, H.; Coddet, C. Mater. Sci. Semi. Proc. 2006, 9, 1006–1013. (49) Kra´lik, B.; Chang, E. K.; Louie, S. G. Phys. ReV. B 1998, 57, 7027. (50) Christensen, A.; Carter, E. A. Phys. ReV. B 1998, 58, 8050–8064. (51) Stapper, G.; Bernasconi, M.; Nicoloso, N.; Parrinello, M. Phys. ReV. B 1999, 59, 797–810. (52) Zhao, X.; Vanderbilt, D. Phys. ReV. B 2002, 65, 075105.

Protonated Forms of Monoclinic Zirconia (53) Ackermann, R.; Rauh, E. G.; Alexander, C. A. High Temp. Sci. 1975, 7, 304. (54) Boysen, H.; Frey, F.; Vogt, T. Acta Crystallogr. 1991, B47, 881– 886. (55) Soriano, L.; Abbate, M.; Faber, J.; Morant, C.; Sanz, J. M. Solid State Commun. 1995, 93, 659–665. (56) Blo¨chl, P. E. J. Chem. Phys. 1995, 103, 7422–7428. (57) Makov, G.; Payne, M. C. Phys. ReV. B 1995, 51, 4014–4022. (58) McComb, D. W. Phys. ReV. B 1996, 54, 7094–7102. (59) Dash, L. K.; Vast, N.; Baranek, P.; Cheynet, M.-C.; Reining, L. Phys. ReV. B 2004, 70, 245116. (60) Medvedeva, J. E.; Freeman, A. J.; Geller, C. B.; Rishel, D. M. Phys. ReV. B 2007, 76, 235115. (61) Van de Walle, C. G. Phys. ReV. Lett. 2000, 85, 1012–1015. (62) Zhang, S. B. J. Phys.: Condens. Matter 2002, 14, R881–R903. (63) Yoshino, M.; Liu, Y.; Tatsumi, K.; Tanaka, I.; Morinaga, M.; Adachi, H. Mater. Trans. 2002, 43, 1444–1450. (64) Bjo¨rketun, M. E.; Sundell, P. G.; Wahnstro¨m, G. Faraday Discuss. 2007, 134, 247–265. (65) Ashcroft, N. W.; Mermin, N. D. Solid State Physics; Brooks/Cole Thomson Learning: Pacific Grove, CA, 1976. (66) Pantelides, S. T.; Mickish, D. J.; Kunz, A. B. Phys. ReV. B 1974, 10, 5203–5212. (67) Noguera, C. Physics and Chemistry at Oxide Surfaces; Cambridge University Press: Cambridge, 1996. (68) Rak, Z. S.; Mahanti, S. D.; Mandal, K. C. J. Electron. Mater. 2009, 38, 1539–1547. (69) Tuckerman, M.; Laasonen, K.; Sprik, M.; Parrinello, M. J. Chem. Phys. 1995, 103, 150–161. (70) Todorova, T.; Hu¨nenberger, P. H.; Hutter, J. J. Chem. Theory Comput. 2008, 4, 779–789. (71) Jeffrey, G. A. An Introduction to Hydrogen Bonding; Oxford University Press: New York, 1997. (72) Bader, R. F. W. Atoms in Molecules - A Quantum Theory; Oxford University Press: Oxford, 1990. (73) Henkelman, G.; Arnaldsson, A.; Jońsson, H. Comput. Mater. Sci. 2006, 36, 354–360.

J. Phys. Chem. C, Vol. 114, No. 17, 2010 8025 (74) Sanville, E.; Kenny, S. D.; Smith, R.; Henkelman, G. J. Comput. Chem. 2007, 28, 899–908. (75) Tang, W.; Sanville, E.; Henkelman, G. J. Phys.: Condens. Matter 2009, 21, 084204. (76) Haiber, M.; Ballone, P.; Parrinello, M. Am. Mineral. 1997, 82, 913– 922. (77) Bjo¨rketun, M. E.; Sundell, P. G.; Wahnstro¨m, G. Phys. ReV. B 2007, 76, 054307. (78) Mantz, Y. A.; Chen, B.; Martyna, G. J. Phys. Chem. B 2006, 110, 3540–3554. (79) Lee, H.-S.; Tuckerman, M. E. J. Phys. Chem. C 2008, 112, 9917– 9930. (80) Lobaugh, J.; Voth, G. A. J. Chem. Phys. 1997, 106, 2400–2410. (81) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1987. (82) Tuckerman, M.; Berne, B. J.; Martyna, G. J. J. Chem. Phys. 1992, 97, 1990–2001. (83) Steinfeld, J. I.; Francisco, J. S.; Hase, W. L. Chemical Kinetics and Dynamics, 2nd ed.; Prentice Hall: Hoboken, NJ, 1999. (84) Bjo¨rketun, M. E.; Sundell, P. G.; Wahnstro¨m, G.; Engberg, D. Solid State Ionics 2005, 176, 3035–3040. (85) Kreuer, K. D.; Mu¨nch, W.; Traub, U.; Maier, J. Ber. Bunsenges. Phys. Chem. 1998, 102, 552–559. (86) Kreuer, K. D. Solid State Ionics 2000, 136-137, 149–160. (87) Kreuer, K. D. Ann. ReV. Mater. Res. 2003, 33, 333–359. (88) Kreuer, K. D.; Mu¨nch, W.; Ise, M.; He, T.; Fuchs, A.; Traub, U.; Maier, J. Ber. Bunsenges. Phys. Chem. 1997, 101, 1344–1350. (89) Giannici, F.; Longo, A.; Balerna, A.; Kreuer, K. D.; Martorana, A. Chem. Mater. 2009, 21, 2641–2649. (90) Merinov, B.; Goddard, W., III J. Chem. Phys. 2009, 130, 194707. (91) Humphrey, W.; Dalke, A.; Schulten, K. J. Mol. Graph. 1996, 14, 33–38. (92) Huber, K. P.; Herzberg, G. Constants of Diatomic Molecules; Van Nostrand Reinhold: New York, 1979. (93) Herbert, J. M.; Head-Gordon, M. Phys. Chem. Chem. Phys. 2005, 7, 3269–3275.

JP810601J

Protonated Forms of Monoclinic Zirconia: A Theoretical Study - The

Recommend Documents