Initial Reduction of the PdO(101) Surface: Role of Oxygen Vacancy

Subscriber access provided by UNIV OF NEW ENGLAND

C: Surfaces, Interfaces, Porous Materials, and Catalysis

Initial Reduction of the PdO(101) Surface: Role of Oxygen Vacancy Formation Kinetics Minkyu Kim, Li Pan, Jason F. Weaver, and Aravind Asthagiri J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.8b08226 • Publication Date (Web): 19 Oct 2018 Downloaded from http://pubs.acs.org on October 21, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Initial Reduction of the PdO(101) Surface: Role of Oxygen Vacancy Formation Kinetics Minkyu Kim†,∥, Li Pan†,∥, Jason F. Weaver‡, and Aravind Asthagiri†* †William

‡

G. Lowrie Department of Chemical and Biomolecular Engineering, The Ohio State University, Columbus, OH 43210

Department of Chemical Engineering, University of Florida, Gainesville, FL 32611, USA

*To whom correspondence should be addressed ([email protected])

∥These authors (M.K. and L.P.) contributed equally to this work.

1 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract We investigated the initial thermal reduction of the PdO(101) surface using density functional theory (DFT) calculations. The defect-free PdO(101) surface consists of parallel rows of undercoordinated Pd and O atoms (Pdcus and Ocus, respectively). DFT was used to map out the elementary processes of Ocus vacancy formation, oxygen surface atom diffusion along the Pdcus rows, and various O2 formation pathways. Because oxygen vacancies occur during the reduction of the surface, the elementary processes were examined with and without the presence of adjacent Ocus vacancies. DFT calculations show that the presence of oxygen vacancies strongly affects the barriers of adjacent surface processes. Barriers for O vacancy formation are reduced by nearly 50% with an adjacent O vacancy. The barrier for O2 formation along the Pdcus row drops from 1.54 eV on the defect-free PdO(101) surface to 0.81 eV in the presence of a Ocus vacancy dimer. However, the presence of Ocus vacancies also leads to more strongly bound O2 on the Pdcus row. The strong neighbor effect of O vacancies on the PdO(101) surface will favor the growth of O vacancy chains and increase the rate of thermal reduction of the surface, both features that have been observed in recent ultrahigh vacuum isothermal reduction experiments of the PdO(101) surface.


Page 2 of 35

Page 3 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1. Introduction Palladium based catalysts have been widely used for oxidation catalysis such as methane oxidation in gas turbines and CO oxidation in catalytic car converters. Despite extensive surface science studies under both ultra-high vacuum (UHV) conditions (Ptotal: 10-10 Torr) and ambient pressures (Ptotal: 10-2 – 102 Torr), there is still ongoing debate on the active phases for Pd single crystal surfaces under oxygen-rich conditions

1–4.

For CO oxidation, experiments on low-Miller

index single crystal Pd surfaces under ambient pressures based on polarization modulation infrared reflection adsorption spectroscopy claimed that metallic Pd is more active toward CO oxidation than Pd oxide phases5–9. However, near-ambient pressure X-ray photoelectron spectroscopy and surface X-ray diffraction experiments on Pd(100) and Pd(111) under similar conditions1,10–15 reported that Pd oxide phases are responsible for the high CO oxidation reactivity. Studies based on UHV conditions also reported facile pathways for CO oxidation on bulk PdO(101)2,16,17. Recently, UHV studies varying CO/O2 ratio at low pressure (Pi ~ 10-9 Torr) and pseudo-steady state conditions proposed that both metallic and oxide phases are active toward CO oxidation, but that metallic Pd shows 2 to 3 times greater reactivity than the oxide phases18. There are several challenges to unambiguously characterize the reactivity of surface phases, including the coexistence of phases and mass transfer effects. Recent planar laser-induced fluorescence (PLIF) experiments for CO oxidation on Pd(100) have shown that the near surface gas concentration in these experiments can be dramatically different than the operating partial pressures and can couple to changes in the surface phases19–21. Furthermore, experiments have also demonstrated that mass transfer limitations can occur where the reactions are faster on the metal and metal oxide surfaces than CO diffusion to the surface. Under this scenario it is nontrivial to resolve differences in the reactivity of the metal versus the oxide phases. 3 ACS Paragon Plus Environment


First principles based kinetic Monte Carlo (kMC) and micro kinetic modeling simulations can be important tools for directly exploring the reactivity-surface phase relationship, since these methods connect elementary steps on the surface to catalyst reactivity. The challenge in using these methods is the uncertainty that all relevant processes are captured (both that the rates are accurate and that important effects/processes are not neglected). This challenge is particularly acute for oxidation reactions since explicitly following large scale changes involved in surface phase transitions is difficult. A common approach to this challenge is to characterize the reactivity of phases separately and to couple this information with ab initio thermodynamics and experimental studies to develop a holistic picture of the relationship between surface oxygen phase and reactivity. First principles kMC simulations have been used to study CO oxidation on Pd(100) and Pd(111) surfaces22–26, RuO2(110)27–29, and single layer (( 5 × 5)𝑅27°) surface oxide on Pd(100)27,30,31. In the work on the surface oxide on Pd(100), a multi-lattice kMC framework was used to allow for the transition between Pd(100) and the surface oxide phase locally. These studies indicate that high CO reactivity on Pd(100) is associated with the coexistence of the metal and surface oxide phase. Grönbeck et al. have developed microkinetic models for methane oxidation on Pd(111), Pd(100), and PdO(101) surfaces based on mapping out the pathways for complete oxidation of methane using density functional theory (DFT)32,33. Duan and Henkelman used DFT in combination with simple microkinetic models to propose that the surface Pd5O4 phase on Pd(111) was active for CO oxidation through a Eley-Rideal reaction mechanism34. Although these modeling studies provide important insights into the differences in reactivity of the oxide versus


Page 4 of 35

Page 5 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


metal phase, prior studies of the oxide phase did not consider neighbor effects due to the presence of oxygen vacancies. Under such an assumption, the elementary rates calculated with DFT are assumed to not change in the presence of adjacent oxygen vacancies. While such an assumption may be valid for many oxide surfaces where the O vacancies do not directly interact with the undercoordinated metal site (e.g. RuO2(110)), on the PdO(101) surface the coordinatively undersaturated (cus) Pd surface atom is bonded directly to the coordinatively undersaturated (cus) O surface atom. Therefore, O vacancies on the PdO(101) surface have potential to effect elementary steps occurring on the adjacent Pd site, but to our knowledge this effect has not been extensively studied33,35. In this study, we focus on the initial stages of the thermal reduction of PdO(101) using density functional theory (DFT). Our goal is to explore the effect of O vacancies on the surrounding process related to PdO(101) surface reduction. Experimentally, Hinojosa and Weaver studied the thermal reduction of PdO(101) using scanning tunneling microscopy (STM) and measurements of the isothermal decomposition rate36. Isothermal reduction of PdO(101) was performed at 720 K, a temperature at the onset of PdO(101) film reduction as observed in temperature programmed desorption (TPD) experiments. The isothermal decomposition rate exhibits several regimes as the surface transforms from a PdO(101) film to a Pd(111) surface. In the early stages of reduction, STM images show the presence of oxygen vacancy chains along the Pd rows on the PdO(101) surface. Corresponding to the presence of these oxygen vacancy chains they observe an increase in O2 desorption rate with time (i.e. O2 desorption rate accelerates). As reduction continues the STM images show co-existing domains and eventually as the metal domain phase grows they observe a drop in the O2 desorption rate. The regime where the O2 desorption rate accelerates can



be defined as an autocatalytic regime since the rate of reduction of the PdO(101) surface increases with increasing reduction events. While the STM studies provide us snapshots of the surface morphology during the initial stages of reduction, they do not provide information on the atomic-level processes that dictate the observed surface evolution. The evolution of PdO(101) decomposition occurs through several distinct atomic-scale processes, such as the hopping of lattice oxygen (i.e. initial oxygen vacancy formation) and the formation and subsequent desorption of O2 molecules. We have applied DFT to examine the energetics of pathways for possible surface processes relevant to the thermal reduction of the PdO(101) surface. As a first step we have focused on processes related to the formation and evolution of oxygen vacancies in the initial stages of PdO(101) surface reduction. By constraining our study to the initial stages of reduction, we avoid processes that may require complex rearrangement of the surface (e.g. small scale rearrangement of the Pd metal atoms on the reduced surface or the evolution of co-existing metal and metal oxide phases). While this choice is partly driven by practical considerations of applying DFT to these more complex processes, we expect that under oxygen-rich conditions the surface would primarily contain surface O vacancies and not undergo extensive reduction. In this present study, we report DFT derived barriers for the elementary processes involved in the thermal reduction of PdO(101). These processes include initial vacancy formation, diffusion of the oxygen adatom, and O2 formation and desorption pathways. The unique aspect of this study is that these pathways are all examined with and without the presence of O vacancies. The DFT calculations show a large neighbor effect associated with O vacancies. Both O vacancy formation and O2 formation barriers are reduced when adjacent O vacancies are present. Based on the pathways examined we find that O vacancy chain formation is favored over random distribution 6 ACS Paragon Plus Environment

Page 6 of 35

Page 7 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


of O vacancies. Furthermore, the O vacancy chains are associated with a more facile sequence of pathways to O2 formation and desorption, which may be associated with the experimentally observed acceleration of thermal reduction with the formation of O vacancy chains. We also find that the Heyd-Scuseria-Ernzerhof (HSE06) functional37,38, which has been shown to provide better agreement with experimental measurements of desorption energies and adsorption sites for CO and O2 on the PdO(101) surface17,39, dramatically overestimates the barriers for O vacancy formation due to a destabilization of O atoms on the Pdcus rows. The PBE functional is more appropriate for most of the processes except for O2 desorption which is more accurately captured by the HSE06 functional39. Overall, if the appropriate combination of PBE and HSE06 DFT barriers are used, we obtain a qualitative explanation of the proposed autocatalytic regime. Since oxidation catalysis over transition-metal oxides typically follows a Mars-van Krevelen mechanism40,41, the strong neighbor effect of the oxygen vacancies reported in this study may potentially play an important role in mediating both CO and alkane oxidation on PdO(101). 2. Computational Details The periodic plane wave DFT calculations reported in this paper were performed using the Vienna ab initio simulation package (VASP)42,43 with projector augmented wave (PAW)44 pseudopotentials provided in the VASP database. The general gradient approximation (GGA) with the Perdew-Burke-Ernzerhof (PBE) form was used for the exchange-correlation functional45 with a plane wave cutoff of 400 eV. We performed test calculations with spin-polarization but the effects on the energetics of surface processes were negligible (< 0.01 eV) (see Sec. S1 in SI). Thus, all the results, except for the O2 desorption processes, reported in this paper are from non-spin polarized calculations. Figure 1 illustrates the stoichiometric PdO(101) surface that is investigated in this study. Bulk crystalline PdO has a tetragonal unit cell and consists of square planar units of 7 ACS Paragon Plus Environment


Pd atoms four-fold coordinated with oxygen atoms. The bulk-terminated PdO(101) surface is defined by a rectangular unit cell, where the a and b lattice vectors coincide with the [010] and [1 01] directions of the PdO crystal, respectively. The stoichiometric PdO(101) surface consists of four alternating rows of threefold Pd, threefold O, fourfold Pd, and fourfold O atoms that run parallel to the a direction shown in Figure 1. Thus, half of the surface O and Pd atoms are coordinatively unsaturated (cus).

Figure 1. Top and side views of the PdO(101) surface. Rows of coordinatively unsaturated (cus) and fourfold-coordinated (4f) Pd or O atoms are indicated. The vertical and horizontal arrows a and b represent the [010] and [101] crystallographic directions of PdO.

Unless otherwise noted, the PdO(101) surface was modeled by a rectangular 4 × 1 unit cell with a corresponding 4 × 2 × 1 Monkhorst-Pack k-point mesh46. For the other supercells, the kpoint mesh is scaled accordingly. As in our prior studies47–49 the PdO(101) film was strained (a = 3.057 Å, b = 6.352 Å) to match the PdO(101) film structure resolved by Kan and Weaver50,51. The PdO(101) slab was represented by four layers resulting in a 9 Å thick slab. The bottom layer is fixed, but other lattice atoms are allowed to relax until the forces are less than 0.03 eV/Å. As in 8 ACS Paragon Plus Environment

Page 8 of 35

Page 9 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


our previous work the underlying Pd(111) surface is not included since this would require a large unit cell due to the registry of the PdO(101) film with the Pd(111) surface. We use a vacuum spacing of 20 Å, which is sufficient to eliminate spurious interactions in the surface normal direction. We determined the barriers and pathway related to the reduction of PdO(101) surface using the climbing nudged elastic band (cNEB) method52. The zero point correction (ZPC) reduces the O2 desorption energy by about 0.06 eV, but subsequent application of the ZPC to other select surface processes results in changes less than 0.03 eV. Because of the negligible effects of ZPC on the non-desorption elementary processes, we report non-zero point corrected energies of all the processes except for the O2 desorption processes. GGA-DFT is known to have limitations with the band gap of transition-metal (TM) oxides because of the self-interaction error, but the incorporation of some fraction of exact exchange through the use of hybrid functionals has been shown to be successful for several TM oxides53–56. Experiments find that PdO is a semiconductor with a small band gap varying from 0 to 2.67 eV using different measurements57–60, whereas DFT-PBE predicts it is a metal54,61,62. The HSE06 functional predicts a band gap of 0.8 eV and a oxidation energy of ~100 kJ/mol for bulk PdO, both of which are in good agreement with experimental results39,54. In our previous studies, we found that PBE overestimates the adsorption energy of O2 on the pristine PdO(101) surface, but the HSE06 results are in a good agreement with desorption energies determined from experimental TPD data39. We have also shown, that unlike the PBE functional, the HSE06 functional is able to capture both the adsorption energy and the preferred site for CO adsorption on PdO(101) in comparison to RAIRS experiments17,63. In addition, PBE underestimates the Pd 3d core level shift (CLS) by ~ 1.0 eV in comparison to high resolution core-level spectroscopy experiments of bulk PdO(101), while the PBE0 hybrid functional shows better agreement (~ 0.3 eV overestimation)64.



Page 10 of 35

It should be noted that HSE06 functional is not universally superior to the PBE functional for surface species on PdO(101); good agreement between PBE-derived adsorption sites and energies are found for H2, H2O, CO2, and alkane adsorption on the PdO(101) surface65–68. Based on the above discussion, we have examined the accuracy of our PBE results in this paper by performing parallel HSE06 calculations. However, because of the much larger computational cost associated with using the hybrid functional, we performed single-point HSE06 calculations on relaxed DFT-PBE configurations and for NEB barriers these single-point calculations are performed on the initial, transition, and final states to get an approximation of the HSE06 barriers. We have performed select tests with fully relaxed HSE06 calculations and these comparisons are discussed in more detail below. 3. Results and Discussion 3.1 Energetics of O vacancy formation To determine the preferred type of oxygen vacancies (Ocus vs. O4f) and the strength of the lateral interactions between oxygen vacancies, we initially examine the vacancy formation energy. The oxygen vacancy formation energy is defined as follows: 1

𝐸𝑓,𝑂𝑣 =

𝐸𝑁𝑂𝑣 + 2(𝑁𝐸𝑂2) ― 𝐸𝑠𝑙𝑎𝑏 𝑁

,

(1)

where 𝐸𝑓,𝑂𝑣 is given in units of eV/O atom and 𝐸𝑁𝑂𝑣 denotes the energy of the 4 layer PdO(101) surface with N oxygen vacancies, 𝐸𝑂2, the energy of one isolated O2 molecule in the triplet ground state, and 𝐸𝑠𝑙𝑎𝑏, the energy of the 4 layer pristine PdO(101) surface. By this definition, the oxygen vacancy formation energy is a positive value and larger values indicate that the configuration is less favorable because more energy is needed to remove the oxygen from the surface. We also 10 ACS Paragon Plus Environment

Page 11 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


found negligible change in the oxygen vacancy formation energy when spin-polarization is included for the surface. The lateral interactions between oxygen vacancies are determined by removing one oxygen atom (Ocus or O4f) from the PdO(101) surfaces using unit cells (1 × 1, 2 × 1, 4 × 1, and 8 × 1) that increase in size along the direction parallel to the Ocus rows (𝑎 in Figure 1). Then the same calculation was performed along the direction normal to the Ocus rows (𝑏 in Figure 1) to study the across-row interactions (on the 2 × 2, 2 × 4, and 4 × 2 supercells). Table 1 shows the results of these calculations using both PBE and single-point HSE06 DFT calculations. Table 1. The oxygen vacancy formation energy (in eV/O atom) for Ocus and O4f in different unit cells. Single-point DFT-HSE06 results are given in parenthesis. Unit Cell

Ocus O4f

1×1 2.04 (2.27) 2.23 (2.77)

Direction 𝒂 (in-row) 2×1 4×1 2.01 (2.26) 1.98 (2.27) 2.35 (2.84) 2.36 (2.86)

8×1 1.97 2.36

Direction 𝒃 (across-row) 2×2 2×4 4×2 2.15 (2.40) 2.14 2.13 2.38 (2.83) 2.36 2.39

From Table 1, the formation of Ocus vacancy is more favored than that of O4f by 0.2 - 0.4 eV for different supercells suggesting that oxygen vacancies generated in the initial stage of PdO(101) reduction will be predominantly Ocus vacancies. When increasing the supercell size along the 𝑎 direction, 𝐸𝑓,𝑂𝑣 for Ocus (O4f) slightly decreases (increases), indicating that the in-row interaction for Ocus vacancies is repulsive but attractive for O4f. However, these interactions are weak (~ 0.03 eV/O atom) and we would not expect a strong thermodynamic preference for oxygen vacancy clustering. To test this explicitly we removed two Ocus atoms from the 4×1 and 8×1 unit cells. In one configuration, the two vacancies are placed adjacently, and in the second configuration the two vacancies have the maximum separation allowed in that unit cell. The oxygen vacancy formation energy varies from 1.98 to 2.01 eV from the adjacent versus maximally-separated O vacancy configuration, which is consistent with the results above. We also calculated the across-



row oxygen vacancy interaction by increasing supercell sizes along the direction 𝑏. As shown in Table 1 the vacancy formation energy increases by about 0.1 eV for Ocus vacancies when we go from a (2×1) to an (2×2) unit cell. This weakly attractive interaction is essentially negligible when we move from a (2×2) to an (4×2) unit cell, indicating that this relatively small interaction occurs only over adjacent parallel oxygen rows. The 0.1 eV difference is about 5% of the energy to form a Ocus vacancy and is not sufficiently large to have an effect on oxygen vacancy arrangement on the PdO(101) surface at higher temperatures. We also evaluated the vacancy formation energy by performing single-point DFT-HSE06 calculations for Ocus and O4f vacancies on the 1×1, 2×1, 2×2, and 4×1 unit cells. As noted above, single-point HSE06 calculations are computationally more tractable than full relaxation, but to check the effect of this approximation we performed a full relaxation HSE06 calculation for an Ocus and O4f vacancy in the 4×1 unit cell. The fully relaxed (single-point) HSE06 calculation gives a Ocus vacancy formation energy of 2.25 (2.27) eV/O atom and a O4f vacancy formation energy of 2.85 (2.86) eV/O atom, a difference of less than 0.02 eV. This result implies that the results from single-point DFT-HSE06 can be reliable for this system. Overall, the single-point DFT-HSE06 results lead to the same conclusions as obtained from DFT-PBE. With HSE06 the Ocus vacancy is favored over O4f by 0.4 to 0.6 eV, which is about 0.2 eV more stable than the PBE results; however, the interactions along and perpendicular to the oxygen rows are essentially the same as found using PBE. From the calculated vacancy formation energies, one would expect to observe a random distribution of Ocus vacancies on the PdO(101) surface. This conclusion conflicts with the observed STM images of clustering of oxygen vacancies along the Ocus rows36. One potential cause of the observed clustering may be due to the kinetics of O vacancy formation. In the next section, we 12 ACS Paragon Plus Environment

Page 12 of 35

Page 13 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


present our results on the detailed kinetic mechanism for O vacancy formation on the PdO(101) surface. 3.2 Kinetics of vacancy formation To investigate the detailed kinetics of O vacancy formation, we started with a pristine PdO(101) surface (Figure 1) and examined a series of processes starting with the initial O vacancy formation. Based on our finding of Ocus vacancies being more favorable than O4f vacancies, we constrain our focus to vacancies from the Ocus row. We will refer to an oxygen atom that moves from the Ocus row onto the Pdcus row as O*. After O vacancy formation, the O* atom may follow several different paths including immediate healing of the vacancy, bonding with an adjacent Ocus atom to form O2* molecule that then desorbs, or diffusion away from the vacancy and eventually meeting another O* atom on the Pdcus row to form a O2* molecule that then desorbs. As these processes (in particular O2 desorption) occur, the PdO surface will become populated with an increasing number of O vacancies and each of the above processes may be affected by the presence of O vacancy clusters. We have used NEB calculations to map out all of the above processes, including modifications due to the presence of O vacancies, and the results are described below. For each process, we first make conclusions based on DFT-PBE results and then compare with results from single-point DFT-HSE06 calculations. Unless otherwise noted, all the results discussed in this section are obtained using a 4 × 1 supercell. We have tested the barriers of oxygen vacancy formation with and without the presence of adjacent oxygen vacancies using a 4 × 2 supercell and the difference in barriers is found to be sufficiently small (below 0.1 eV) that it will not affect our conclusions.



Initial vacancy formation To map out the mechanism and barrier for initial O vacancy formation on a pristine PdO(101) surface, we evaluated the energetics of O* being placed on top of and in between the Pdcus or Pd4f atoms adjacent to the Ov site. The O* atoms prefers to reside between two Pdcus atoms (See Figure 2(FS)) and is more stable by 0.3 eV than the second stable site (on top of Pdcus atom). Using these plausible initial (pristine PdO(101)) and final (O* sitting in between Pdcus atoms) states, the barrier to the initial O vacancy formation was calculated using NEB and the pathway is shown in Figure 2. The energy barrier for the initial vacancy formation is 1.71 eV. This high barrier is consistent with TPD experiments which show that the O2 desorption peak starts around 700 K36, which can be associated with Ov formation on the surface. On the other hand, the energy barrier for the reverse process, the healing of the O vacancy, is only 0.47 eV. This finding implies that a single O vacancy is readily refilled unless another O atom can be provided to form an O2 molecule, which would subsequently desorb leaving behind two O vacancy sites. Further below we examine the kinetics of some possible processes that can lead to O2 formation.


Page 14 of 35

Page 15 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2. The energy diagram of the initial Ocus vacancy formation using DFT-PBE. The orange and dark blue spheres denote the O* atom and the new formed Ov, respectively. Values obtained from single-point HSE06 calculations are indicated in parenthesis.

We also investigated the barrier for vacancy formation on the pristine PdO(101) surface with the HSE06 functional. As for earlier calculations, we perform HSE06 single-point calculations but for the barrier we have evaluated the energy for the IS, TS, and FS and the results are also shown in Fig. 2. The resulting HSE06 barrier is 3.07 eV versus the PBE value of 1.28 eV. As shown in Fig. 2, HSE06 significantly destabilizes both the transition and final states with O* on the Pdcus row (by nearly 1.4 eV), which in turn leads to the very large HSE06 barrier for O vacancy formation of 3.07 eV. Because both the TS and FS are destabilized in a similar fashion the reverse barrier for O vacancy healing shows little change from PBE to HSE06. To ensure that this result is not an artifact of the single-point HSE06 calculations, we have performed full relaxation within the HSE06 functional for the initial and final state. PBE (fully relaxed HSE06) gives a reaction 15 ACS Paragon Plus Environment


energy (Erxn), defined as the energy of the final state minus the energy of the initial state, for oxygen vacancy formation of 1.28 (2.51) eV. The large Erxn from the fully relaxed HSE06 calculation confirms that the hybrid functional does indeed destabilize the O* configuration. The high barrier with HSE06 suggests that PdO(101) would begin to thermally reduce at much higher temperatures than the observed ~700 K in the TPD experiments in UHV36, which in turn implies that the HSE06 destabilization of O* on the Pdcus row is fundamentally incorrect. We have performed further testing to confirm that HSE06 destabilizes adsorbed O* on the Pdcus rows by examining adsorption of O* on pristine PdO(101) and in the presence of O vacancies. The details of the calculations can be found in the Supporting Information (SI). Table S2 shows that HSE06 destabilizes O* independent of the presence of O vacancies. It is important to note that HSE06 is relatively accurate for O vacancy formation energy (see Section 3.1), O2 and CO adsorption on the PdO(101) surface17,39, and also for bulk PdO oxidation energy39,54. The failure with HSE06 occurs only when the O atom is adsorbed on the PdO(101) surface. We have explored oxygen atom adsorption on several other transition metal oxide surfaces (RuO2(110), IrO2(110), TiO2(110)) and the results are reported in Figure S1 in the SI. Similar to PdO(101), large destabilization of O* is found using HSE06 on the other oxide surfaces. The destabilization is mainly attributed to a downward shift in the d band center with HSE06. To our knowledge, the large destabilization of adsorbed O on oxide surfaces within HSE06 has not been extensively reported. The majority of studies applying HSE to examine O vacancies focus on the oxygen vacancy formation energies and not on barriers to surface oxygen vacancy formation69–71. Santos and co-workers did study the adsorption energy of O on metal-doped LaMO3(001) surfaces using both PBE and HSE0672. They also find large destabilization for HSE versus PBE (on the order of 1-2 eV) for select transition metal dopants. More work examining O atom adsorption on oxides 16 ACS Paragon Plus Environment

Page 16 of 35

Page 17 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


and comparisons to experiment are required to better characterize the error in HSE for such systems. In the discussion that follows on pathways subsequent to initial Ov formation, we will present PBE results and then discuss the accuracy of single-point HSE06 results. However, our general conclusion is that for processes that involve O* in only the IS or FS (i.e. not in both) the PBE results provide more reasonable barriers. O* diffusion along the Pdcus row Facile diffusion of O* away from the Ov site would open the possibility of O* atoms diffusing along the Pdcus row and reacting with other O* atoms to form O2. Such processes would favor the formation of uncorrelated oxygen vacancies, which would be reflected in a more random Ov distribution on the surface. To test the importance of O* mobility, all the diffusion barriers for O* near the Ov site were calculated and the resulting energy diagram is shown in Figure 3. The numerical values for the diffusion barriers of O* using both PBE and single-point HSE06 are also reported in Table 2. For different diffusion paths, the predicted energy barriers are not significantly different with the two functionals, therefore we focus on the PBE results. We note that HSE06 still destabilizes the O atoms on the Pdcus rows; however, since O* diffusion involves only states along the Pdcus row there is a cancellation of error for the diffusion barriers. After initial oxygen vacancy formation, O* resides on a site adjacent to Ov (labeled site c in Fig. 3). The neighboring O vacancy stabilizes O* and decreases its diffusion barrier (0.47 eV) to the equivalent site d along the Pdcus row. The barrier is 1.03 eV for O* to escape from the initial vacancy position (cb or de), but the reverse step (bc or ed) has a smaller barrier of 0.71 eV. Note that after the O* atom diffuses one site away from the O vacancy the diffusion barrier is essentially the same as for O* on a pristine PdO(101) surface (see Table 2). This result indicates



that the effect of an O vacancy on O* diffusion is local (only to the nearest neighbor). From the results in Fig. 3, we would expect O* to be trapped adjacent to the vacancy site due to the kinetics of O* diffusion. However, the well depth of ~0.3 eV is not sufficient to prevent escape from these sites at the higher temperatures (~720 K) of the thermal reduction isothermal experiment.

Figure 3. Diffusion barrier diagram of O* using DFT-PBE. Different sites for O* are labeled (a-f) and the energies of the corresponding configurations are shown in the right panel.

Table 2. Calculated diffusion barriers (in eV) of O* using PBE and single-point HSE06 calculations. Primary results were obtained using a 4 × 1 unit cell but values using a 8 × 1 unit cell are shown in parenthesis. diffusion path a ↔ b, e ↔ f b → c, e → d c → b, d → e c↔d Perfect (See Figure 3) (equivalent) (toward Ov) (away from (equivalent) Surface Ov) PBE 0.65 (0.61) 0.71 (0.65) 1.03 (0.97) 0.47 (0.42) 0.64 (0.59) Single-point 0.58 0.59 1.11 0.69 0.58 HSE06


Page 18 of 35

Page 19 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


O2 formation pathways and desorption As shown above, after the initial oxygen vacancy formation, O* tends to go back to refill the vacancy with a small barrier (0.47 eV) or diffuse along the Pdcus row with a higher barrier. However, if O* can bind with another O atom to form O2 that then desorbs, O vacancies will begin to form on the surface. Therefore, the O2 formation and desorption pathways are critical to the evolution of O vacancies during thermal reduction. The adsorbed O* atom has two pathways to form O2 molecules: it can abstract an adjacent lattice Ocus atom (O* + Ocus  O2*) or react with another O* on the Pdcus row (O* + O*  O2*). In both of the two mechanisms, the barriers may also depend on the presence of neighboring Ocus vacancies. As a result, we have explored five distinct O2 formation pathways that incorporate adjacent oxygen vacancies and the schematic images are shown in Figure 4 along with the resulting PBE barriers. The single point HSE06 barriers associated with the same processes are shown in Table 3 in Section 3.3.

(a)

(b)

(c)

(d)

(e)

Figure 4. The schematic images of five possible O2 formation pathways investigated in DFT. The associated PBE barrier for the process is shown along with the reverse barrier in square brackets. Solid red circles correspond to O atoms (both O* and Ocus). Hollow black circles and solid blue squares indicate Ocus vacancies and Pdcus atoms, respectively. Two reaction mechanisms were studied: O* + Ocus  O2* with 0 (a) or 1 (b) neighboring Ov and O* + O*  O2* with 0 (c), 1 (d), or 2 (e) neighboring Ov. 19 ACS Paragon Plus Environment


Focusing on the first mechanism for O2 formation (O* + Ocus  O2*), we observe a decrease in the barrier to O2 formation from 2.50 to 1.87 eV due to the presence of an O vacancy. These values can be compared with the initial barrier for Ov formation of 1.71 eV. This comparison suggests that an O* atom is unlikely to abstract an Ocus atom in the absence of an adjacent Ov species being present. Furthermore, the sequence of initial Ov formation (barrier of 1.71 eV) followed by the new O* abstracting an adjacent Ocus atom (barrier of 1.87 eV) involves two steps with relatively large barriers. When added to the smaller reverse barriers for these two steps, these results would suggest that O2 formation through the first mechanism is difficult. For the second mechanism for O2 formation (O* + O*  O2*), the effect of adjacent Ov sites is even more significant. When the number of neighboring Ocus vacancies increases from 0 to 2, the association reaction barrier declines from 1.54 to 0.81 eV. Generally, the O* + O*  O2* step is more facile than O* + Ocus  O2* with a barrier lower by ~ 0.7 eV. Furthermore, the barrier for O2 dissociation on the Pdcus row is larger than O2 formation (by ~ 0.5 to 0.6 eV – see Fig. 4(c) and 4(d)) when there are one or less adjacent Ov. In the presence of 2 Ov (Fig. 4(e)) the forward and reverse barriers are relatively close. Both the lower association (i.e. formation) O2 barriers and relative O2 formation/dissociation barriers should favor O2 formation on the Pdcus rows over the process of O* abstracting Ocus to form O2. Based on the results in Fig. 4, we can envision a sequence where initial O* atoms are formed separately through the Ov formation pathway (Figure 2) and then forming O2* through diffusion along the Pdcus row (i.e. the second O2 formation pathway O* + O*  O2*). However, this scenario requires two separate Ov formation events, which with no neighbor vacancy effects would involve barriers of 1.71 eV each. In the next section, we will examine the effect of adjacent Ov sites on the Ov formation barriers.


Page 20 of 35

Page 21 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The single-point HSE06 results (see Table 3 in Sec 3.3) for the O2 pathways in Fig. 4 show a similar trend and therefore the Ov neighbor effects remain the same. As we have discussed earlier, HSE06 predicts that the transfer of an Ocus atom to the bridging Pdcus site has a larger barrier than found with PBE. As a result, HSE06 reduces (increases) the barriers of O* + O*  O2* (O* + Ocus  O2*), so the former association reaction becomes even more favorable with HSE06. However, we would expect that the PBE results are more accurate for the O2 formation pathways. Finally, we note that while O vacancies facilitate O2 formation, they also stabilize the final O2 state and the desorption barrier for O2 increases from 1.09 to 1.56 to 2.20 eV as we go from 0, 1, and 2 adjacent O vacancies. Therefore, while oxygen vacancies aid in the formation of surface O2 they also hinder the desorption of O2 from the surface. The O2 desorption energy is reduced by 0.35-0.45 eV when using HSE06 (see Table 3), but the trend of increasing O2 desorption barrier with adjacent Ov remains. We will discuss the impact of HSE06 O2 desorption energies to the overall picture of thermal reduction of the PdO(101) in more detail further below. Kinetics of Ov chains formation As shown above, the presence of Ov sites significantly affects subsequent processes that involve the neighboring atoms. Experimental STM images also show the formation of Ov chains during thermal reduction of PdO(101) suggesting that oxygen vacancies promote adjacent oxygen vacancy formation. Hence, we evaluated the barriers for further Ov formation in the presence of existing oxygen vacancies under various scenarios shown in Figure 5. Comparing Fig. 5(a)-(c), we find a decrease in the oxygen vacancy barrier due to the presence of adjacent oxygen vacancies. The barrier declines by nearly 45% when adding one adjacent vacancy (1.71 eV to 0.96 eV from Fig. 5(a) to 5(b)) and drops an additional ~ 20% with a second adjacent vacancy (Fig. 5(c)). This decrease in barrier is associated with the added stability for the TS for O* formation as it moves 21 ACS Paragon Plus Environment


from the lattice site to the Pdcus row due to the removal of repulsive interactions with the neighboring Ocus atom. To test this interpretation, we have examined two other scenarios. Figure 5(d) shows an O vacancy pathway moving away from the existing O vacancy and the resulting barrier is very similar to what is observed on a perfect surface (Fig. 5(a)). Similarly, when an O* atom remains next to an oxygen vacancy we find the barrier to O vacancy formation of the adjacent Ocus atom to also have a high barrier as on the perfect surface (see Fig. 5(e)). We note that in the presence of the O* atom the reverse barrier to Ov healing is larger than what is found on the defectfree PdO(101). This increase in the reverse barrier is due to the extra stability of a pair of O* atoms adjacent to two Ov sites. In fact, from the final state in Fig. 5(e) the pair of O* atoms can form an O2 molecule through the pathway shown in Fig. 4(e) with a barrier of 0.81 eV. Based on the results in Fig. 5, we conclude that existing Ov sites can facilitate neighboring Ov formation (1.71 eV0.96 eV for 1 adjacent Ov and 0.96 eV  0.74 eV with 2 adjacent Ov). Our expectation is that this effect will taper off after 2 Ov sites based on our findings for 1 and 2 Ov (i.e. we do not expect a further decrease in barrier for Ov formation with 3 Ov) but we have not explicitly tested this hypothesis. The presence of O* next to an Ov site can hinder this reduction in barriers (see Fig. 5(e)). However, when comparing the barriers for O* diffusion away from an Ov site (~ 1 eV) with that for forming Ov sites (~1.7 eV) we expect that scenarios where the O* atom diffuses away from the Ov site will occur frequently enough to provide opportunities for the growth of Ov chains.


Page 22 of 35

Page 23 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(a)

(b)

(c)

(d)

(e)

Figure 5. Schematic images of oxygen vacancy pathways in the presence of oxygen vacancies. The barrier to the process is indicated with the barrier for the reverse barrier in square brackets. The key follows Fig. 4 except the solid orange circle represents the final position of the Ocus atom as it moves from the initial location onto the Pdcus row (the arrow is provided to aid in understanding this diffusion). The oxygen vacancy formation pathway with (a) zero (b) one and (c) two adjacent O vacancies, along with (d) the oxygen vacancy formation pathway moving away from an existing O vacancy and (e) in the presence of both an O* and an adjacent oxygen vacancy.

We also investigated the pathways in Fig. 5 using single-point HSE06 method (see Table 3 in Section 3.3 for the HSE06 values). The results show that the overall trend for oxygen vacancy formation barriers in the presence of oxygen vacancies are consistent with the PBE results (i.e. adjacent O vacancies lower the barriers); however, the barriers are ~ 0.6 eV higher than the results obtained from DFT-PBE. As mentioned earlier, we expect that the PBE results are more accurate based on the experimental temperature range of the thermal reduction of the PdO(101) surface. Regeneration of Ocus atoms by subsurface O atoms In an earlier study, we identified a pathway for subsurface O to diffuse into a surface Ocus vacancy that has a barrier for 1.61 eV17. We have also examined this mechanism in the presence of an Ov dimer and we find that the barrier is nearly the same (1.67 eV with 2 Ov). The energy 23 ACS Paragon Plus Environment


barrier for Ov healing from the subsurface is ~ 1.6 eV, which is similar to the formation energy of a single vacancy (1.71 eV) but substantially larger than that of vacancy dimer formation (0.96 eV) and subsequent Ov chain growth (0.74 eV). This finding also supports the clustering of vacancies because an isolated single vacancy can be eliminated by subsurface O atom diffusion. We also investigated if this transfer could happen from deeper within the PdO(101) surface by examining O atom diffusion in a 6 layer thick slab and in bulk PdO. We used the same 4 × 1 supercell for the 6 layer thick PdO film and a 2 × 2 × 2 supercell model for bulk PdO. The O diffusion barrier from a third layer to a second layer is found to be 1.90 eV which is ~ 0.3 eV higher than the barrier from second layer to first layer. This barrier dramatically increases to 2.60 eV when O atom diffuses from a fourth layer to a third layer. This value is very similar to our calculated barrier of O diffusion in bulk PdO (~ 2.8 eV). The results imply that there is a near surface region (second and third layers) where subsurface O atoms can more readily diffuse into the surface and that below this near surface region subsurface O atoms are almost immobile at temperatures near 700 K. However, we also expect that other processes involving additional rearrangement of surface atoms could occur as the surface reduction proceeds. Such complex processes that occur later in the reduction processes are difficult to predict a priori and beyond the scope of this paper. 3.3 Discussion of Ov neighbor effect on thermal reduction of PdO(101) Here, we connect the elementary steps that we have studied in the previous sections to observations from the experimental isothermal reduction of the PdO(101) film in UHV36. The PBE barriers of the elementary processes with or without neighboring Ov are summarized in Table 3. Single-point HSE06 derived barriers for all the processes are also shown in the table. We note that the negative values for some HSE06 barriers in Table 3 (i.e. O*  Ocus and O* + O*  O2* in the presence of 2 Ov) are an artifact of the single point HSE06 approach to estimate barriers. In the 24 ACS Paragon Plus Environment

Page 24 of 35

Page 25 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


presence of 2 Ov sites the structure undergoes more changes and the use of single-point calculations will be less accurate. Table 3. The summary of all the processes involved in the evolution of surface Ov and O2 desorption. The single-point HSE06 results are put into parenthesis. Process

Eb (No Ov) / eV

Eb (1 Ov) / eV

Eb (2 Ov) / eV

Ocus  O*

1.71 (3.07)

0.96 (1.58)

0.74 (1.25)

O*  Ocus

-

0.47 (0.43)

0.33 (-0.25)

O* + Ocus  O2*

2.50 (2.61)

1.87 (2.31)

-

O2*  O* +Ocus

-

1.60 (1.71)

1.22 (1.58)

O* + O*  O2*

1.54 (1.51)

1.20 (1.08)

0.81 (-0.13)

O2*  O* + O*

2.11 (3.03)

1.70 (2.60)

0.70 (0.59)

O2*  O2 (gas)

1.09 (0.64)

1.56 (1.12)

2.20 (1.85)

O*  O* (diffusion)

0.64 (0.58)

toward Ov: 0.71 (0.59)

-

away from Ov: 1.03 (1.11)

-

Osub  Ocus

-

1.61 (2.02)

1.67 (2.20)

Ocus  Osub

2.18 (2.52)

2.34 (2.83)

-

To understand the impact of the barriers reported in Table 3 to the early stages of reduction of the PdO(101) surface, it is helpful to focus on the pathways to O2 formation and desorption, since this is the mechanism that leads to irreversible thermal reduction in the experiment. As discussed earlier, there are two main pathways to O2 formation: (1) O* from the Pdcus row directly reacts with Ocus atoms and (2) the association reaction occurs on the Pdcus row. In general, the barriers are substantially lower for the second pathway and our expectation is the O2 formation occurs primarily through O* association on the Pdcus row. If we ignore the effect of neighbor oxygen vacancies, the difference between the two O2 formation pathways is nearly 1 eV. Based on that large energy difference, we would expect that the primary pathway would involve O vacancy 25 ACS Paragon Plus Environment


formation and then O* diffusion and association reaction on the Pdcus row before O2 desorption. Among those steps the largest barrier is associated with Ov formation (1.71 eV), which as we have noted earlier would be a feasible value in comparison to the onset of thermal reduction in the TPD experiments of ~ 700 K. Comparing this barrier to O* diffusion (~ 0.6 eV), we can assume that O* diffusion is much faster than O* formation. Therefore, under the scenario where we ignore any Ov neighbor effects, O2 formation and desorption will not be spatially correlated to the location of Ov sites (i.e. where O* is initially formed). This process will result in oxygen vacancies randomly distributed on the PdO(101) surface, which is clearly not observed in the STM images36. We now consider thermal reduction incorporating the neighbor effect due to O vacancies. When an oxygen vacancy forms, the resulting O* has a barrier to diffusing away from the initial location (about a 0.3 eV bias to stay adjacent to the Ov), but the impact of this barrier will be relatively small in comparison to the high temperatures required to initiate Ov formation. While O* trapping by adjacent Ov (or Ov chains) may occur, we would expect it to be minor under thermal reduction conditions (this conclusion may be different for lower temperature conditions for CO oxidation). The much larger impact of the oxygen vacancies is on the Ov formation barrier. After a few vacancies form on the surface, the existing vacancies activate neighboring Ocus atoms to form Ov with much lower barriers (0.96 eV versus 1.71 eV ), which generates Ov dimers. The barrier to forming additional Ov sites next to Ov dimers is much more favorable than isolated Ov formation and subsequently long chains of Ov should start to form as observed in STM images. Note that the much smaller reverse barrier to Ov healing is only weakly dependent on the presence of adjacent oxygen vacancies (0.47 to 0.33 eV for 1 Ov versus 2 Ov). The small barrier to Ov healing will also enhance Ov chain formation over isolated Ov sites, since any isolated Ov site will likely become healed. 26 ACS Paragon Plus Environment

Page 26 of 35

Page 27 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


As Ov chains grow on the PdO(101) surface, the probability of O2 formation events occurring in their presence will increase. Such O2 formation processes have lower barriers than when they occur on the defect-free PdO(101) surface, which should lead to an acceleration of the reduction of the PdO(101) surface. This acceleration is observed in the isothermal O2 desorption traces experimentally, where there is an initial rise of the O2 desorption rate as expected but instead of leveling off quickly there is a dramatic increase in the rate of O2 desorption, which is referred to as the autocatalytic regime36. This observation is generally in agreement with our DFT findings. As the reduction proceeds more O* are favorably generated on the surface and the kinetics of O2 formation becomes more facile thereby systematically promoting O2 desorption. As a result, the O2 desorption rate increases as the surface reduces, which is the signature of the autocatalytic behavior. As we have noted throughout the paper, the HSE06 results generally follow the trend from PBE with one significant exception where the barrier of Ov formation (OcusO*) becomes substantially larger with HSE06 (by ~ 1.3 eV). This process is a key step in initiating PdO(101) thermal reduction and such a large barrier is not reconcilable with the experimental decomposition temperature of the PdO(101) surface. However, we have also established in an earlier DFT and TPRS study of O2 desorption on PdO(101) that the HSE06 values for O2 desorption are closer to experimental values39,73. As shown in Table 3, the PBE O2 desorption energy rises from 1.09 eV on the defect-free PdO(101) to 2.20 eV in the presence of a Ov dimer. In fact, comparing all of the values in Table 3 suggests that O2 desorption may become the rate limiting step (it will at least have the largest barrier in the presence of oxygen vacancies) as the surface becomes populated by long Ov chains. The difference between PBE and single-point HSE06 O2 desorption energy drops from 0.45 eV for the defect-free PdO(101) to 0.35 eV for a Ov dimer, but this value is still 27 ACS Paragon Plus Environment


sufficiently large to impact surface reduction kinetics. Our expectation is that HSE06 values will be more accurate for O2 desorption energies but future work comparing DFT-based kinetic Monte Carlo (kMC) simulations directly to isothermal desorption spectra will be needed to confirm this expectation. Overall, the DFT calculations clearly indicate that the presence of oxygen vacancies promotes the kinetics of further adjacent Ov formation, which leads to the observed long Ov chains and likely triggers the experimentally observed autocatalytic behavior. However, because the formation of Ov chains is linked to the formation and desorption of O2 which can occur in a myriad of ways, kMC simulations that integrate all the surface processes are needed to confirm the dominant processes that mediate autocatalytic behavior. Such kMC simulations could also be compared to experimental desorption curves (over the regime of early state thermal reduction) to resolve if HSE06 O2 desorption energies are more accurate for the PdO(101) surface. 4. Summary We studied the initial reduction of the PdO(101) surface using DFT calculations with the goal of understanding UHV isothermal reduction experiments that show the formation of Ov chains and corresponding acceleration of the rate of surface reduction. Both PBE and single-point HSE06 calculations were performed based on an earlier study which indicates that PBE overestimates the O2 binding energy while HSE06 values are close to binding energies derived from O2 TPD experiments39. An evaluation of the formation energy of Ocus and O4f surface atom vacancies indicates that Ocus vacancies are more favorable. For the Ocus vacancies, both the in-row and acrossrow interactions are negligible, which suggests that any clustering of O vacancies is not driven by a thermodynamic preference. We systematically evaluated the barriers for Ocus vacancy formation, 28 ACS Paragon Plus Environment

Page 28 of 35

Page 29 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


diffusion of the resulting O*, O2* formation through a myriad of pathways, and O2* desorption. All of these processes were evaluated on the defect-free PdO(101) surface and in the presence of Ocus vacancies. Our DFT calculations reveal that neighboring Ocus vacancies have a strong effect on the barriers of processes occurring on adjacent sites. Vacancy formation (Ocus  O*) on a defect-free PdO(101) surface has a PBE barrier of 1.71 eV with a reverse barrier (vacancy healing) of 0.47 eV. These barriers decrease to 0.74 and 0.33 eV, respectively, in the presence of two adjacent oxygen vacancies. The DFT results indicate that there is a kinetic preference for Ov dimers and chains to form on the PdO(101) surface, as observed in STM experiments. While PBE and HSE06 provide qualitatively similar results for the dependence of barriers on oxygen vacancies, HSE06 predicts unrealistically large barriers for O vacancy formation. Comparison between PBE and HSE06 values suggests that HSE06 is not accurate for configurations containing O*. Furthermore, we have identified an O2 formation pathway on the Pdcus row (O* + O*  O2*) that also becomes more facile in the presence of oxygen vacancies. The barrier for O2 formation in this mechanism drops from 1.54 eV on a defect-free surface to 0.81 eV with a pair of adjacent Ov sites. Therefore, the preference for Ov chain formation also leads to more facile O2 formation, which should lead to an acceleration of the reduction of the surface as observed in isothermal desorption measurements in UHV at 720 K. In addition, to the main conclusion above we propose that HSE06 barriers for O2 desorption are more accurate than PBE for both defect-free and reduced PdO(101) surfaces. Future kMC simulations and additional surface science experiments will be required to confirm the differences between HSE06 and PBE functionals for these two important processes (i.e. Ov formation and O2 desorption). Overall, our DFT study shows the importance of incorporating neighbor effects of oxygen vacancies to accurately capture the reduction of the 29 ACS Paragon Plus Environment


PdO(101) surface. The findings in this paper may also have important consequences in understanding oxidation catalysis on this surface, where the Mars-van Krevelen process can be affected by the neighbor effects of oxygen vacancy chains. Supporting Information Testing of spin polarization for select processes on the PdO(101) surface. Energy of O* using single point HSE06 functional in comparison to PBE functional on PdO(101) and other TM oxides.

Acknowledgements We acknowledge the Ohio Supercomputing Center for providing computational resources. We gratefully acknowledge financial support for this work provided by the Department of Energy, Office of Basic Energy Sciences, Catalysis Science Division through Grant DE-FG0203ER15478.

References (1) (2) (3) (4) (5) (6) (7)

Spronsen, M. A. Van; Frenken, J. W. M.; Groot, I. M. N. Surface Science under Reaction Conditions : CO Oxidation on Pt and Pd Model Catalysts. Chem. Soc. Rev. 2017, 46, 4347–4374. Weaver, J. F.; Zhang, F.; Pan, L.; Li, T.; Asthagiri, A. Vacancy-Mediated Processes in the Oxidation of CO on PdO(101). Acc. Chem. Res. 2015, 48, 1515–1523. Weaver, J. F.; Hakanoglu, C.; Antony, A.; Asthagiri, A. Alkane Activation on Crystalline Metal Oxide Surfaces. Chem. Soc. Rev. 2014, 43, 7536. Matera, S.; Blomberg, S.; Hoffmann, M. J.; Zetterberg, J.; Gustafson, J.; Lundgren, E.; Reuter, K. Evidence for the Active Phase of Heterogeneous Catalysts through In Situ Reaction Product Imaging and Multiscale Modeling. ACS Catal. 2015, 5, 4514–4518. Chen, M. S.; Cai, Y.; Yan, Z.; Gath, K. K.; Axnanda, S.; Goodman, D. W. Highly Active Surfaces for CO Oxidation on Rh, Pd, and Pt. Surf. Sci. 2007, 601, 5326–5331. McClure, S. M.; Goodman, D. W. New Insights into Catalytic CO Oxidation on Pt-Group Metals at Elevated Pressures. Chem. Phys. Lett. 2009, 469, 1–13. Gao, F.; McClure, S. M.; Cai, Y.; Gath, K. K.; Wang, Y.; Chen, M. S.; Guo, Q. L.; Goodman, D. W. CO Oxidation Trends on Pt-Group Metals from Ultrahigh Vacuum to near Atmospheric Pressures: A Combined in Situ PM-IRAS and Reaction Kinetics Study. 30 ACS Paragon Plus Environment

Page 30 of 35

Page 31 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(8) (9) (10)

(11) (12) (13)

(14) (15) (16) (17) (18) (19) (20) (21) (22)

Surf. Sci. 2009, 603, 65–70. Chen, M.; Wang, X. V.; Zhang, L.; Tang, Z.; Wan, H. Active Surfaces for CO Oxidation on Palladium in the Hyperactive State. Langmuir 2010, 26, 18113–18118. Gao, F.; Wang, Y.; Cai, Y.; Goodman, D. W. CO Oxidation on Pt-Group Metals from Ultrahigh Vaccum to Near Atmospheric Pressures. 2.Palladium and Platinum. J. Phys. Chem. C 2009, 113, 174–181. Blomberg, S.; Hoffmann, M. J.; Gustafson, J.; Martin, N. M.; Fernandes, V. R.; Borg, A.; Liu, Z.; Chang, R.; Matera, S.; Reuter, K.; et al. In Situ X-Ray Photoelectron Spectroscopy of Model Catalysts: At the Edge of the Gap. Phys. Rev. Lett. 2013, 110, 117601. Onderwaater, W. G.; Taranovskyy, A.; Baarle, G. C. van; Frenken, J. W. M.; Groot, I. M. N. In Situ Optical Reflectance Difference Observations of CO Oxidation over Pd(100). J. Phys. Chem. C 2017, 121, 11407–11415. van Rijn, R.; Balmes, O.; Resta, A.; Wermeille, D.; Westerström, R.; Gustafson, J.; Felici, R.; Lundgren, E.; Frenken, J. W. M. Surface Structure and Reactivity of Pd(100) during CO Oxidation near Ambient Pressures. Phys. Chem. Chem. Phys. 2011, 13, 13167. Van Rijn, R.; Balmes, O.; Felici, R.; Gustafson, J.; Wermeille, D.; Westerström, R.; Lundgren, E.; Frenken, J. W. M. Comment on “CO Oxidation on Pt-Group Metals from Ultrahigh Vacuum to near Atmospheric Pressures. 2. Palladium and Platinum.” J. Phys. Chem. C 2010, 114, 6875–6876. Toyoshima, R.; Yoshida, M.; Monya, Y.; Kousa, Y.; Suzuki, K.; Abe, H.; Mun, B. S.; Mase, K.; Amemiya, K.; Kondoh, H. In Situ Ambient Pressure XPS Study of CO Oxidation Reaction on Pd (111) Surfaces. J. Phys. Chem. C 2012, 116, 18691–18697. Gustafson, J.; Olivier, B.; Zhang, C.; Shipilin, M.; Merte, L. R.; Schaefer, A.; Hagman, B.; Merte, L. R.; Marin, N. M.; Carlsson, P.; et al. The Role of Oxides in Catalytic CO Oxiation over Rhodium and Palladium. ACS Catal. 2018, 8, 4438–4445. Zhang, F.; Li, T.; Pan, L.; Asthagiri, A.; Weaver, J. F. Catalysis Science & Technology CO Oxidation on Single and Multilayer Pd Oxides on Pd(111): Mechanistic Insights from RAIRS. Catal. Sci. Technol. 2014, 4, 3826–3834. Zhang, F.; Pan, L.; Li, T.; Diulus, J. T.; Asthagiri, A.; Weaver, J. F. CO Oxidation on PdO(101) during Temperature-Programmed Reaction Spectroscopy: Role of Oxygen Vacancies. J. Phys. Chem. C 2014, 118, 28647–28661. Weaver, J. F.; Choi, J.; Mehar, V.; Wu, C. Kinetic Coupling among Metal and Oxide Phases during CO Oxidation on Partially Reduced PdO(101): Influence of Gas-Phase Composition. ACS Catal. 2017, 7, 7319–7331. Zhou, J.; Blomberg, S.; Gustafson, J.; Lundgren, E.; Zetterberg, J. Simultaneous Imaging of Gas Phase over and Surface Reflectance of a Pd(100) Single Crystal during CO Oxidation. J. Phys. Chem. C 2017, 121, 23511–23519. Zhou, J.; Blomberg, S.; Gustafson, J.; Lundgren, E.; Zetterberg, J. Visualization of Gas Distribution in a Model AP-XPS Reactor by PLIF: CO Oxidation over a Pd(100) Catalyst. Catalysts 2017, 7, 29. Blomberg, S.; Zetterberg, J.; Gustafson, J.; Zhou, J.; Brackmann, C.; Lundgren, E. Comparison of AP-XPS and PLIF Measurements during CO Oxidation over Pd Single Crystals. Top. Catal. 2016, 59, 478–486. Liu, D. J.; Zahariev, F.; Gordon, M. S.; Evans, J. W. Predictive Beyond-Mean-Field Rate Equations for Multisite Lattice− Gas Models of Catalytic Surface Reactions: CO 31 ACS Paragon Plus Environment


(23) (24) (25) (26) (27) (28) (29) (30) (31) (32) (33) (34) (35) (36) (37) (38) (39) (40)

Oxidation on Pd(100). J. Phys. Chem. C 2016, 120, 28639–28653. Cisternas, J.; Karpitschka, S.; Wehner, S. Travelling Fronts of the CO Oxidation on Pd(111) with Coverage-Dependent Diffusion. J. Chem. Phys. 2014, 141, 164106. Liu, D. J.; Evans, J. W. Atomistic Lattice-Gas Modeling of CO Oxidation on Pd(100): Temperature-Programed Spectroscopy and Steady-State Behavior. J. Chem. Phys. 2006, 124, 154705. Piccinin, S.; Stamatakis, M. Steady-State CO Oxidation on Pd(111): First-Principles Kinetic Monte Carlo Simulations and Microkinetic Analysis. Top. Catal. 2017, 60, 141– 151. Piccinin, S.; Stamatakis, M. CO Oxidation on Pd(111): A First-Principles-Based Kinetic Monte Carlo Study. ACS Catal. 2014, 4, 2143–2152. Reuter, K.; Scheffler, M. First-Principles Kinetic Monte Carlo Simulations for Heterogeneous Catalysis, Applied to the CO Oxidation at RuO2(110). Phys. Rev. B 2006, 73, 045433. Hess, F.; Farkas, A.; Seitsonen, A. P.; Over, H. “First-Principles” Kinetic Monte Carlo Simulations Revisited: CO Oxidation over RuO2(110). J. Comput. Chem. 2012, 33, 757– 766. Sutton, J. E.; Lorenzi, J. M.; Krogel, J. T.; Xiong, Q.; Pannala, S.; Matera, S.; Savara, A. Electrons to Reactors Multiscale Modeling: Catalytic CO Oxidation over RuO2. ACS Catal. 2018, 8, 5002–5016. Hoffmann, M. J.; Scheffler, M.; Reuter, K. Multi-Lattice Kinetic Monte Carlo Simulations from First Principles: Reduction of the Pd(100) Surface Oxide by CO. ACS Catal. 2015, 5, 1199–1209. Rogal, J.; Reuter, K.; Scheffler, M. First-Principles Statistical Mechanics Study of the Stability of a Subnanometer Thin Surface Oxide in Reactive Environments: CO Oxidation at Pd(100). Phys. Rev. Lett. 2007, 98, 046101. Jørgensen, M.; Grönbeck, H. First-Principles Microkinetic Modeling of Methane Oxidation over Pd(100) and Pd(111). ACS Catal. 2016, 6, 6730–6738. Van Den Bossche, M.; Grönbeck, H. Methane Oxidation over PdO(101) Revealed by First-Principles Kinetic Modeling. J. Am. Chem. Soc. 2015, 137, 12035–12044. Duan, Z.; Henkelman, G. CO Oxidation on the Pd(111) Surface. ACS Catal. 2014, 4, 3435–3443. Chin, Y. C.; Garcia-Dieguez, M.; Iglesia, E. Dynamics and Thermodynamics of Pd − PdO Phase Transitions : Effects of Pd Cluster Size and Kinetic Implications for Catalytic Methane Combustion. J. Phys. Chem. C 2016, 120, 1446–1460. Hinojosa, J. A.; Weaver, J. F. Surface Structural Evolution during the Thermal Decomposition of a PdO(101) Thin Film. Surf. Sci. 2011, 605, 1797–1806. Heyd, J.; Scuseria, G. E.; Ernzerhof, M. Hybrid Functionals Based on a Screened Coulomb Potential. J. Chem. Phys. 2003, 118, 8207–8215. Heyd, J.; Scuseria, G. E. Efficient Hybrid Density Functional Calculations in Solids: Assessment of the Heyd-Scuseria-Ernzerhof Screened Coulomb Hybrid Functional. J. Chem. Phys. 2004, 121, 1187–1192. Pan, L.; Weaver, J. F.; Asthagiri, A. First Principles Study of Molecular O2 Adsorption on the PdO(101) Surface. Top. Catal. 2017, 60, 401–412. Doornkamp, C.; Ponec, V. The Universal Character of the Mars and Van Krevelen Mechanism. J. Mol. Catal. A Chem. 2000, 162, 19–32. 32 ACS Paragon Plus Environment

Page 32 of 35

Page 33 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(41) (42) (43) (44) (45) (46) (47) (48) (49) (50) (51) (52) (53) (54) (55) (56) (57) (58) (59) (60) (61)

Weaver, J. F. Surface Chemistry of Late Transition Metal Oxides. Chem. Rev. 2013, 113, 4164–4215. Kresse, G.; Hafner, J. Ab Initio Molecular Dynamics for Liquid Metals. Phys. Rev. B 1993, 47, 558(R). Kresse, G. Ab Initio Molecular Dynamics for Liquid Metals. J. Non. Cryst. Solids 1995, 192&193, 222–229. Blöchl, P. E. Projector Augmented-Wave Method. Phys. Rev. B 1994, 50, 17953–17979. Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865–3868. Monkhorst, H. J.; Pack, J. D. Special Points for Brillouin-Zone Integrations. Phys. Rev. B 1976, 13, 5188–5192. Weaver, J. F.; Hankanoglu, C.; Antony, A.; Astharigi, A. High Selectivity for Primary CH Bond Cleavage of Propane σ-Complexes on the PdO(101) Surface. J. Am. Chem. Soc. 2011, 133, 16196–16200. Hakanoglu, C.; Hawkins, J. M.; Asthagiri, A.; Weaver, J. F. Strong Kinetic Isotope Effect in the Dissociative Chemisorption of H2 on a PdO (101) Thin Film. J. Phys. Chem. C 2010, 114, 11485–11497. Antony, A.; Asthagiri, A.; Weaver, J. F. Pathways for C–H Bond Cleavage of Propane σComplexes on PdO(101). Phys. Chem. Chem. Phys. 2012, 14, 12202–12212. Kan, H. H.; Weaver, J. F. A PdO(101) Thin Film Grown on Pd(111) in Ultrahigh Vacuum. Surf. Sci. 2008, 602, L53–L57. Kan, H. H.; Weaver, J. F. Mechanism of PdO Thin Film Formation during the Oxidation of Pd(111). Surf. Sci. 2009, 603, 2671–2682. Henkelman, G.; Uberuaga, B. P.; Jónsson, H. A Climbing Image Nudged Elastic Band Method for Finding Saddle Points and Minimum Energy Paths. J. Chem. Phys. 2000, 113, 9901–9904. Hu, Z.; Metiu, H. Choice of U for DFT+U Calculations for Titanium Oxides. J. Phys. Chem. C 2011, 115, 5841–5845. Bruska, M. K.; Czekaj, I.; Delley, B.; Mantzaras, J.; Wokaun, A. Electronic Structure and Oxygen Vacancies in PdO and ZnO: Validation of DFT Models. Phys. Chem. Chem. Phys. 2011, 13, 15947–15954. Kurzman, J. A.; Miao, M. S.; Seshadri, R. Hybrid Functional Electronic Structure of PbPdO2, a Small-Gap Semiconductor. J. Phys. Condens. Matter 2011, 23, 465501. Gopal, P.; De Gennaro, R.; Gusmao, M. S.; Orabi, R. A. R. Al; Wang, H.; Curtarolo, S.; Fornari, M.; Nardelli, M. B. Improved Electronic Structure and Magnetic Exchange Interactions in Transition Metal Oxides. J. Phys. Condens. Matter 2017, 29, 444003. Okamoto, H.; Asô, T. Formation of Thin Films of PdO and Their Electric Properties. Jpn. J. Appl. Phys. 1967, 6, 779. Rey, E.; Kamal, M. R.; Miles, R. B.; Royce, B. S. H. The Semiconductivity and Stability of Palladium Oxide. J. Mater. Sci. 1978, 13, 812–816. Nilsson, P. O.; Shivaraman, M. S. Optical Properties of PdO in the Range of 0.5-5.4 EV. J. Phys. C Solid State Phys. 1979, 12, 1423–1427. Pillo, T.; Zimmermann, R.; Steiner, P.; Hüfner, S. The Electronic Structure of PdO Found by Photoemission (UPS and XPS) and Inverse Photoemission (BIS). J. Phys. Condens. Matter 1997, 9, 3987–3999. Park, K.; Novikov, D. L.; Gubanov, V. A.; Freeman, A. J. Electronic Structure of Noble33 ACS Paragon Plus Environment


(62) (63)

(64) (65) (66) (67) (68) (69) (70) (71) (72) (73)

Metal Monoxides: PdO, PtO, and AgO. Phys. Rev. B 1994, 49, 4425–4431. Blanco-Rey, M.; Jenkins, S. J. Methane Dissociation and Methyl Diffusion on PdO(100). J. Chem. Phys. 2009, 130, 014705. Mehar, V.; Kim, M.; Shipilin, M.; Van Den Bossche, M.; Gustafson, J.; Merte, L. R.; Hejral, U.; Gronbeck, H.; Lundgren, E.; Asthagiri, A.; et al. Understanding the Intrinsic Surface Reactivity of Single and Multilayer PdO(101) on Pd(100). ACS Catal. 2018, 8, 8553–8567. Van Den Bossche, M.; Martin, N. M.; Gustafson, J.; Hakanoglu, C.; Weaver, J. F.; Lundgren, E.; Grönbeck, H. Effects of Non-Local Exchange on Core Level Shifts for GasPhase and Adsorbed Molecules. J. Chem. Phys. 2014, 141, 034706. Zhang, X.; Ptasinska, S. Dissociative Adsorption of Water on an H2O / GaAs(100) Interface : In Situ Near-Ambient Pressure XPS Studies. J. Phys. Chem. C 2014, 118, 4259–4266. Hakanoglu, C.; Hawkins, J. M.; Asthagiri, A.; Weaver, J. F. Strong Kinetic Isotope Effect in the Dissociative Chemisorption of H2 on a PdO( 101 ) Thin Film. J. Phys. Chem. C 2010, 114, 11485–11497. Hinojosa, J. A.; Antony, A.; Hakanoglu, C.; Asthagiri, A.; Weaver, J. F. Adsorption of CO2 on a PdO(101) Thin Film. J. Phys. Chem. C 2012, 116, 3007–3016. Antony, A.; Hakanoglu, C.; Asthagiri, A.; Weaver, J. F. Dispersion-Corrected Density Functional Theory Calculations of the Molecular Binding of n-Alkanes on Pd(111) and PdO(101). J. Chem. Phys. 2012, 136, 054702. Janotti, A.; Varley, J. B.; Rinke, P.; Umezawa, N.; Kresse, G.; Van de Walle, C. G. Hybrid Functional Studies of the Oxygen Vacancy in TiO2. Phys. Rev. B 2010, 81, 085212. Krcha, M. D.; Janik, M. J. Examination of Oxygen Vacancy Formation in Mn-Doped CeO2(111) Using DFT+U and the Hybrid Functional HSE06. Langmuir 2013, 29, 10120– 10131. Sun, L.; Huang, X.; Wang, L.; Janotti, A. Disentangling the Role of Small Polarons and Oxygen Vacancies in CeO2. Phys. Rev. B 2017, 95, 245101. Santos, E. J. G.; Nørskov, J. K.; Vojvodic, A. Screened Hybrid Exact Exchange Correction Scheme for Adsorption Energies on Perovskite Oxides. J. Phys. Chem. C 2015, 119, 17662–17666. Hinojosa, J. A.; Kan, H. H.; Weaver, J. F. Molecular Chemisorption of O2 on a PdO(101) Thin Film on Pd(111). J. Phys. Chem. C 2008, 112, 8324–8331.


Page 34 of 35

Page 35 of 35

TOC graphic

Ov formation on perfect PdO(101) Ov formation on reduced PdO(101)

2

44 % 0.96

1

Reduced PdO(101) 0.47

0.5 0

Perfect PdO(101)

1.71

1.5

Energy (eV)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0.0

0.33


Initial Reduction of the PdO(101) Surface: Role of Oxygen Vacancy

Recommend Documents