Predictive Beyond-Mean-Field Rate Equations for ... - ACS Publications

Article pubs.acs.org/JPCC

Predictive Beyond-Mean-Field Rate Equations for Multisite Lattice− Gas Models of Catalytic Surface Reactions: CO Oxidation on Pd(100) Published as part of The Journal of Physical Chemistry virtual special issue “Mark S. Gordon Festschrift”. Da-Jiang Liu,*,† Federico Zahariev,† Mark S. Gordon,†,‡ and James W. Evans*,†,§,∥ †

Ames Laboratory − USDOE, ‡Department of Chemistry, §Department of Mathematics, and ∥Department of Physics & Astronomy, Iowa State University, Ames, Iowa 50011, United States S Supporting Information *

ABSTRACT: Tailored multisite lattice−gas (msLG) models are developed for CO oxidation on Pd(100) at low-pressures. These models include multiple adsorption site types and superlattice adlayer ordering due to short-range exclusion for highly mobile reactant adspecies. However, they are simplified to neglect longer-range weaker adspecies interactions, so that the key energetic parameters are the CO desorption barrier and the reaction barrier. We discuss existing density functional theory results for these energies and present additional analysis for CO adsorption. After also including an appropriate nontrivial specification of the dynamics of adsorption onto mixed reactant adlayers, we develop rate equations for the reaction kinetics. Our formulation goes beyond traditional mean-field (MF) Langmuirian treatments by accounting for multiple adsorption sites and for the strong spatial correlations associated with superlattice ordering. Specifically, we utilize factorization approximations based on appropriate site motifs, and also Padé resummation of exact low-coverage expansions for sticking coefficients. Our beyond-MF rate equations are successful in accurately predicting key aspects of reactive steady-state behavior, and thus expand the utility of rate equation formulations in surface chemistry. This is confirmed by comparison with precise kinetic Monte Carlo simulation results. Specifically, we not only assess bistability and criticality observed for CO oxidation but also find more complex multistability associated with symmetry-breaking transitions in high-coverage CO adlayers.

1. INTRODUCTION Many recent theoretical studies of heterogeneous catalytic reactions on crystalline metal surfaces have focused on determination of key energetics and barriers exploiting density functional theory (DFT) calculations.1 Issues remain regarding the accuracy of such energetics. However, beyond reliable determination of energetics, an ultimate goal of theoretical studies is predictive description of the overall catalytic reaction kinetics and turnover frequency (TOF) often for flow reactors as a function of experimental control parameters (reactant partial pressures and surface temperature, T).2,3 This requires molecular-level modeling of all the relevant adsorption, diffusion, desorption, and reaction processes, accounting for the diverse local environments in which they can occur. Also, such modeling must have the capability to track behavior on the appropriate time- and length-scales.4−6 With regard to elucidation of reaction kinetics and associated assessment of steady-state bifurcation behavior in surface reactions, it should be emphasized that traditional mean-field (MF) rate equations for chemical kinetics have proved immensely valuable.2,3,7 The prescription of rates in these equations can be guided by Langmuirian formulations for adsorption and incorporates transition-state expressions for activated processes such as desorption and reaction. Selection of © XXXX American Chemical Society

activation barriers is often guided by theoretical analysis or by experiment but often incorporates an heuristic treatment of coverage dependence. However, these MF rate equations generally assume a single type of adsorption site, and it is well recognized that they neglect spatial correlations in the reactant distribution associated with ordering of the reactant adlayer. Thus, they can at best heuristically capture the local environment dependence of rates. To address these short-comings, realistic stochastic multisite lattice−gas (msLG) reaction models have been increasingly utilized over the past decade, where behavior of these models is precisely analyzed by kinetic Monte Carlo (KMC) simulations.4−6 These msLG models attempt to provide a detailed molecular-level description of overall catalytic reaction processes on metal surfaces: adsorbed reactants are located on periodic lattices representing adsorption sites; adsorption, desorption, reaction, and surface diffusion are implemented stochastically by adding and removing species from sites, or by hopping between sites, with probabilities proportional to the physical rates. As reactant adspecies can occupy distinct types of adsorption sites Received: October 6, 2016 Revised: November 29, 2016 Published: November 29, 2016 A

DOI: 10.1021/acs.jpcc.6b10102 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

2. CO OXIDATION KINETICS AND STEADY STATES: OVERVIEW If α denotes CO or oxygen (O or O2 depending on the context), let θα denote the coverage of adspecies α in monolayers (ML), with a maximum allowed value of θα(max). Below, we set θα(max) = 1/mα for an appropriate integer mα. Let pα denote the impingement rate per surface unit cell which is determined from the partial pressure, Pα, via pα ≈ ΩPα/(2πMαkTg)1/2 where Ω ∝ a2 is the unit cell area for surface lattice constant a, Mα is the mass of the adsorbing species, and Tg is the gas temperature. Thus, one has pα ≈ 105Pα at 300 K where the units of pα are ML per sec and Pα is in mbar. We consider the low-pressure regime where pα = O(1) ML/s. Catalytic CO oxidation involves adsorption of CO, dissociative adsorption of O2, nonreactive desorption of CO, and reaction of CO and O and is described by the rate equations2,3,6

(hence the multisite aspect of the models), a separate lattice is required for each site type. Also, the mixed reactant adlayers can display complex superlattice ordering primarily due to strong short-range repulsive lateral interactions between adspecies. At lower pressures (which correspond to lower adspecies coverages), the high effective adspecies mobility means that these adlayers generally exhibit near-equilibrium ordering. However, this high mobility also makes KMC simulations computationally demanding.4−6 DFT analysis can facilitate a detailed description of system thermodynamics in the msLG models (assessing site-specific adsorption energies, as well as adspecies interactions), and thus of adlayer ordering.4−6 However, limitations to DFT energetics mean that development of higher-level techniques and also input from experiment are valuable. Furthermore, we emphasize that the TOF and other features of reactive steady states are determined not just by thermodynamic parameters but also by the detailed kinetics of adsorption, desorption, and reaction, which depends on the local mixed reactant adlayer environment.6 Again, DFT analysis can supplement the typically limited information available from experiment on kinetics. Despite the increasing potential to develop detailed and realistic msLG models, we remark that their complexity and uncertainties in the thermodynamics and kinetics can limit insight into which features or parameters control key aspects of reaction behavior. Thus, an appealing alternative is to develop somewhat simplified or tailored msLG reaction models, although these should still incorporate essential features of the catalytic reaction process. We note that there exists some limited tailored msLG modeling for CO oxidation on Pd(100) by our group,8,9 and also by Reuter and co-workers.10−12 We emphasize that KMC simulation of these models provides the most precise assessment of behavior. Indeed, such simulation will be used to provide benchmark results in this study. However, our focus is on the development of a predictive beyond-MF analytic rate equation treatment motivated by the following observations. First, such analytic modeling can provide deeper insight into reaction behavior than comes from “black-box” KMC simulations. Indeed, reliable analytic treatment and ideally exact solution of statistical mechanical models has been a longstanding goal for equilibrium systems,13 and the current work extends this philosophy to nonequilibrium reaction systems. Second, analytic treatments are far more computationally efficient than KMC simulation. This is particularly the case given the extremely high hop rates for reactant adspecies (especially CO) under reaction conditions of interest.6 In fact, it is not even viable due to computational expense to perform simulations for physical high values of these rates, so typically lower values are used assuming that this does not impact behavior. Third, availability of an analytic formulation enables implementation of so-called analytic continuation methods to efficiently map out the entire bifurcation diagram for steady-state reaction behavior.2,3,14 The content of subsequent sections is as follows. In section 2, we review traditional rate equation descriptions for CO oxidation kinetics and also basic concepts of bistability and multistability. In section 3, we present analysis of key energetics for the system. Our tailored models for CO oxidation on Pd(100) are described in section 4. An analytic description of the associated kinetics is developed in section 5. We compare predictions of our beyondMF analytic treatment with KMC simulation results in section 6 and also demonstrate the limitations of traditional MF treatments. Conclusions are provided in section 7.

d/dtθCO = pCO SCO − dCOθCO − R CO + O d/dtθO = 2pO SO2 − R CO + O 2

and (1)

Here SCO (SO2) denote the sticking probabilities for CO (O2), dCO denotes the CO desorption rate, and RCO+O denotes the reaction rate. In principal, these equations are exact if sticking coefficients and desorption and reaction rates account for the dependence of these rates on local environment and on adlayer ordering as prescribed in the msLG model. The equations are only closed when both CO and O are highly mobile so that the adlayer is near-equilibrated. Then, Sα, dCO, and RCO+O are completely but nontrivially determined by the coverages and surface temperature, T. As noted above, standard MF rate equation expressions for rates fail to capture this nontrivial coverage dependence. We consider only this regime of equilibrated adlayers, which should apply for low pα and typical T from 400 to 600 K due to high adspecies mobility. In addition to the TOF, a central goal of rate equation modeling is to predict the variation with Pα (or pα) and T of the reactant coverages in the steady state(s), i.e., to describe the “reactive equation-of-state”. For systems exhibiting bistability or multistability, a key additional goal is to predict the steady-state bifurcation diagram usually upon varying one partial pressure and T.2,3,6 For our models, the form of the reaction rate does not differ from the classic MF form RCO+O = zkθCOθO, where k is the microscopic reaction rate for a CO−O pair in the appropriate configuration and z is a coordination number. Thus, the critical factor in determining the reactive equation-of-state is the nontrivial dependence of the Sα on coverages. Traditional treatments assume a single type of adsorption site (and below we set Sα = 1 for a clean surface). For O2 adsorption, which occurs on a pair of empty adsorption sites, a classic Langmuirian form for SO2(lang) = (1 − mCOθCO − mOθO)2 is adopted.7 For CO adsorption, one might also adopt the classic Langmuirian form, SCO(lang) = 1 − mCOθCO − mOθO, reflecting the feature that CO requires a single empty site.7 Alternatively, one could consider an asymmetric form, SCO(asym) = 1 − (mCOθCO)q, reflecting the experimental observation that adsorbed O does not strongly inhibit CO adsorption. In addition to the natural choice q = 1, one might also select q = 2 or 3 on the basis of a Gasser−Smith precursor theory.15 We caution that none of these simplistic forms for Sα will apply to realistic msLG reaction models. Next, we describe MF-type bistability in CO oxidation reactions for small dCO wherein a reactive state (with high θO B


Article

The Journal of Physical Chemistry C and low θCO) and a near-CO-poisoned state (with high θCO and low θO) coexist for a range of pCO at fixed pO2. This range is terminated by upper and lower spinodals, pCO±. In particular, our main goal here is to clarify how bistability reflects or is dependent on specific behavior of the Sα. Existence and behavior of the reactive (R) steady state, particularly for small pCO, is most easily understood by regarding CO desorption as negligible in this regime. Then there is a balance between adsorption and reaction, i.e., pCOSCO ≈ 2pO2SO2 ≈ RCO+O. Thus, it is clear that SO2 must decrease faster to zero than SCO as pCO → 0 and θO → θO(max) to balance the adsorption rates for CO and O2. Imposing the above rate balances, for small pCO, one finds that θCO ≈ pCO ( 1 2 pCO /pO2 )n /[zkθO(max)]

adsorption energetics on extended metal surfaces uses the planewave code VASP with a periodic slab geometry.19 Our analysis uses the projector augmented-wave (PAW) formulation20,21 with the Methfessel−Paxton method.22 This analysis is typically performed with paired electron spin (i.e., for the singlet state) noting that bulk Pd is nonmagnetic. For more details, see the Supporting Information. We perform such an analysis using the PBE functional23 for a (2√2×2√2)R45° unit cell corresponding to CO coverage of 1/8 ML. We also average over Pd slab thickness from 7 to 12 layers to minimize quantum size effects for finite slab thickness.24 We find that EaCO = −1.94, −1.88, and −1.51 eV for br, 4fh, and top sites, respectively. This analysis does correctly predict the preference for the br site in contrast to an earlier study using a smaller unit cell that predicted a slight preference for the 4fh site.25 We do note, however, that this PBE prediction of EaCO = −1.94 eV for binding at the preferred br adsorption site significantly overestimates the true value. We take the true value to be EaCO ≈ −1.64 eV based on a correction to the DFT analysis described below,25 and also based on modeling of CO temperature-programmed desorption spectra.6,26 As an aside, we have also performed limited plane wave calculations for a slab geometry using the hybrid HSE06 functional.27 Given the expense of such analysis in VASP, we restrict attention to thin slabs using a small c(2×2) unit cell with sparse k points and use frozen PBE optimized geometries. We obtain br site adsorption energies close to PBE values, e.g., EaCO = −2.10 eV (−1.92 eV) for a 4-layer slab and EaCO = −1.89 eV (−1.90 eV) for a 5-layer slab for HSE06 (PBE). 3.2. CO Adsorption Energy: Paired Spin Analysis for Cluster Geometries. Additional analyses of EaCO using higherlevel methods, and also at the DFT level for various hybrid functionals, naturally utilize localized atomic orbitals basis sets. Such analyses are generally performed on clusters rather than on slab geometries. Thus, for comparison with such analyses presented below, we have also performed a plane-wave VASP study using the PBE functional for CO adsorption on Pd clusters with one cluster and a single adsorbed CO per unit cell, and also initially considering the singlet state. This analysis uses Gaussian smearing with σ = 0.2 eV as the default. The choice of cluster geometry is described in the Supporting Information and also below. These analyses predict a preference for the br site in the limit of large cluster size, consistent with the above analysis for slab geometries. Furthermore, we find that limiting large-size values of EaCO for br and top sites are achieved relatively quickly with increasing cluster size compared to that for 4fh site. See the Supporting Information. This observation supports the value of further analysis of the EaCO for br sites based on relatively small Pd clusters, noting again that EaCO is a key parameter for our msLG models. Specifically, to mimic br sites on extended (100) surfaces with a relatively small cluster, we consider primarily an 8 atom Pd cluster with 3 × 2 = 6 atoms in the top layer and 2 atoms beneath, as shown in Figure 1a. One could also consider a 20 atom Pd cluster with a 3 × 4 layer of atoms on top of those in the 8 atom cluster, a 40 atom Pd cluster with a 4 × 5 layer of atoms on top of those in the 20 atom cluster, etc., as illustrated in the Supporting Information. These clusters are constructed to have a br site exactly in the center of the top layer at which the adsorbed CO is always located. We now discuss in some detail several analyses of EaCO for the 8-atom Pd cluster. For reference, our plane-wave PBE estimate for CO binding at br sites is EaCO ≈ −2.04 eV (versus −1.94 eV for the slab geometry), and plane-wave DFT analysis using the

(2)

vanishing as pCO → 0, where n = 0 (1) for Langmuirian (asymmetric) SCO. Existence of the near CO-poisoned (P) steady state for small dCO requires that SO2 decreases faster to zero than SCO as θCO → θCO(max). This feature produces an adsorption−desorption equilibrium of CO with θO ≈ 0. Because SO2 ≪ SCO as θCO → θCO(max), where θO ≈ 0, a rough balance between CO adsorption and desorption yields θCO(max) − θCO ≈ dCO[θCO(max)]2 /(n′pCO )

(3)

vanishing as dCO → 0, where n′ = 1 (q) for Langmuirian (asymmetric) SCO. A more comprehensive analysis reveals the presence of an unstable (U) steady state connecting the reactive and near COpoisoned steady states (cf. results in section 6). This analysis also shows that as T and thus dCO increases, the bistability region narrows and disappears, i.e., pCO± merge, at a nonequilibrium critical point,3,6 Tc. This critical point is in the mean-field universality class for high adspecies mobility.16 Above Tc, there is a unique steady state with θCO increasing monotonically with pCO. Finally, we note that MF-type bistability and criticality are expected to be preserved for more realistic msLG models with high surface mobility at least of CO.6,17,18 See section 6. However, we shall see that MF bistability can be lost if the adsorption mechanisms and associated Sα are not suitably prescribed. Behavior can also be complicated by the feature that the Sα can adopt distinct forms in extreme regimes of coverage close to maximum values due to symmetry-breaking order− disorder transitions in the adlayer. This feature can induce more complex multistability of steady states.

3. KEY ENERGETICS FOR CO OXIDATION Two key energies controlling behavior in the CO oxidation reaction are the barrier, Edes, for CO desorption which determines the desorption rate, dCO, and the barrier, Erxn, for the CO + O reaction which controls the reaction rate, k. In particular, Edes or dCO will be a particularly important parameter in our modeling as it determines the critical temperature, Tc, in the bifurcation diagram for the CO oxidation reaction. Our modeling will assume that CO adsorption is not activated, so then Edes corresponds to the CO adsorption energy, EaCO. Thus, in section 3.1 and 3.2, we review previous analysis of EaCO and perform additional studies. In section 3.3, we discuss analysis of Erxn. 3.1. CO Adsorption Energy: Paired Spin Analysis for Slab Geometries. The standard strategy for DFT analysis of C


Article


Supporting Information presents additional analyses for localized atomic orbital basis sets. 3.3. CO Adsorption Energy: Unpaired Spin Analysis. It is well recognized that, although the macroscopic bulk crystalline form of certain transition metals is nonmagnetic, nanosized clusters35,36 or strained heteroepitaxial thin films37 of these same metals can become magnetic. This is the case for Pd. One can justify the paired spin singlet-state analysis of section 3.1 and 3.2 because we are interested in behavior in the limit of large clusters mimicking the extended surface. This limiting behavior should be the same for both paired and unpaired spin analyses given that bulk Pd is nonmagnetic (and assuming that the surface layers of a macroscopic Pd(100) crystal are nonmagnetic). Nonetheless, it is instructive to present selected results of an analysis with unpaired spin for small Pd clusters to compare with the above paired spin calculations. Plane-wave PBE using VASP predicts EaCO ≈ −1.87 eV (versus −2.04 eV) for the analysis based on unpaired (paired) spins for CO + 8Pd, and EaCO ≈ −2.23 eV (versus −2.14 eV) for unpaired (paired) spins for CO + 20Pd. These calculations also reveal a magnetic moment, μ (reported in units of 1/2 spin) of 4.18 (2.23) for 8Pd (CO + 8Pd), and 9.99 (9.78) for 20Pd (CO + 20Pd). Using the same Gaussian basis set as in section 3.4 for a PBE calculation, one obtains EaCO ≈ −1.55 eV (versus −1.98 eV) for unpaired (paired) spin analysis for CO + 8Pd, and EaCO ≈ −2.02 eV (versus −1.87 eV) for unpaired (paired) spins for CO + 20Pd. These calculations reveal spin multiplicities M = 5 (3) for the ground state of 8Pd (CO + 8Pd), and M = 7 (11) for the ground state of 20Pd (CO + 20Pd). These predictions of spin multiplicity are reasonably consistent with the DFT VASP results noting again that VASP reports the total spin, S, in units of 1/2 spin, so that S = μ/2 and M = 2S + 1 = μ + 1. Again, we expect that results for EaCO for unpaired and paired spins will merge in the limit of large cluster size. 3.4. CO + O Reaction Barrier. There have been multiple previous estimates based on DFT of the activation barrier, Erxn(√5/2) for the CO + O reaction on Pd(100) starting with a br CO and 4fh O separated by d = √5a/2.10,38−40 These estimates range from Erxn(√5/2) = 0.75 to 1.05 eV. Our own DFT PBE analysis gives estimates at the lower end of this range, which we believe are lower than the true value. We attribute this result to the feature that such a DFT PBE analysis6 also appears to underestimate the barrier for O diffusion.41,42 This is relevant because O diffusion corresponds to moving the O atom from a 4fh to a br site, which also occurs in the initial stage of the CO + O reaction process. We also remark that previous stochastic modeling26 successfully recovered temperature-programmed reaction spectra for CO oxidation on Pd(100) selecting a barrier of Erxn(√5/2) = 1.0 eV. As an aside, we have performed a VASP analysis which predicts an O diffusion barrier on Pd(100) of 0.23 eV from PBE versus 0.47 eV from the meta-GGA revTPSS functional43 (from a unit cell with θO = 0.2 ML). These values should be compared with the rough experimental estimate of ∼0.6 eV.41,42 Thus, we anticipate that revTPSS will give a more accurate estimate of the reaction barrier than PBE but do not pursue that analysis in the current study.

Figure 1. (a) Top-view of the 8-atom Pd cluster. (b) CO br site adsorption energy, EaCO, versus CO singlet−triplet energy split, ΔET‑S, for various analyses. The experimental ΔET‑S (and the corresponding EaCO) are indicated by the blue dashed lines.

HSE06 functional obtained EaCO ≈ −1.99 eV (versus −1.89 eV for the slab geometry). Thus, the HSE06 values are slightly weaker than PBE values. Next, we present the results of additional analyses utilizing localized atomic orbital bases and various methods. These were performed using primarily NWChem28 (but with some CR-CC(2,3) calculations using GAMESS29,30). We use a LANL2DZ basis set with effective core potentials (ECP)31 for Pd, this choice being effectively utilized for other transition metal clusters,32 and the 6-311++G(d,p) basis set33 for CO. Various levels of theory are considered including Hartree−Fock (HF), PBE,23 and PBE034 DFT functionals, and CR-CC(2,3). In all cases, the cluster geometry is fixed to that obtained from the PBE VASP analysis. We obtain EaCO ≈ −3.19 eV for HF, −1.98 eV for PBE (quite consistent with the VASP value), −1.89 eV for PBE0 (close to the VASP value from the similar HSE06 functional), and −1.01 eV for CRCC(2,3). For each of these treatments, we have also determined the singlet−triplet energy splitting, ΔET‑S, for CO yielding 5.417 eV for HF, 5.905 eV for PBE, 5.894 eV for PBE0, and 6.207 eV for CR-CC(2,3). These should be compared with the experimental value of ΔET‑S = 6.095 eV.25 In Figure 1b, we plot the predicted EaCO versus ΔET‑S for these various methods showing a strong linear correlation. This correlation has been noted previously and used to correct PBE predictions for EaCO.25 However, the previous study explored this correlation only for the PBE functional with different choices of pseudopotentials, as opposed to our study using different levels of theory. One reasonable strategy25 to estimate the actual CO adsorption energy is to use the above linear relationship to extract the value of EaCO corresponding to the experimental ΔET‑S. Applying this strategy suggests that EaCO ≈ −1.5 eV for the 8 atom Pd cluster, significantly below the PBE and PBE0 estimates. This assessment is consistent with our proposal that the value EaCO ≈ −1.64 eV for the extended Pd(100) surface is significantly below the PBE slab geometry estimate of EaCO ≈ −1.94 eV. There is, however, a caveat with regard to the reliability of the above results obtained using localized atomic orbital basis sets. CO + 8Pd has many nearly degenerate energy states, so the system can potentially get trapped in an excited state during energy minimization, a feature that would corrupt ground-state energy estimates. This behavior may become more prevalent for larger clusters. We have identified this potential problem by comparing GAMESS and NWChem results, where these codes generally have somewhat different self-consistent field (SCF) algorithms and convergence criteria. For plane-wave DFT VASP analysis, the extent of this problem appears to be reduced, perhaps ameliorated by smearing and partial occupancy of states.

4. TAILORED MODELS FOR CO OXIDATION ON PD(100) 4.1. Overview of Tailored Models. The key ingredients of our tailored msLG models for (CO + O)/Pd(100) are as follows (see Figure 2 for relevant features of the surface geometry): D


Article


influence of finite CO−O lateral interactions. Finally, we comment that any prescription of CO adsorption−desorption must be consistent with the prescribed adlayer thermodynamics through detailed-balance constraints. Sometimes desorption rates are written as the product of a coverage-dependent sticking probability times a factor dependent on adlayer thermodynamics.48 However, in the kinetic prescription of our model, it is more natural to directly specify the (constant) CO desorption rate, so that all coverage dependence is incorporated into the sticking probabilities (still consistent with detailed balance). 4.2. CO Adsorption Dynamics. First, we describe a conventional Langmuir-type model based on the idea that CO adsorption requires a single available br site. Then, we describe a revised model more compatible with experimental CO sticking behavior. Langmuir-Type Model. For pure CO adlayers, this model specifies that CO adsorption occurs when the selected br site, as well as all 4 NN br sites separated by d = a/√2 and all 4 2NN br sites separated by d = a, are unoccupied by CO. For mixed reactant adlayers, one specifies that in addition to the above 9 br site ensemble being unoccupied, both 4fh sites separated from the selected br site by d = a/2 are free of O. See Figure 3a. This model was recently utilized by Reuter and co-workers.10−12

Figure 2. Adsorption site geometry on a fcc(100) surface. NN br CO pairs, vicinal and geminal 2NN br CO pairs, NN 4fh O pairs, and NN 4fh-br CO−O pairs are excluded. Gray circles (Pd surface atoms) will be omitted in some of the following schematics, but these schematics retain the square grid.

(i) CO resides exclusively on bridge (br or b) sites. (ii) O resides exclusively on 4-fold hollow (4fh or h) sites. (iii) Adspecies pairs with separation d ≤ a are excluded, so that nearest-neighbor (NN) 4fh O (with d = a), NN and 2NN br CO (with d = a/√2 and d = a, respectively), and NN 4fh-br O−CO (with d = a/2) configurations are blocked. This implies that θα(max) = 1/2ML (mα = 2) for both α = CO and O. “Available sites” are defined as those consistent with these blocking constraints. (iv) No finite lateral adspecies interactions are included for separations d > a (where some caveats are discussed below). (v) Dissociative adsorption of O2 and nondissociative adsorption of CO occur according to rules depending on the local adlayer environment and consistent with the exclusion rules (iii) described below. There are two br sites per unit cell, each with a CO impingement rate of 1/2pCO. There are two orientations for the dissociating O2 per unit cell, each with an impingement rate of 1/2pO2. (vi) CO desorbs nonreactively from br sites with rate dCO. (vii) Each 2NN br-4fh CO-O pair with separation d = √5a/2 reacts at rate k to form the product CO2, which immediately desorbs from the surface. (viii) Extremely rapid hopping of CO occurs to available nearest-neighbor (NN) br sites (a distance d = a/√2) with rate hCO. Rapid hopping of O occurs to available NN 4fh sites (a distance d = a) at rate hO. There have been comprehensive DFT-based studies of lateral interactions for CO/Pd(100)6,26 and O/Pd(100)6,26,44 and for the mixed system6,26 (as well as earlier estimates from statistical mechanical modeling of experimental data for the single component systems42,45). There do exist finite interactions for d > a. However, the above prescription of adsorption sites and short-range exclusion suffices to recover key features of experimentally observed adlayer ordering in this system: c(2×2) ordering of O at 4fh sites (but not more subtle p(2×2) ordering);46 c(2×2) and local c(2×2√2)R45° ordering of CO at br sites (but not long-range c(2×2√2)R45° ordering).47 As indicated in section 2, a critical ingredient in modeling is the prescription of adsorption onto mixed reactant adlayers which accounts for the “ensemble” of available sites required for adsorption to occur. Thus, a detailed prescription is provided below for CO in section 4.2 and for O2 in section 4.3. Again, we mention that our tailored model incorporates two key activation barriers, Edes for CO desorption and Erxn for CO + O reaction. In section 4.4, we discuss further treatment of these thermally activated processes as well as diffusion and also comment on the

Figure 3. CO adsorption dynamics: (a) Langmuir model where red (blue) open circles indicate br (4fh) site required to be empty for adsorption; (b) schematic of the steering + funneling model. Dashed arrows indicate br sites at which CO cannot adsorb (funneling is blocked). Solid arrows indicate available sites, where CO selects the shorter funneling path.

Steering to Top Sites and Subsequent Funneling. For CO adsorption on pure CO adlayers, the experimental SCO versus θCO is initially flat at 350 K,47 a feature usually attributed to precursor effects. However, motivated by a DFT analysis of adsorption energetics,6,49 to capture this feature we adopt an alternative picture involving steering of CO to (unpopulated) top sites and subsequent funneling to nearby available br sites.6,26 Specifically, in our model, CO is steered to top sites49 and then funnels at first a distance d = a/2 to one of the four neighboring br sites should they be available.6 Without further refinement of this model, just one preadsorbed CO can block adsorption resulting in a linear decrease in SCO with θCO. Thus, to obtain the above-mentioned initially flat SCO versus θCO, we also allow funneling to one of the eight br sites a distance d = √5a/2 away to attempt to find an available adsorption site should none of the d = a/2 br sites be available.6 See Figure 3b. 4.3. Dissociative Oxygen Adsorption Dynamics. The classic Langmuir model, which assumes dissociative adsorption on NN pairs of (4fh) adsorption sites, is not compatible with the exclusion of NN 4fh O pairs in our model. In fact, experimental studies of oxygen sticking on Pd(100) indicate that SO2 saturates around 0.34 ML at 200 K and 0.38 at 300 K due to ensemble or blocking effects not included in the classic model.41,46 This feature and NN site exclusion must be incorporated into any E


Article


νr exp[−Erxn/(kBT)] where νr = 1013/s and Erxn = 1.00 eV.6 It is appropriate to note that previous DFT analysis revealed a significant d = √5a/2 br CO−4fh O repulsion, ωCO−O, above 0.1 eV. Then, at least for lower reactant coverages, the effective reaction barrier is given by Erxn = Erxn(√5/2) + ωCO−O as the diffusing CO must surmount the repulsive barrier to reach the initial d = √5a/2 br CO−4fh O reaction configuration.6,26 In modeling including the repulsion ωCO−O, the CO desorption barrier at br sites separated by d = √5a/2 from 4fh O would also be lowered. However, the CO population of such sites would also be reduced from other sites, thus counterbalancing the effect of the lower barrier. Finally, we note that barriers for thermally activated diffusion are expected to be low, resulting in very high hop rates (especially for CO). Thus, to make KMC simulation computationally viable, we cap these rates below their true values but still well above other rates, which should suffice to achieve adlayer equilibration. Specifically, we choose hCO ∼ 104/s and hO ∼ 102/s.

realistic model. In addition, analysis of experimental data for the coverage dependence of sticking indicates that SO2 ∝ 1 − cθO, for low O-coverage θO, with c ≈ 3−4 at around 200 K.50 BBB 8-Site Model Adapted to Mixed Adlayers. In the Brundle−Behm−Barker (BBB) 8-site model for oxygen adsorption on pure O adlayers, O2 adsorbs on second NN (2NN) empty 4fh sites provided 6 additional 4fh neighbors are also empty.51,52 Thus, an ensemble of 8 empty 4fh sites is required for adsorption. This model was originally developed for O2 adsorption on Ni(100), as highlighted in Zangwill’s monograph.52 However, it has been adopted and applied extensively to describe O2 adsorption on Pd(100).6,10,11,26,41,42,53 This 8-site model is extended to mixed CO + O adlayers on Pd(100) by demanding that the 8 br sites within a distance d = a/2 of the target 2NN 4fh sites are free of CO (as well as the 8 4fh sites being free of O).10 See Figure 4a.

5. ANALYTIC EXPRESSIONS FOR KINETICS IN TAILORED MODELS Our goal in this section is to provide reliable analytic expressions for SCO, SO2, dCO, and RCO+O appropriate for our msLG models that are applicable for all coverages, θCO and θO. These expressions can then be incorporated into the rate eqs 1 to provide a reliable analytic description of reaction kinetics and steady-state behavior. The success of this analytic approach will be demonstrated in section 6, where we compare analytic predictions with precise results for model behavior from KMC simulation. Here, we first define a few key quantities and discuss relationships between them. The probability that a specific br site is occupied by CO is denoted by θCObr. Then, the CO coverage is given by θCO = 2θCObr because there are 2 br sites per surface unit cell. It is instructive to also consider various compact motifs of n br sites in which at most one br site can be occupied by CO due to short-range CO−CO exclusion. Then, the probability Pnb that all n br sites in the motif are unoccupied or empty is given exactly by Pnb = 1 − nθCObr. This applies, e.g., for single br sites (n = 1), for pairs separated by d = a/√2 (n = 2) or d = a (n = 2′), triangular trios (n = 3), and diamond-shaped quartets with sidelength d = a/√2 (n = 4). See the Appendix. The O coverage, θO, equals the probability that a specific 4fh site is occupied. Because adjacent 4fh sites cannot be occupied by O, it follows that the probability of an adjacent empty NN 4fh pair (Figure 5c) is given exactly by P2h = 1 − 2θO. See the Appendix. Also in the Appendix, we provide a detailed discussion of related conditional probabilities which will be utilized below. 5.1. Langmuir-Type Model for CO-Adsorption from Section 4.2. First, considering adsorption on a pure CO adlayer,

Figure 4. (a) 8-site and (b) 9-site models for dissociative adsorption of oxygen adapted to treat adsorption on mixed reactant adlayers. Open circles indicate empty sites.

9-Site Model Adapted to Mixed Adlayers. Despite the prominence of the 8-site model for describing oxygen adsorption on pure O adlayers on metal(100) surfaces, DFT analysis actually suggests a different picture. Dissociative adsorption tends to occur via O2 impinging on a 2NN vicinal br pair with d = a and then separating to third NN (3NN) 4fh sites54,55 (i.e., adsorption does not occur directly onto 2NN 4fh sites as prescribed in the 8site model). This observation motivates adoption of an alternative 9-site model50 for O2 adsorption on pure O adlayers wherein the O atoms from the dissociating molecule are placed on an empty 3NN 4fh pair, provided that the 7 additional NN 4fh sites are also empty.6,50 Extending the 9-site model to treat mixed CO + O adlayers on Pd(100), in addition we demand that the 8 br sites within a distance d = a/2 of the target 3NN 4fh sites are also free of CO. See Figure 4b. Refined 8-Site and 9-Site Models with Reorientation. It was recently recognized that the initial rate of decrease of SO2 with θO in the above 8-site or 9-site models is significantly greater than the experimentally observed behavior mentioned above.50 Thus, we are motivated to consider “more flexible” adsorption rules to boost SO2 at least for lower θO. One possibility is to allow reorientation of the adsorbing molecule. For example, O2 impinging parallel to the surface with its center of mass above a top site for the 8-site model or a 4fh site for the 9-site model could sample two orthogonal orientations to assess if one of these provides an available adsorption site ensemble. If so, O2 adsorption then randomly selects from available orientations. 4.4. Thermally Activated Processes: CO Desorption, CO + O Reaction, CO and O Diffusion. CO desorption, CO + O reaction, and adspecies diffusion are all thermally activated processes, and their rates are taken to have an Arrhenius form. Given the lack of finite interactions in our model, all CO molecules desorb with the same rate, dCO = νd exp[−Edes/(kBT)], where we select νd = 1016/s and Edes = −EaCO = 1.64 eV.6 Each br CO−4fh O pair with separation d = √5a/2 reacts with rate k =

Figure 5. Langmuir-type CO adsorption on the central br site. Open red (blue) circles denote empty br (4fh) sites. (a) SCO = P9b for a pure CO adlayer. (b) SCO = P9b+2h for a mixed adlayer. F


Article

The Journal of Physical Chemistry C denote the probability of finding the configuration of 9 empty br sites required for adsorption (Figure 3a and also Figure 5a) as P9b. Thus, one has SCO = P9b. Then, extending Kirkwood-type approximations, we factorize the probabilities for multisite configurations as a product of probabilities for natural constituent subunits selected as the largest compact motifs of br sites, no more than one of which can be occupied. We also compensate for overcounting due to overlapping or shared sites between multiple motifs. This ensemble of 9 empty br sites is composed of 4 empty diamonds each of 4 br sites, where adjacent diamonds share a diagonal empty NN br pair. All of these motifs share a central empty br site, as shown in Figure 5a. Recalling that Pnb = 1 − nθCObr gives the probability for a compact ensemble of n empty br sites, one writes SCO = P9b = (P4b)4 P1b/(P2b)4 = (1 − 2θCO)4 (1 − 1 2 θCO)/(1 − θCO)4

Figure 6. Configurations blocking adsorption for CO adsorption with steering to top sites and subsequent funneling to br sites a distance d = a/2 or √5a/2 away. Red (blue) solid circles denote CO (O) adsorbed on br (4fh) sites. In the configurations with 3 CO, smaller red circles indicate alternative blocking positions for the third CO. We indicate the number of distinct configurations accounting for symmetry.

An extra quartic factor with parameter α is included to better control behavior when SCO vanishes at θCO = 1/2 ML. We select α = 16 to impose nonlinear decay as θCO → 1/2 ML. 5.3. BBB 8-Site Model for Oxygen Adsorption from Section 4.3. First, for pure O adlayers, we denote the probability of finding the configuration of 8 empty br sites required for adsorption (Figure 4a and Figure 7a) as P8h. Then, one has SO2 =

(4)

For mixed reactant adlayers, we denote the probability of finding the 9 empty br sites and the 2 additional 4fh sites configuration required for adsorption (Figure 3a and Figure 5b) as P9b+2h. Now, one has SCO = P9b+2h and it is instructive to write SCO = P9b + 2h = P9bQ 2h | 9b

(5)

where Q2h|9b is the conditional probability that the central pair of 4fh sites in Figure 5b is unoccupied by O given that the 9 br sites are unoccupied by CO. We adopt the approximation Q2h|9b ≈ 1 − 2qO where qO = θO/(1 − 4θCObr) is the conditional coverage for 4fh site to be populated by O given that the 4 NN br sites are free of CO. See the Appendix. Thus, we conclude that SCO = [(1 − 2θCO − 2θO)(1 − 2θCO)3 (1 − 1 2 θCO)]/(1 − θCO)4 (6)

This expression satisfies SCO ≈ 1 − 9θCObr − 2θO ≈ 1 − 4.5θCO − 2θO, for low coverages, and SCO appropriately vanishes when θCO + θO = 1/2. 5.2. Steering + Funneling Model for CO-Adsorption from Section 4.2. First, we determine various configurations with a minimal number of preadsorbed species which block CO adsorption according to the prescription of section 4.2, and thereby develop a low-coverage expansion for SCO. Just 1 or 2 preadsorbed CO molecules cannot block adsorption. However, there exist 76 distinct configurations with 3 preadsorbed CO molecules which block adsorption. Likewise, there are 4 configurations with 1 preadsorbed O and 2 preadsorbed CO molecules which block adsorption, and 2 configurations with 6 preadsorbed O atoms. See Figure 6. Consequently, we conclude that for low coverages

Figure 7. Refined 8-site (a,b) and 9-site (c,d) models for O2 adsorption. Open red (blue) circles denote empty br (4fh) sites. The diagonal ovals indicate the orientation of the adsorbing O2. Schematic of configurations enabling adsorption: (a, c) pure O adlayers; (b, d) mixed reactant adlayers.

P8h. This 8-site ensemble can be regarded as composed of 8 empty pairs of NN 4fh sites, where each two overlapping pairs share a single site. Recall that Pnh = 1 − nθO gives the probability of an empty site (pair) for n = 1 (n = 2). Thus, utilizing a standard Kirkwood approximation, one writes50 SO2 = P8h ≈ (P2h)8 /(P1h)8 = (1 − 2θO)8 /(1 − θO)8

SCO ≈ 1 − 76(θCO br)3 − 4(θCO br)2 θO − 2(θO)6 3

2

6

= 1 − 9.5(θCO) − (θCO) θO − 2(θO)

(9)

For mixed CO + O adlayers, we denote the probability of finding the 8 empty 4fh sites and the additional 8 empty br sites (Figure 4a and Figure 7b) as P8h+8b, so now SO2 = P8h+8b. It is instructive to write

(7)

A more useful and global expression for SCO is obtained by converting the above Taylor expansion into a Padé-type approximant56 which not only retains the low-coverage behavior but also ensures that SCO vanishes at θCO + θO = 1/2 ML.9 Thus, one writes

SO2 = P8h + 8b = P8hQ 8b | 8h

(10)

where Q8b|8h is the conditional probability that the eight br sites in Figure 7b are unoccupied by CO given that the set of 8 4fh sites are unoccupied by O. This ensemble of 8 empty br sites can be regarded as composed of 3 diamonds of 4 br sites each of which is unoccupied with a conditional probability denoted by Q4b. Adjacent diamonds share an overlapping diagonal pair of br sites

SCO = [1 − 2θCO − 2θO][1 − α(θCO)4 ] /[1 − 2θCO − 2θO + 9.5(θCO)3 + (θCO)2 θO + 2(θO)6 ] (8) G


Article

The Journal of Physical Chemistry C that is unoccupied with conditional probability denoted by Q2b. Thus, one writes Q 8b | 8h ≈ (Q 4b)3 /(Q 2b)2

orientations are available requires an ensemble of 12 (13) empty 4fh sites and 12 (16) empty br sites, so we write P(1∧2) = P12h+12b (P13h+16b). The latter quantities can be factorized analogous to the former. See the Appendix for a detailed analysis. For the 8-site model with reorientation, one finds that SO2 = 2P8h+8b − P12h+12b ≈ 1−4θO − 4θCObr ≈ 1 − 4θO − 2θCO, for low coverages. For the 9-site model, one has that SO2 = 2P9h+8b − P13h+16b ≈ 1 − 5θO − O(θCO2) for low coverages reflecting the feature that two CO are required to block adsorption. Thus, introducing reorientation produces the desired slower initial decrease in S O 2 with θ O (relative to models without reorientation). 5.6. CO Desorption and CO + O Reaction. As noted in section 4.4, all CO molecules desorb with the same rate, dCO = νd exp[−Edes/(kBT)], where νd = 1016/s, and Edes = 1.64 eV.6 Thus, the rate of loss of CO due to desorption is simply dCOθCO with this environment-independent dCO [cf. eq 1]. As also noted in section 4.4, for each CO, reaction is possible with an O sitting on one of four 4fh sites a distance d = √5a/2 away with rate k = νr exp[−Erxn/(kBT)] where νr = 1013/s and Erxn = 1.00 eV.6 Accounting for all possible local environments of the CO and for the associated differing reaction rates (Figure 8), as well as for the

(11)

We adopt the approximation Qnb ≈ 1 − nqCO where qCObr = θCObr/(1 − 2θO) is the conditional coverage for br site to be populated by CO given that the 2 NN 4fh sites are free of O. See the Appendix. Thus, we conclude that br

SO2 = P8h + 8b ≈ [(1 − 2θO)7 (1 − 2θCO − 2θO)3 ] /[(1 − θO)8 (1 − θCO − 2θO)2 ]

(12)

This expression satisfies SO2 ≈ 1 − 8θO − 8θCO ≈ 1 − 8θO − 4θCO, for low coverages, and SO2 = 0 when θCO + θO = 1/2. 5.4. 9-Site Model for Oxygen Adsorption from Section 4.3. For pure O adlayers, adsorption requires an ensemble of 9 empty 4fh sites which occurs with probability denoted P9h. See Figure 4b and Figure 7c. This 9-site ensemble can be regarded as composed of 8 pairs of NN 4fh sites, so from a standard Kirkwood analysis, one writes49 br

SO2 = P9h ≈ (P2h)8 /(P1h)7 = (1 − 2θO)8 /(1 − θO)7

(13)

which satisfies SO2 ≈ 1 − 9θO for low θO and vanishes when θO = 1/2. For mixed CO + O adlayers, an additional 8 br sites must also be empty which occurs with probability denoted P9h+8b. See Figure 4b and Figure 7d. Thus, now we write SO2 = P9h + 8b = P9hQ 8b | 9h

Figure 8. Different possible local environments for a CO adsorbate on Pd(100) reacting with nearby O a distance d = √5a/2 away. Red + smaller blue (blue) solid circles denote CO (O) adsorbed on br (4fh) sites. Open blue circles denote empty 4fh sites. The overall reaction rate in (a) is k and in (b, c) is 2k.

(14)

where Q8b|9h is the conditional probability that the eight br sites in Figure 7d are unoccupied by CO given that the set of 9 4fh sites are unoccupied by O.9 This ensemble of 8 empty br sites can be regarded as composed of 2 diamonds of 4 br sites (each with probability Q4b) connected by a 2NN br pair (probability Q2b′) sharing two single br sites (probability Q1b). Then, from an analysis analogous to that for the 8-site model, it follows that Q 8b | 9h ≈ (Q 4b)2 (Q 2b ′)/(Q 1b)2

feature that there are 2 br sites per unit cell, it follows that the rate of loss of CO or O due to reaction is R CO + O = 2 × [4θCO brθOP2h × k + 4θCO br(θO)2 × 2k]

(15)

= 4kθCOθO

and SO2 = P9h + 8b

(17)

6. RESULTS OF MSLG MODELS FOR CO OXIDATION ON PD(100) We will not describe in detail the KMC simulation algorithm used below to precisely assess behavior of our tailored msLG reaction models. However, it is appropriate to mention that we implement a “constant CO-coverage” simulation algorithm (versus the standard “constant CO partial pressure” algorithm).6 This constant-coverage algorithm is particularly useful as it automatically maps out both stable and unstable steady states in the case of bistability or multistability (whereas unstable states are not readily assessed with the standard simulation algorithm).6,17,18 6.1. Model of Hoffmann et al. (Langmuir-Type CO and 8-Site O2 Adsorption). The model of Hoffmann and Reuter10,11 for CO oxidation on Pd(100) utilizes a Langmuirtype prescription of CO adsorption (section 5.1), and the BBB 8site model adapted to mixed adlayers for dissociative O2 adsorption (section 5.3). Basic results from KMC simulation and our analytic beyond-MF rate equation treatment are shown in Figure 9 for variation of steady-state values of coverages and of

≈ [(1 − 2θO)7 (1 − 2θCO − 2θO)2 (1 − θCO − 2θO)] /[(1 − θO)7 (1 − 1 2 θCO − 2θO)2 ] (16)

See the Appendix for further discussion. Thus, one finds that SO2 ≈ 1 − 9θO − 8θCObr ≈ 1 − 9θO − 4θCO, for low coverages, and SO2 = 0 when θCO + θO = 1/2. 5.5. Refined 8- and 9-Site Models with Reorientation. To reduce the initial rate of decrease of SO2 with θO in the 8-site or 9-site models to better match experiment, we can allow reorientation of the adsorbing molecule. As described in section 4.3, the impinging O2 samples two orthogonal orientations to find an available adsorption site ensemble, and then randomly selects from available orientations. If 1 (2) denotes the case where the first (second) orientation is available, then using set notation, one has SO2 = P(1∨2) = P(1) + P(2) − P(1∧2) where ∨ (∧) denotes “or” (“and”). For the 8-site (9-site) model, one has that P(1) = P(2) = P8h+8b (P9h+8b), and the case where both H


Article


Figure 10. Behavior of sticking probabilities versus both θCO and θO for the model of Hoffman et al.10,11 Plotted are contours for ln S at increments of −0.5 decreasing from the lower left where S = 1 (and ln S = 0) to the upper right where S → 0 (and ln S → −∞) for θCO + θO → 1/2.

Figure 9. Behavior of the (CO + O)/Pd(100) model of Hoffman et al.10,11 with pO2 = 1 ML/s. Steady-state coverages for CO (a) and O (b); SCO and SO2 (c); and TOF (d) in the steady state vs θCO at 460 K. Dots (solid curves) denote KMC (analytic rate equations or RE) results.

θCO, we emphasize that our beyond-MF analytic treatment of sticking accurately reproduces precise model behavior. Again, an appealing feature of this analytic formulation is that it avoids the computational demands of KMC simulations for high surface mobility. 6.2. Refined Model (CO Steering + Funneling; 9-Site O2 Adsorption + Reorientation). This model incorporates steering and funneling for CO adsorption (section 5.2), and a 9-site model with reorientation for O2 adsorption (section 5.5), thereby capturing key features observed in experimental sticking studies. Basic results from KMC simulation and our beyond-MF rate equations are shown in Figure 11 for variation of steady-state

the TOF as a function of pCO at fixed pO2 = 1 ML/s. The corresponding variation of sticking probabilities, but plotted as a function of θCO, is also shown. Consistently in both the KMC and rate equation analysis, we find a lack of MF-type bistability for typical reaction temperatures around 400−500 K despite experimental observation of such bistability below a critical temperature of Tc ≈ 530 K.6,10,57 We note, however, that the model was developed to reasonably capture behavior for higher T around 600 K and succeeds in this respect.10,11 As indicated in section 2, the presence of bistability is controlled by the behavior of the sticking probabilities as a function of both reactant coverages. A comprehensive assessment of such behavior is rarely provided, although in principle it can be extracted from KMC simulations.6 However, our beyondMF analytic formulation provides deeper and more direct insight into relevant features. Figure 10 shows the behavior of SCO from (6) and SO2 from (12) as a function of both θCO and θO. From these plots and from the associated analytic expressions, one finds that the lack of MF bistability derives from the feature that Langmuir-type CO adsorption is too restrictive for higher θCO. Specifically, one has that SCO ∼ 192(1/2 − θCO)4 from (6) is well below SO2 ∼ 32(1/2 − θCO)3 from (12) for higher θCO close to θCO(max) = 1/2 ML (where θO ≈ 0). From section 2, this implies a lack of a near CO-poisoned state. We remark that this behavior of SCO (SO2) reflects the feature that CO (O2) adsorption requires 4 (3) empty diamonds of br sites. On the contrary, SO2 ∼ 216(1/2 − θO)8 is far below SCO ∼ 2(1/2 − θO) for higher θO (and θCO ≈ 0) consistent with the requirement from section 2 for existence of a conventional reactive steady state. Interestingly, the model does display what we characterize as non-MF bistability reflected in the multivalued θCO very close to θCO(max) = 1/2 ML. This feature derives from modified CO adsorption behavior due to symmetry breaking in the CO adlayer in this regime, behavior that is not captured in our analytic treatment. This modified behavior does ensure that SO2 ≪ SCO for θCO sufficiently close to 1/2 ML. See section 6.3 for further discussion. Apart from this non-MF bistability for near-saturation

Figure 11. Behavior of our refined (CO + O)/Pd(100) model with pO2 = 1 ML/s. Steady-state coverages for CO (a) and O (b); SCO and SO2 (c); and TOF (d) in the steady state vs θCO at 460 K. Symbols (solid curves) denote KMC (analytic RE) results. R, U, P indicate reactive, unstable, and near CO-poisoned steady states, respectively. In (a), (b), and (d), we indicate the regime of MF bistability as predicted by the analytic RE. I


Article

The Journal of Physical Chemistry C coverages and TOF versus pCO at fixed pO2 = 1 ML/s. The corresponding variation of sticking probabilities, but plotted as a function of θCO, is also shown. This model does recover the experimentally observed MF-type bistability. In fact, our model reasonably reproduces the entire experimental bifurcation diagram57 including the critical temperature for disappearance of bistability, Tc ≈ 515 K (versus Tc ≈ 535 K from analytic theory), as shown below. However, this model also displays additional non-MF multistability close to θCO(max) = 1/2 ML analogous to the model of Hoffmann et al.,10 and as discussed further in section 6.3. Apart from this non-MF multistability, our analytic treatment based on appropriate factorization approximations and Padé-type treatments of sticking accurately reproduces precise model behavior. Elucidation of the appearance of MF-type bistability in this model is readily provided by our analytic beyond-MF analysis of sticking probability behavior. Figure 12 shows the behavior of

Figure 13. (a) Bifurcation diagram showing the bistable region for ECO+O = 1.0 eV in our refined model. Dots (curves) show KMC (analytic RE) predictions. The inset shows experimental data from ref 52. (b) Analytic RE predictions for the bistable region for ECO+O = 1.1 eV (yellow), 1.0 eV (black), 0.9 eV (purple), and 0.8 eV (blue), the region shifting from left to right.

Figure 12. Behavior of sticking probabilities versus both θCO and θO for our refined model. Plotted are contours for ln S at increments of −0.5 decreasing from the lower left where S = 1 (and ln S = 0) to the upper right where S → 0 (and ln S → −∞) for θCO + θO → 1/2.

Figure 13a compares KMC simulation results with the analytic theory for the choice Erxn = 1.0 eV. There is reasonable agreement in the value of the critical temperature, Tc. Significantly, neither the experimental data nor the modeling display the expected and/or often assumed V-shaped bistable region9,57 in these standard log pCO or log PCO versus 1000/T plots. The apex of such a V corresponds to the critical point terminating the bistable region on the right (for higher T), and the V opens to the left (lower T). Rather, after opening near the critical point, we find that the bistable region closes (i.e., becomes narrower as pCO+ − pCO− decreases) for lower T. Our analytic rate equation formulation combined with a bifurcation analysis also allows ready assessment of other aspects of behavior. We have already indicated in section 4.1 that the critical temperature, Tc, is largely independent of the choice of reaction barrier, Erxn, but rather controlled by Edes (prompting extensive analysis of EaCO = −Edes in section 3). Indeed, Figure 12b reveals that the critical point is effectively invariant as Erxn varies between 0.8 and 1.1 eV (although the shape of the bistability region changes significantly). We have performed additional analysis of the dependence of Tc on Edes selecting Erxn = 1.00 eV. This analysis reveals that Tc = 488.2, 520.2, 553.5,

SCO from (8) and SO2 from (16) and (24) as a function of both θCO and θO. Here, CO adsorption is quite facile even for higher θCO. Specifically, one has that SO2 ∼ (8/3)(1/2 − θCO)2 from (16) and (24) is well below SCO ∼ (512/19)(1/2 − θCO)2 from (8) for higher θCO close to θCO(max) = 1/2 ML (where θO ≈ 0), which from section 2 produces a near CO-poisoned state. In addition, SO2 ∼ 215(1/2 − θO)8 is far below SCO ∼ 64(1/2 − θO) for higher θO (and θCO ≈ 0) consistent with the requirement from section 2 for existence of a conventional reactive steady state. We have performed a more extensive analysis of model behavior mapping out the complete bifurcation diagram both using KMC simulations, and also based on our analytic beyondMF rate equations utilizing AUTO continuation and bifurcation software.14 The results are shown in Figure 13. Experimental data of Vogel et al.57 for the (CO+O)/Pd(100) system are also shown as an inset in Figure 13a. The latter were obtained using local photoemission electron microscopy (PEEM) of (100) facets to assess the state of the system, specifically cycling PCO for fixed PO2 = 1.3 × 10−5 mbar to assess hysteresis in the PEEM intensity associated with the presence of bistability. The main frame in J


Article

The Journal of Physical Chemistry C 586.2, and 618.8 K for Edes = 1.5, 1.6, 1.7, 1.8, and 1.9 eV, respectively, so that Tc/Edes ≈ 325.5(±0.3) K/eV is roughly constant over this range. This latter relation is useful allowing estimation of Edes from experimental observation of Tc. 6.3. Discussion. We have emphasized in section 1 that traditional MF rate equation treatments fail to quantitatively describe behavior observed in models with multiple adsorption sites and strong adlayer ordering (e.g., induced by short-range exclusion). This feature is explicitly illustrated here using as an example the model of Hoffman et al.10,11 for CO oxidation on Pd(100). First, we consider the classic MF site approximation where all spatial correlations are ignored and thus multisite probabilities are factorized into a product of constituent single site probabilities. In this approximation, with a Langmuirian treatment of CO-adsorption and an 8-site model for O2 adsorption, one has SCO(site) = P9b + 2h ≈ (1 − θCO br)9 (1 − θO)2 br 8

and 8

SO2 (site) = P8h + 8b ≈ (1 − θCO ) (1 − θO)

θCObr

Figure 14. Comparison of MF site-approximation (solid curves) and pair-approximation (dashed curves) predictions for steady-state values of (a) θCO, (b) θO, and TOF (inset), for the model of Hoffman et al.10,11 with our precise KMC results (dots).

Finally, we discuss the failure of our refined factorization approximations (and also the Padé-resummation formulation) utilized in our beyond-MF rate equation analysis to capture the non-MF bistability or multistability observed for θCO very close to θCO(max) = 1/2 ML. The reason for failure is simply that these treatments do not account for order−disorder symmetry breaking in the CO adlayer in this regime. Such symmetrybreaking transitions can occur for adpsecies with short-range exclusion interactions when the coverage approaches the saturation value.45 One expects somewhat different forms of sticking probabilities induced by such symmetry breaking. For example, for strongly symmetry-broken CO-rich states for θCO very close to θCO(max) = 1/2 ML with c(2×2) ordering, a simple model for adlayer statistics is to assume that CO adspecies are randomly removed or missing from otherwise perfect c(2×2) ordering. This model implies that SCO ∼ 1 − 2θCO decays linearly as θCO → 1/2 ML (versus faster decay in the Langmuirian and steering + funneling formulations). It is this special behavior that leads to the non-MF type bistability or multistability described in sections 6.1 and 6.2. It is also appropriate to mention that this type of subtle non-MF type behavior is no doubt difficult to observe experimentally, and indeed it has not yet been observed in the (CO + O)/Pd(100) system. However, it has been recently identified experimentally in the similar (CO + O)/Pt(111) system.58 It should also been noted that our formulation does not account for symmetry breaking in c(2×2)-O adlayers.59 However, this feature does not qualitatively impact steady-state behavior.

(18)

Recalling that = /2θCO, this treatment produces θCO(max) = 2 ML (4 times the correct value) and θO(max) = 1 ML (2 times the correct value). Because SCO(site) ≪ SO2(site) for θCO approaching θCO(max) (with θO ≈ 0), we do not expect MFbistability due to the lack of a near CO-poisoned state. Employing instead a Kirkwood-type pair approximation yields 1

SCO(pair) = P9bQ 2h | 9b

with

br 12

P9b ≈ (1 − 2θCO ) /(1 − θCO br)15 and Q 2h | 9b ≈ 1 − 2qO(pair)

(19)

with the conditional coverage for O surrounded by 4 NN empty br sites now given by qO(pair) ≈ θO(1 − 2θCObr)4/(1 − θCObr)4. See the Appendix. One also obtains SO2(pair) = P8hQ 8b | 8h

where

Q 8b | 8h(pair) ≈ (Q 2b)10 /(Q 1b)12

(20)

Qnb = 1 − nqCObr with the CO qCObr = θCObr/(1 − 2θO) as in

where P8h is given by (9), and conditional coverage given by section 5.3. In these pair factorizations, we account for NN pairs of br sites for CO, and NN 4fh sites for O, in factorizing multisite motifs. An alternative also accounts for 2NN pairs of br sites.50 See the Appendix. Thus, from (19) and (20), the pair approximation produces θCO(max) = 1 ML (twice the correct value) and θO(max) = 1/2 ML (the correct value). Because SCO(pair) ≪ SO2(pair) for θCO approaching θCO(max) (with θO ≈ 0), we do not expect MF-type bistability. Predictions for steady-state behavior from the site- and pair-approximations are shown in Figure 14 and compared against our precise KMC results, which are accurately recovered by our more reliable higher-order beyond-MF rate equation treatment as shown in Figure 9. As expected, predictions for θCO versus pCO are poor (but those for θO and TOF are better). From the above, all treatments of the model of Hoffmann et al.10,11 consistently recover the lack of MF-bistability. However, we should note that various modifications to incorporate more flexible CO adsorption do produce MF-type bistability. See the Supporting Information. Thus, it is not just our refined model of section 5.2 with CO steering and funneling that recovers the experimentally observed behavior.

7. CONCLUSIONS Traditional MF rate equations for the kinetics of catalytic surface reactions are based on a picture involving a single type of adsorption site and neglecting spatial correlations in the adlayer. However, for most reactions on crystalline metal surfaces, different reactants often occupy distinct types of adsorption sites, and sometimes a single reactant occupies multiple site types. This is the case for CO oxidation on Pd(100), where O exclusively occupies hollow sites, but CO prefers to occupy br sites. Furthermore, MF rate equations neglect the ubiquitous superlattice adlayer ordering due to strong short-range adspecies repulsions. For various tailored msLG reaction models incorporating these features, we are able to develop analytic beyond-MF rate equations that reliably recover key aspects of KMC simulation results for these models in the regime of high surface mobility where KMC simulation is computationally demanding. Our approach extends Kirkwood-type approximaK


Article

The Journal of Physical Chemistry C tion ideas and also exploits Padé-type approximants of exact lowcoverage Taylor expansions for multisite probabilities which describe sticking on mixed adlayers. This approach not only provides deeper insight into behavior than KMC simulation, e.g., revealing the importance of facile CO adsorption for MF bistability, but also allows efficient mapping of bifurcation behavior versus experimental control parameters using analytic continuation software. Our beyond-MF rate equation approach applies beyond CO oxidation to thus significantly extend the utility of rate equation formulations for surface reactions. Next, we note that there exists more detailed modeling for (CO + O)/Pd(100).6,26 This modeling includes longer-range weak adspecies interactions stabilizing subtle c(2×2√2)R45° ordering and also p(2×2) ordering for O, although this does not significantly affect reaction kinetics at least for higher T approaching Tc.6,26 Also included is CO population of both br and 4fh adsorption sites. This detailed modeling incorporates dissociative oxygen adsorption via a NN pair of 4fh sites with barrier Eact = Eads + ΔETS.6 Here, Eads is the adsorption energy incorporating lateral adspecies interactions at these NN 4fh sites, and ΔETS is an additional barrier. Direct adsorption primarily at a 2NN vicinal br pairs of sites is quickly followed by thermally activated hopping separating constituent O adatoms typically to third NN 4fh sites. This two-stage process is simplified to a onestep process in our tailored modeling. On the basis of new KMC simulations, we suggest an optimum choice of ΔETS ≈ 0.75 eV (versus ΔETS ≈ 1.0 eV used previously54), which yields a realistic saturation O coverage between 0.3 and 0.4 ML at 300 K46,60 corresponding to that in our tailored modeling.50 Finally, it is natural to consider extension of our analytic beyond-MF approach to more complex and realistic models. We consider three related issues. First, the treatment of multiple adsorption sites for CO is straightforward, although now more complex multisite configuration probabilities describe sticking probabilities and need to be factorized. Second, incorporation of the effect of finite longer-range adspecies interactions is also possible at least with additional MF or Kirkwood type approximation. Such CO−CO interactions are quite weak, so instead consider as an example the effect of a significant 2NN O− O and d = √5a/2 br-4fh CO−O replusions.6,26 Retaining the current prescription of empty site ensembles for sticking, the associated probabilities will be effected by these interactions should they include such site pairs. One could implement a refined Kirkwood approximation including such pairs, and use Virial expansion to represent associated pair probabilities. Detailed-balance means that CO desorption will now have rates dependent on the presence of nearby O separated by d = √5a/2. Prescription of reaction is also impacted by the CO−O repulsion. Third, another challenge is development of analytic descriptions for high-coverage regimes where there exist symmetry-breaking transitions in adlayer ordering. However, this is possible by incorporating an appropriate order parameter into factorization approximations.61

Figure 15. Configuration and conditional probabilities. Solid red (blue) circles denote CO (O) on br (4fh) sites. Open red (blue) circles denote unoccupied sites. Dashed open circles denote sites that are specified to be unoccupied. (a) Compact configurations of empty br sites. (b) Conditional coverage, qCObr, and conditional probabilities, Qnb ≈ 1 − nqCObr, for compact motifs of n empty br sites. (c) Empty NN pair of 4fh sites. (d) Conditional coverage, qO, and conditional probability, Q2h ≈ 1 − 2qO, for an empty NN 4fh pair.

ment of the latter quantities is presented here. Their application to the derivation of sticking probability expressions is described with additional results presented for O2 adsorption models with reorientation. As in section 5, let qCObr denote the conditional coverage for CO to occupy a br site where both NN 4fh sites separated by d = a/2 are specified to be unoccupied by O. See Figure 15b. Then, qCObr equals the probability of the configuration with a CO between two unoccupied 4fh sites (which just equals θCObr) divided by the probability of an empty NN 4fh pair (P2h = 1 − 2θO). Thus, one has the exact relation qCObr = θCObr/(1 − 2θO). In our analysis of sticking, we will need to assess the conditional probabilities, Qnb, of various compact motifs of n unoccupied br sites of which at most one can be occupied by CO, such that each br site is between a specified empty pair of 4fh sites. See Figure 15b. We will apply the reasonable approximation that Qnb ≈ 1 − nqCObr for n > 1 with Q2b ≈ Q2b′. To illustrate the nature of the approximation, note that using conservation of probability, one can exactly write Q2b = 1 − 2q*CObr where q*CObr is the conditional coverage for CO on a br site given that the two NN 4fh sites are unoccupied as well as one 4fh site separated by d = √5a/2. Thus, we approximate q*CObr by the similar qCObr. Note also that Q1b = 1 − qCObr is exact. In section 5.3 (section 5.4), we made an additional approximation factorizing Q8b|8b (Q8b|9h) in terms of the Qnb. Also as in section 5, let qO denote the conditional probability for O on a 4fh site where all 4 of the neighboring br sites separated by d = a/2 are specified to be unoccupied by CO. See Figure 15d. Then, qO equals the probability of the configuration with an O with 4 NN unoccupied br sites (which just equals θO) divided by the probability of a diamond-shaped quartet of 4 empty br sites (P4b). Our high-level beyond MF treatment used the approximation P4b = 1 − 4θCObr = 1 − 2θCO. The lower-level Kirkwood type pair approximation in section 6.3 used the approximation P4b ≈ (P2b)4/(P1b)4 where P2b = 1 − 2θCObr and

■

APPENDIX: CONDITIONAL PROBABILITIES AND DERIVATION OF STICKING PROBABILITIES In the introduction to section 5, we have introduced reactant adspecies coverages θCO = 2θCObr for CO and θO for O and presented exact relations for probabilities, Pnb (Pnh), of empty br (4fh) site motifs at most one of which can be populated by CO (O). See Figure 15a and Figure 15c. We have also indicated in section 5 the utility of introducing various conditional coverages and related conditional probabilities. A more detailed developL


Article

The Journal of Physical Chemistry C P1b = 1 − θCObr. Next, consider the conditional probability, Q2h, of a NN pair of unoccupied 4fh sites each of which is surrounded by a specified empty quadruple of br sites. See Figure 15d. Analogous to the discussion of Q2b above, one can use conservation of probability to exactly write Q2h = 1 − 2q*O where q*O is the conditional coverage for O on a 4fh site given that the 7 br sites shown in Figure 15d are all empty. Then, we make the reasonable approximation q*O ≈ qO. In addition, in section 5.1, we use the additional reasonable approximation that Q2h|9b ≈ Q2h. As indicated in section 5.5, for oxygen adsorption models where one of two adsorption orientations (1 or 2) are possible, the sticking probability has the form, SO2 = P(1) + P(2) − P(1∧2), where P(1) = P(2) = P8h+8b (P9h+8b) and P(1∧2) = P12h+12b (P13h+16b) for the 8-site (9-site) model. We have analyzed P8h+8b (P9h+8b) in section 5.3 (section 5.4). The configurations associated with P12h+12b (P13h+16b) that have 12 empty 4fh sites and 12 empty br sites (13 empty 4fh sites and 16 empty br sites) are shown in Figure 16. We now describe the appropriate factorization of these latter quantities.

P13h + 16b ≈ (1 − 2θO)15 (1 − 2θCO − 2θO)9 × (1 − 1 2 θCO − 2θO)4 (1 − θO)−19

■

× (1 − θCO − 2θO)−12

(24)

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.6b10102. Cluster geometries for analysis of CO adsorption at br sites (S1); additional plane-wave DFT VASP analysis (S2); results for various levels of theory with localized atomic bases for CO+8Pd (S3); additional results with unrestricted spin for PBE and PBE0 (S4); modifications to the model of Hoffmann et al. producing MF-type bistability (S5) (PDF)

■

AUTHOR INFORMATION

Corresponding Authors

*D.-J. Liu. E-mail: [email protected]. *J. W. Evans. E-mail: [email protected]. ORCID

Da-Jiang Liu: 0000-0002-3019-9247 Federico Zahariev: 0000-0002-3223-3576 James W. Evans: 0000-0002-5806-3720 Notes

Figure 16. Configurations for the 8-site and 9-site models for O2 adsorption with reorientation that allow adsorption in either of the two possible orthogonal orientations.

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank M.W. Schmidt for useful suggestions regarding quantum chemistry analysis. This work was supported by the U.S. Department of Energy (USDOE), Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences, and Biosciences through the Ames Laboratory Chemical Physics program. We acknowledge the use of NERSC computational resources for analysis of CO adsorption energetics. The work was performed at Ames Laboratory which is operated for the USDOE by Iowa State University under Contract No. DE-AC0207CH11358.

For the 8-site model with reorientation, P12h+12b (Figure 16a) is determined from P12h + 12b = P12hQ 12b | 12h

with

P12h = (P2h)16 /(P1h)20

and Q 12b | 12h = (Q 4b)5 (Q 3b)4 (Q 1b)4 /(Q 2b)12

(21)

The factorization of P12h uses that the ensemble of 12 4fh sites involves 16 NN pairs that overlap and share many single 4fh sites. The factorization of Q12b|12h uses that the ensemble of 12 br sites can be regarded as consisting of 5 quartets and 4 triangles that overlap and share many NN pairs of br sites. From (21), one finds that

■

P13h + 16b ≈ (1 − 2θO)15 (1 − 2θCO − 2θO)5 × (1 − 3 2 θCO − 2θO)4 (1 − 1 2 θCO − 2θO)4 × (1 − θO)−20 (1 − θCO − 2θO)−12

(22)

For the 9-site model with reorientation, P13h+16b (Figure 16b) is determined from P13h + 16b = P13hQ 16b | 13h

with

P13h = (P2h)16 /(P1h)19

and Q 16b | 13h = (Q 1b)4 (Q 4b)9 /(Q 2b)12

REFERENCES

(1) Norskov, J. K.; Abild-Pederesen, F.; Studt, F.; Bligaard, T. Density functional theory in surface chemistry and catalysis. Proc. Natl. Acad. Sci. U. S. A. 2011, 108, 937−943. (2) Imbihl, R.; Ertl, G. Oscillatory Kinetics in Heterogenous Catalysis. Chem. Rev. 1995, 95, 697−733. (3) Imbihl, R. Non-linear Dynamics in Catalytic Reactions. In Handbook in Surface Science; Hasselbrink, E., Lundqvist, B. I., Eds.; Elsevier: Amsterdam, 2008; Vol. 3, Chapter 9. (4) Reuter, K. In Modelling Heterogeneous Catalytic Reactions: From the molecular process to the technical system; Deutschmann, O., Ed.; WileyVCH: Weinheim, 2011. (5) Stamatakis, M.; Vlachos, D. G. Unraveling complexity of catalytic reactions via kinetic Monte Carlo simulation: Current status and frontiers. ACS Catal. 2012, 2, 2648−2663. (6) Liu, D.-J.; Evans, J. W. Realistic Multisite Lattice-gas Modeling and KMC Simulation of Catalytic Surface Reactions: Kinetics and Multiscale Spatial Behavior for CO-oxidation on Metal(100) Surfaces. Prog. Surf. Sci. 2013, 88, 393−521. (7) van Santen, R. A.; Niemantsverdriet, J. W. Chemical Kinetics and Catalysis; Plenum: New York, 1995.

(23)

The factorization of P13h uses the feature that the ensemble of 13 4fh sites involves 16 NN pairs that overlap and share many single 4fh sites. The factorization of Q16b|13h uses the feature that the ensemble of 16 br sites can be regarded as consisting of 9 quartets that overlap and share many pairs of NN br sites. One finds that M


Article

The Journal of Physical Chemistry C (8) Liu, D.-J.; Evans, J. W. Fronts and fluctuations in a tailored model for CO oxidation on metal(100) surfaces. J. Phys.: Condens. Matter 2007, 19, 065129. (9) Liu, D.-J.; Garcia, G.; Wang, J.; Ackerman, A. M.; Wang, C.-J.; Evans, J. W. Kinetic Monte Carlo simulation of statistical mechanical models and coarse-grained mesoscale descriptions of catalytic reactiondiffusion processes: 1D nanoporous and 2D surface systems. Chem. Rev. 2015, 115, 5979−6050. (10) Hoffmann, M. J.; Reuter, K. CO Oxidation on Pd(100) Versus PdO(101)-(5√ × 5√)R27°(5 × 5)R27°: First-Principles Kinetic Phase Diagrams and Bistability Conditions. Top. Catal. 2014, 57, 159− 170. (11) 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. (12) Reuter, K. Ab Initio Thermodynamics and First-Principles Microkinetics for Surface Catalysis. Catal. Lett. 2016, 146, 541−563. (13) Baxter, R. J. Exactly Solved Models in Statistical Mechanics; Academic Press: London, 1982. (14) AUTO software: http://www.macs.hw.ac.uk/~gabriel/auto07/ auto.html. (15) Bar, M.; Zulicke, Ch.; Eiswirth, M.; Ertl, G. Theoretical modeling of spatiotemporal self-organization in a surface catalyzed reaction exhibiting bistable kinetics. J. Chem. Phys. 1992, 96, 8595−8604. (16) Liu, D.-J; Pavlenko, N.; Evans, J. W. Crossover between MeanField and Ising Critical Behavior in a Lattice-Gas Reaction-Diffusion Model. J. Stat. Phys. 2004, 114, 101−114. (17) Tammaro, M.; Sabella, M.; Evans, J. W. Hybrid Treatment of Spatio-Temporal Behavior in Surface Reactions with Coexisting Immobile and Highly Mobile Reactants. J. Chem. Phys. 1995, 103, 10277−10285. (18) Evans, J. W.; Liu, D.-J.; Tammaro, M. From Atomistic Lattice-Gas Models for Surface Reactions to Hydrodynamic Reaction-Diffusion Equations. Chaos 2002, 12, 131−143. (19) Kresse, G.; Furthmüller, J. Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. Comput. Mater. Sci. 1996, 6, 15−50. (20) Blöchl, P. E. Projector augmented-wave method. Phys. Rev. B: Condens. Matter Mater. Phys. 1994, 50, 17953−17979. (21) Kresse, G.; Joubert, D. From ultrasoft pseudopotentials to the projector augmented-wave method. Phys. Rev. B: Condens. Matter Mater. Phys. 1999, 59, 1758−1775. (22) Methfessel, M.; Paxton, A. T. High-precision sampling for Brillouin-zone integration in metals. Phys. Rev. B: Condens. Matter Mater. Phys. 1989, 40, 3616−3621. (23) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1997, 78, 1396. (24) Liu, D.-J. Density functional analysis of key energetics in metal homoepitaxy: quantum size effects in periodic slab calculations. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 81, 035415. (25) Mason, S. E.; Grinberg, I.; Rappe, A. M. First-principles extrapolation method for accurate CO adsorption energies on metal surfaces. Phys. Rev. B: Condens. Matter Mater. Phys. 2004, 69, 161401. (26) Liu, D.-J; Evans, J. W. Atomistic lattice-gas modeling of COoxidation on Pd(100): Temperature-programmed spectroscopy and steady-state behavior. J. Chem. Phys. 2006, 124, 154705. (27) Heyd, J.; Scuseria, G. E.; Ernzerhof, M. Hybrid functionals based on a screened Coulomb potential. J. Chem. Phys. 2003, 118, 8207−8215. (28) Valiev, M.; Bylaska, E.; Govind, N.; Kowalski, K.; Straatsma, T.; van Dam, H.; Wang, D.; Nieplocha, J.; Apra, E.; Windus, T.; et al. NWChem: a comprehensive and scalable open-source solution for large scale molecular simulations. Comput. Phys. Commun. 2010, 181, 1477− 1489. (29) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; et al. J. Comput. Chem. 1993, 14, 1347−1363. (30) Gordon, M. S.; Schmidt, M. W. In Theory and Applications of Computational Chemistry, the first forty years; Dykstra, C. E., Frenking, G.,

Kim, K. S., Scuseria, G. E., Eds.; Elsevier: Amsterdam, 2005; pp 1167− 1189. (31) Hay, P. J.; Wadt, W. R. Ab initio effective core potentials for molecular calculations. Potentials for the transition metals atoms Sc to Hg. J. Chem. Phys. 1985, 82, 270−283. (32) Boschen, J. S.; Lee, J.; Windus, T. L.; Evans, J. W.; Liu, D.-J. Size Dependence of S-bonding on (111) Facets of Cu Nanoclusters. J. Phys. Chem. C 2016, 120, 10268−10274. (33) Clark, T.; Chandrasekhar, J.; Spitznagel, G. W.; von Ragué Schleyer, P. Efficient diffuse function-augmented basis sets for anion calculations. III. The 3-21+G basis set for first-row elements, Li-F. J. Comput. Chem. 1983, 4, 294−301. (34) Adamo, C.; Barone, V. Toward reliable density functional methods without adjustable parameters: The PBE0 model. J. Chem. Phys. 1999, 110, 6158−6170. (35) Oba, Y.; Okamoto, H.; Sato, T.; Shinohara, T.; Suzuki, J.; Nakamura, T.; Muro, T.; Osawa, H. X-ray magnetic circular dichroism study on ferromagnetic Pd nanoparticles. J. Phys. D: Appl. Phys. 2008, 41, 134024. (36) Moseler, M.; Häkkinen, H.; Barnett, R. N.; Landman, U. Structure and Magnetism of Neutral and Anionic Palladium Clusters. Phys. Rev. Lett. 2001, 86, 2545−2548. (37) Sakuragi, S.; Sakai, T.; Urata, S.; Aihara, S.; Shinto, A.; Kageshima, H.; Sawada, M.; Namatame, H.; Taniguchi, M.; Sato, T. Thicknessdependent appearance of ferromagnetism in Pd(100) ultrathin films. Phys. Rev. B: Condens. Matter Mater. Phys. 2014, 90, 054411. (38) Zhang, C. J.; Hu, P. CO Oxidation on Pd(100) and Pd(111): A Comparative study of reaction pathways and reactivity for low and medium coverage. J. Am. Chem. Soc. 2001, 123, 1166−1172. (39) Hammer, B. The NO+CO Reaction Catalyzed by Flat, Stepped, and Edged Pd Surfaces. J. Catal. 2001, 199, 171−176. (40) Eichler, A. CO oxidation on transition metal surfaces: reaction rates from first principles. Surf. Sci. 2002, 498, 314−320. (41) Chang, S.-L.; Thiel, P. A. Formation of a metastable ordered surface phase due to competitive diffusion and adsorption kinetics: Oxygen on Pd(100). Phys. Rev. Lett. 1987, 59, 296−299. (42) Liu, D.-J.; Evans, J. W. Lattice-Gas modeling of the formation and ordering of oxygen adlayers on Pd(100). Surf. Sci. 2004, 563, 13−26. (43) Sun, J.; Marsman, M.; Csonka, G. I.; Ruzsinszky, A.; Hao, P.; Kim, Y.-S.; Kresse, G.; Perdew, J. P. Self-consistent meta-generalized gradient approximation within the projector-augmented-wave method. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 84, 035117. (44) Zhang, Y.; Blum, V.; Reuter, K. Accuracy of first-principles lateral interactions: Oxygen on Pd(100). Phys. Rev. B: Condens. Matter Mater. Phys. 2007, 75, 235406. (45) Liu, D.-J. Lattice-gas modeling of CO adlayers on Pd(100). J. Chem. Phys. 2004, 121, 4352−4357. (46) Chang, S.-L.; Thiel, P. A. Oxygen on Pd(100): Order, Reconstruction, and Desorption. J. Chem. Phys. 1988, 88, 2071−2082. (47) Behm, R. J.; Christmann, K.; Ertl, G.; Van Hove, M. A. Adsorption of CO on Pd(100). J. Chem. Phys. 1980, 73, 2984−2996. (48) Payne, S. H.; McEwen, J. S.; Kreuzer, H. J.; Menzel, D. Adsorption and desorption of CO on Ru(0001): a comprehensive analysis. Surf. Sci. 2005, 594, 240−262. (49) Eichler, A.; Hafner. Adsorption of CO on Pd(100): Steering into less favored adsorption sites. Phys. Rev. B: Condens. Matter Mater. Phys. 1998, 57, 10110−10114. (50) Evans, J. W.; Liu, D.-J. Statistical mechanical models for dissociative adsorption of O2 on metal(100) surfaces with blocking, steering, and funneling. J. Chem. Phys. 2014, 140, 194704. (51) Brundle, C. R.; Behm, R. J.; Barker, J. A. Summary Abstract: How many metal atoms are involved in dissociative chemisorption? J. Vac. Sci. Technol., A 1984, 2, 1038−1039. (52) Zangwill, A. Physics at Surfaces; Cambridge University Press: Cambridge, U.K., 1988. (53) Evans, J. W. Non-Equilibrium Percolative c(2 × 2) Ordering: Oxygen on Pd(100). J. Chem. Phys. 1987, 87, 3038−3048. N


Article

The Journal of Physical Chemistry C (54) Liu, D.-J.; Evans, J. W. Dissociative adsorption of O2 on unreconstructed metal(100) surfaces: Pathways, energetics, and sticking kinetics. Phys. Rev. B: Condens. Matter Mater. Phys. 2014, 89, 205406. (55) Meyer, J.; Reuter, K. Modeling Heat Dissipation at the Nanoscale: An Embedding Approach for Chemical Reaction Dynamics on Metal Surfaces. Angew. Chem., Int. Ed. 2014, 53, 4721−4724. (56) Evans, J. W. Comment on “Kinetics of Random Sequential Adsorption. Phys. Rev. Lett. 1989, 62, 2642. (57) Vogel, D.; Spiel, C.; Suchorski, Y.; Trinchero, A.; Schloegl, R.; Groenbeck, H.; Rupprechter, G. Local catalytic ignition during CO oxidation on low-index Pt and Pd surfaces: a combined PEEM, MS, and DFT study. Angew. Chem., Int. Ed. 2012, 51, 10041−10046. (58) Wrobel, R. J.; Becker, S.; Weiss, H. J. Second/additional Bistability in a CO Oxidation Reaction on Pt(111): An Extension and Compilation. J. Phys. Chem. C 2012, 116, 22287−22292. (59) Liu, D.-J.; Evans, J. W. Symmetry-Breaking and Percolation Transitions in a Surface Reaction Model with Superlattice Ordering. Phys. Rev. Lett. 2000, 84, 955−958. (60) den Dunnen, A.; Jacobse, L.; Wiegman, S.; Berg, O. T.; Juurlink, L. B. F. Coverage-dependent adsorption and desorption of oxygen on Pd(100). J. Chem. Phys. 2016, 144, 244706. (61) Fowler, R.; Guggenheim, E. A. Statistical Thermodynamics; Cambridge University Press: Cambridge, U.K., 1965; Chapter XIII, p 1322.

O


Predictive Beyond-Mean-Field Rate Equations for ... - ACS Publications

Recommend Documents