We explore simple lattice-gas reaction models for CO-oxidation on 1D and 2D periodic arrays of surface adsorption sites with CO adsorption and desorption, dissociative O2 adsorption and recombinative desorption (at low rate), and CO + O reaction to form CO2. Adspecies interactions are neglected, and adspecies diffusion is effectively absent. The models are motivated by studies of CO-oxidation on RuO2(110) at high-pressures. Despite the lack of adspecies interactions, negligible adspecies diffusion results in kinetically induced spatial correlations. A transition occurs from a random primarily CO-populated steady-state at high CO-partial pressure, pCO, to a strongly correlated near-O-covered steady-state for low pCO as noted by Matera et al. [J. Chem. Phys. 134, 064713 (2011)]. In addition, we identify a second transition to a random near-O-covered steady-state at very low pCO. Furthermore, we identify and analyze the slow “diffusive dynamics” for very low pCO and provide a detailed characterization of the crossover to the strongly correlated O-covered steady-state as well as of the spatial correlations in that state.

Traditional mean-field (MF) rate equations of chemical kinetics1 have proven immensely useful in elucidating catalytic surface reaction processes. However, it is well-recognized standard MF rate expressions neglect spatial correlations in the reactant and/or product distribution. For example, traditional kinetics ignores lateral correlations in mixed chemisorbed reactant adlayers in surface reactions such as catalytic CO-oxidation,2,3 the process of interest here. Generically, there are two sources of such correlations. They can derive from either (i) thermodynamics where the reactant and product distribution are locally equilibrated, and the associated Gibbs distribution reflects interactions between molecules or (ii) kinetics for unequilibrated distributions where correlations induced by the reaction kinetics are not erased due to low molecular mobility. The former situation generally applies for surface reactions under lower-pressure (P) conditions with low surface coverages and high adspecies surface diffusivity.3 In contrast, for high-P conditions,2 adspecies diffusivity is strongly inhibited at near-saturation or jamming coverages, and kinetically induced correlations can reflect the details of adsorption, desorption, and reaction processes.

Spatial correlations in equilibrated adlayers can be determined precisely by Monte Carlo simulation, finite-size scaling, and transfer matrix techniques. If such correlations are not too strong, they can be effectively characterized via Ursell-Mayer cluster expansions, or Kirkwood or quasi-chemical approximations.4,5 These correlations also satisfy a spatial Markov property.6 Realistic modeling of mixed reactant adlayer structure requires accurate determination of the site-specific adsorption energies (noting that, in general, adspecies can occupy multiple distinct site types) and of the numerous adspecies interactions between the same and distinct adspecies.3 For low-P CO-oxidation, this type of realistic modeling was first implemented for metal(100) surfaces,3,7–9 and more recently for metal(111) surfaces.10,11

However, our focus in this work is instead on elucidating the kinetically induced spatial correlations in surface reactions at high-P, and specifically in CO-oxidation on oxide surfaces.12–16 Such correlations are less well characterized than in the equilibrium case. Unfortunately, non-equilibrium analogues of cluster expansions and concepts related to Markovian spatial statistics do not generally apply. However, for reactive non-equilibrium steady states (ss) in an open system (reactants in, products out) such as CO-oxidation in a flow reactor, one can make some general observations about their behavior. Such a bimolecular reaction will lead to depletion of the population of associated nearby reactant pairs (i.e., CO—O pairs), and boost the population of nearby vacancy pairs relative to a corresponding equilibrium state. Note again that these kinetically induced correlations would be erased by sufficiently high surface mobility.

Also, at least for simpler surface reaction processes with adspecies at well-defined adsorption sites, it is straightforward to develop an associated set of exact hierarchical master equations for the evolution of the probabilities of various configurations of populated or empty sites.17–22 The lowest-order equations in this hierarchy describe the evolution of standard (single-site) adspecies concentrations or coverages. However, these equations couple to the probabilities of multi-site configurations, e.g., for nearby reactant pairs (which control the reaction rates), and for ensembles of empty sites which might be required for dissociative adsorption. Probabilities of even larger configurations are involved if the rates for adsorption, reaction, etc., depend on the local environment. It is thus necessary to develop equations for the probabilities of pairs and any larger configurations, the evolution of which couples to probabilities of even larger configurations, thus generating an infinite hierarchy.

The simplest site-approximation treatment of this hierarchy factorizes the probabilities of larger configurations as a product of single-site quantities, thus neglecting all spatial correlations and recovering MF kinetics. An improved Kirkwood-type or pair approximation (PA) retains the probabilities of adjacent pairs and factorizes those of larger configurations in terms of pair and single-site quantities.17–22 In a triplet approximation, one retains probabilities of triples and factorizes those of larger configurations in terms of these, pairs, and single-site probabilities. These and higher-order approximations do provide insight into kinetically induced (as well as thermodynamic) correlations. This approach has been implemented for various simple surface reaction models with limited or no mobility.18–22 We caution that there are limitations. For example, any finite-order truncation approximation will produce asymptotically exponential kinetics of approach to the steady-state, as the corresponding finite-dimensional evolution operator exhibits a spectral gap. However, the actual dynamics described by the infinite hierarchy might correspond to slower gapless diffusive asymptotic evolution. We provide an example of this behavior below.

The simple surface reaction model considered in this work is motivated by a tailored model of Reuter and Scheffler12,13 which was crafted to capture the essential features of CO-oxidation on the RuO2(110) oxide surface under high-P conditions. The RuO2(110) surface consists of a rectangular grid of adsorption sites which consist of alternating rows of so-called coordinatively unsaturated (cus) sites and bridge (br) sites.(See Appendix  A). Processes involved in the reaction are adsorption and desorption of CO at single sites of either type; dissociative adsorption O2 and recombinative desorption at nearest-neighbor (NN) pairs of sites of any type; irreversible reaction of CO and O on NN pairs of sites of any type to form the product CO2 which immediately desorbs from the surface; and very limited surface diffusion. Detailed density functional theory (DFT) analysis was employed to determine the barriers for all these desorption and reaction processes, the rates adopting an Arrhenius form. Adsorption rates were determined from CO and O2 partial pressures. The modeling neglected adspecies lateral interactions and also suggested that diffusion had a negligible effect under these high-P high-coverage conditions.12,13 Additional work by these authors together with Metiu applied a “degree of rate control” approach to show that the dominant processes occurred on cus sites or adjacent cus pairs, which effectively reduce the model on a 2D rectangular grid of sites to one on a 1D linear lattice of sites.23 As an aside, separate studies by Over and coworkers suggested a modified set of barriers and rates for this system which could significantly impact behavior.16 

However, in the current study, we instead focus on elucidating distinct generic regimes of behavior and the associated spatial correlations in the mixed reactant adlayer. Of specific interest is the observation by Reuter and coworkers24,25 of a transition from primarily CO-populated surface for higher CO-partial pressure, pCO, with a random distribution of CO (or of isolated vacancies) to a near O-covered surface for lower pCO but with strong correlations in the distribution of O adatoms (or of vacancies). These correlations invalidate a treatment based on MF kinetics. Our study of a simple model mimicking the high-P (CO + O)/RuO2(110) system will focus on transitions between correlated and uncorrelated steady-state regimes, in particular identifying a second transition to a random near oxygen-covered surface for very low pCO.

In Sec. II, we describe our reaction model and indicate its steady-state behavior. In Sec. III, we provide a detailed analysis of the pure oxygen adsorption-desorption model for pCO = 0, as its spatial and dynamical features impact the reaction model for low pCO. In Sec. IV, we provide a detailed analysis of steady-state behavior of the reaction model in distinct regimes, including an analysis of crossover behavior and a detailed elucidation of the nature of strong correlations in the near-O-covered surface low pCO. Conclusions are provided in Sec. V.

We consider a simple reaction model for both a 1D linear (d = 1) and 2D square (d = 2) grid of adsorption sites which incorporates the same processes as in the model of Reuter and Scheffler:12,13 reversible CO adsorption on single sites, dissociative O2 adsorption on adjacent sites and the reverse recombinative desorption process, and irreversible reaction of CO and O on adjacent sites. Our rates are chosen to reflect the relative magnitudes of those in the (CO + O)/RuO2(110) system. The adsorption rates are chosen so that pCO + pO2 = 1 setting the time scale. Here, pCO is the adsorption rate of CO per site per unit time. Also, pO2 (pO2/2) is the dissociative adsorption rate and dO2 (dO2/2) is the recombinative desorption rate for oxygen per pair of NN sites for d = 1 (d = 2). Note that the 2D system has two NN pairs per site. A key parameter in the model is ε = dO2/pO2 and we will focus in our simulation study on the range from 10−6 (a typical value in experiment) up to 10−2. The desorption rate of CO is set to dCO = 0.1, and the reaction rate per NN CO—O pair is set to k = 0.01.

Below, we let [X] denote the concentration of sites in state X = CO, O, or E (empty) on the surface, e.g., the probability that a specific site is occupied by X = CO or O, or empty for X = E. Then, one has that [CO] + [O] + [E] = 1. Likewise, [X ⋅ Y] denotes the probability that a specific NN pair is occupied by X and Y, [X ⋅ Y ⋅ Z] that a specific linear triple is occupied by X, Y, and Z, etc. It is also convenient to introduce conditional probabilities or populations, where [X|Y] ≡ [X ⋅ Y]/[Y] is the conditional probability that a site is occupied by X given a NN Y site, [X|Y ⋅ Z] ≡ [X ⋅ Y ⋅ Z]/[Y ⋅ Z] is the conditional probability of an X given a neighboring Y ⋅ Z pair, etc. Behavior of the reaction model is then described exactly by the hierarchical rate equations

(1a)

and

(1b)

The positive gain terms are associated with adsorption, and the negative loss terms are associated with desorption and reaction. CO adsorption occurs at rate pCO [E], as a single empty site (which occurs with probability [E]) is required, and impingement occurs with rate pCO. O2 adsorption occurs at rate 2pO2 [E ⋅ E] reflecting the feature that an empty pair of sites (which occurs with probability [E ⋅ E]) is required, and that the impingement rate per site is 2pO2. Similar consideration applies for the corresponding desorption processes. For reaction, consider the rate at which a specific CO reacts. Since there are 2d neighboring sites which could be populated by O, the corresponding rate is 2dk [O⋅CO] where the pair probability reflects the requirement of an adjacent CO—O pair. Similar considerations yield the same reaction rate for O.

However, Eq. (1) are not closed as the presence of spatial correlations means that pair probabilities cannot be factorized as a product of single-site probabilities. Equations for the evolution of pair quantities couple to triples, etc., generating an infinite hierarchy. The turn-over-frequency (TOF) is defined as the reaction rate per site. The TOF can be regarded as the rate at which CO on a specific site reacts with a NN⋅O to form CO2. Equivalently, the TOF is the rate at which O on a specific site reacts with NN⋅CO. From (1), one has that TOF = 2dk [O ⋅ CO].

Model behavior is analyzed precisely by standard Kinetic Monte Carlo (KMC) simulation procedures. The surface is represented by a finite periodic lattice of adsorption sites with periodic boundary conditions. A site is selected at random, and one attempts to implement the various possible processes selected with probabilities proportional to their rates. Thus, the algorithm includes rejection, e.g., when adsorption is attempted on a site which is already populated. Below, we refer to Monte Carlo Steps (MCS), where 1 MCS corresponds selecting each site in the system on average once. Thus, 1 MCS corresponds to a physical time equaling the reciprocal of the sum of the rates for all distinct adsorption, desorption, and reaction processes.

Steady-state model behavior for ε ∼ 10−6 is shown in Fig. 1 for both 1D and 2D models. One complication is that even these extended simulations for 106 MCS still cannot access the true steady-state for very low pCO (below about 10−5 in 1D and 10−4 in 2D) for the large system sizes chosen. (Simulation can be achieved for longer times for smaller system sizes, but we wish to avoid any subtle finite-size effects for these non-equilibrium models.) Specifically, in this regime of very low pCO, Fig. 1 reveals a dependence on initial condition where we compare results starting from an initially empty and an initially O-covered surface. To obtain a more reliable and comprehensive picture of behavior, it is instructive to examine behavior for higher dO2. Thus, Fig. 2 shows steady-state behavior for the 1D model for dO2 = 10−4, 10−3, and 10−2 where simulation can more readily access the true steady-state. Trends in behavior for the 2D model with increasing dO2 (not shown) are similar. These data suggest three distinct regimes of behavior which we will describe immediately below, and analyze in more detail in Sec. IV.

FIG. 1.

Steady-state behavior for the reaction model with dCO = 0.1, k = 0.01, and pCO + pO2 = 1 from KMC simulation for 106 MCS. 1D model (d = 1, 216 site lattice) with dO2 = 10−6: (a) [E]; (b) TOF. 2D model (d = 2, 512 × 512 site lattice) with dO2 = 2 × 10−6: (c) [E]; (d) TOF. Solid (open) symbols correspond to an initially O-covered (empty) surface. Solid (dashed) curves denote the site (pair) approximation (see Sec. IV D).

FIG. 1.

Steady-state behavior for the reaction model with dCO = 0.1, k = 0.01, and pCO + pO2 = 1 from KMC simulation for 106 MCS. 1D model (d = 1, 216 site lattice) with dO2 = 10−6: (a) [E]; (b) TOF. 2D model (d = 2, 512 × 512 site lattice) with dO2 = 2 × 10−6: (c) [E]; (d) TOF. Solid (open) symbols correspond to an initially O-covered (empty) surface. Solid (dashed) curves denote the site (pair) approximation (see Sec. IV D).

Close modal
FIG. 2.

Steady-state behavior for the 1D model (d = 1) from KMC simulation (symbols) for higher dO2 values. Solid (dashed) curves show MF (pair) approximation predictions.

FIG. 2.

Steady-state behavior for the 1D model (d = 1) from KMC simulation (symbols) for higher dO2 values. Solid (dashed) curves show MF (pair) approximation predictions.

Close modal

Regime I for high pCO ≈ 1 corresponds to a primarily CO-populated surface. The dominant processes are single-site CO adsorption and desorption, thus producing a random CO adlayer with coverage determined by the CO adsorption-desorption balance. This simple behavior, for which MF kinetics naturally applies, was noted by Reuter et al. in their modeling,24,25 and will not be the focus of our analysis. We find a fairly sudden crossover through a peak in the TOF to a Regime II for lower pCO below about 0.1 which corresponds to a near O-covered surface with strong spatial correlations. Specifically, the conditional probability for an empty site given a NN empty site, [E|E], is of order unity (see Sec. IV), well above the MF value of [E] ≪ 1. We will also propose that a broad crossover occurs for very low pCO to a Regime III which is dominated by oxygen dimer adsorption-desorption. This regime will be shown to consist of a near O-covered surface, but now with a random distribution of O (and of vacancies and of CO with very low concentration). KMC simulation can readily access all of these regimes for higher dO2 of about 10−3 or above where they are compressed and shifted to higher pCO. See Fig. 2.

For a comprehensive understanding of the reaction model of Sec. II for low pCO, it is invaluable to develop a detailed characterization of limiting behavior of the pure oxygen adsorption-desorption limit with pCO = 0. We now provide this analysis.

While random adsorption-desorption of monomers does not induce spatial correlations and exhibits trivial exponential kinetics, this is not the case for adsorption-desorption of dimers on adjacent sites.26,27 The kinetics of oxygen adsorption-desorption is described by the exact evolution Eq. (1b) for the O-coverage, [O], but where now the reaction term is absent since pCO = 0 (so pO2 = 1) and there is no CO on the surface. Again, pair quantities cannot be factorized as products of single-site quantities due to spatial correlations, so these equations cannot be closed exactly.

Additional basic insight into the dimer adsorption-desorption model comes from dividing the lattice into two alternating sublattices labeled + and − as shown in Fig. 3. For a finite system with periodic boundary conditions and O sublattice populations N±, the difference M = N+ − N- is a conserved quantity. Thus, the system is not ergodic and exhibits disconnected dynamics in different “sectors” corresponding to different M-values.26,27 Just the sector M = 0 is relevant for our studies. This conservation constraint will significantly impact the dynamics, as discussed further below. An exact analysis of the time autocorrelation function was performed for the 1D version of this model showing slow diffusive decay.26,27 However, these previous studies did not analyze the evolution of the coverage or spatial correlations which are of primary interest for our study.

FIG. 3.

Sublattices on (a) a 1D linear lattice (d = 1) and (b) a 2D square lattice (d = 2).

FIG. 3.

Sublattices on (a) a 1D linear lattice (d = 1) and (b) a 2D square lattice (d = 2).

Close modal

A valuable alternative perspective on this dimer adsorption-desorption model and above-mentioned conservation law comes from noting that the model can be mapped onto a particle diffusion model.27 One simply interchanges occupied and empty sites on only the + sublattice. Then, one obtains a pure particle diffusion model where particles hop with rate pO2/d from + to NN − sites, and with rate dO2/d from − to NN + sites for the d = 1 or 2 dimensional model. Thus, conservation of M now translates into conservation of particle number which is given by Ld − M on a lattice with Ld sites. Also, it is clear that this diffusion model is described by a Hamiltonian with differing site energies for + and − sites. Thus, in the equilibrium steady-state, all configurations with the same population of the + sublattice (and thus the same population of the − sublattice) have equal probability. Translating this observation result back to the dimer adsorption-desorption model on a finite lattice, it implies that all ss configurations with the same coverage have the same probability.

The key conclusion of the above analysis is that for the relevant M = 0 sector of the dimer-adsorption-desorption model in the limit of an infinite system is that the steady-states involve a random distribution of O-populated (or empty) sites, i.e., there are no spatial correlations. This feature of the steady-state implies that [EE]ss = [E]ss2 and [OO]ss = [O]ss2, so together with the constraint [O]ss + [E]ss = 1, one can simply balance adsorption and desorption rates in (1b) to obtain the exact relation

(2)

and thus [O]ss = 1 − [E]ss = 1/(1 + ε1/2). This recovers a result from Ref. 27.

Returning to the issue of kinetics, Fig. 4 shows simulation results for the approach to the steady state when dO2 = 10−4 in 1D and 2D both starting from an O-covered surface (filled symbols) and starting from an empty surface (empty symbols). In both cases, there is a short transient regime followed by a very long quasi-steady-state (qss) regime before the final slow evolution to the true random steady-state.

FIG. 4.

KMC results for dimer adsorption-desorption kinetics for pO2 = 1 and dO2 = 10−4 where [E]ss = 0.0099. 1D model (d = 1, 216 site lattice) for (a) [E]; (b) [E|E]; (c) |[E] − [E]ss|. 2D model (d = 2, 512 × 512 site lattice) for (d) [E]; (e) [E|E]; (f) |[E] − [E]ss|. Solid (open) symbols correspond to an initially O-covered (empty) surface. Fits in (f) of the form et reveal that τ ≈ 1.8 × 105 (τ ≈ 0.83 × 105) for an initially O-covered (empty) surface. Dotted-dashed curves in (a), (b), (d), and (e) correspond to the pair approximation. Solid curves in (a) and (d) are a heuristic rate equation (HRE) treatment (see Sec. III E).

FIG. 4.

KMC results for dimer adsorption-desorption kinetics for pO2 = 1 and dO2 = 10−4 where [E]ss = 0.0099. 1D model (d = 1, 216 site lattice) for (a) [E]; (b) [E|E]; (c) |[E] − [E]ss|. 2D model (d = 2, 512 × 512 site lattice) for (d) [E]; (e) [E|E]; (f) |[E] − [E]ss|. Solid (open) symbols correspond to an initially O-covered (empty) surface. Fits in (f) of the form et reveal that τ ≈ 1.8 × 105 (τ ≈ 0.83 × 105) for an initially O-covered (empty) surface. Dotted-dashed curves in (a), (b), (d), and (e) correspond to the pair approximation. Solid curves in (a) and (d) are a heuristic rate equation (HRE) treatment (see Sec. III E).

Close modal

Starting from an O-covered surface, we propose that the initial transient states and quasi-steady state include strong correlations wherein most vacancies are incorporated in vacancy pairs. (We will see that these correlations are similar to those in the reaction model for low pCO.25) For dimension d = 1 or 2, this implies that

(3)

In addition, utilizing the exact relation [O ⋅ O] = 1 − 2[E] + [E ⋅ E], we conclude that [O ⋅ O] ≈ 1 − (4d − 1)[E]/(2d). Incorporating these relations into (1) yields

(4)

Thus, the quasi-steady-state value of [E] for large t is given by

and

(5)

This result consistent with our simulation analysis for ε = 10−4 which indicates that [E|E]qss ≈ 0.5 (0.25) and [E]qss ≈ 0.0002 (0.0004) for 1D (2D). We will elucidate the slow evolution from the quasi-steady-state to the true random steady-state in Sec. III E below.

Starting from an empty surface, formation of a quasi-steady-state corresponds roughly to Random Sequential Adsorption (RSA) of dimers.28–30 Here, the quasi-steady-state corresponds to the jammed state with respect to dimer adsorption and includes only isolated empty sites (with no empty pairs). The dynamics of exponential approach to this state for the 1D problem is given by28,30

so that

(6)

For the 2D problem, again there is an exponential approach to the quasi-steady-state where [E]qss ≈ 0.090 32.29,30 Simulations yield values of [E]qss consistent with the above results, and also reveal extremely small [E|E]qss ≈ 0.0006 (0.001) in 1D (2D).

To elucidate the role of the conservation law, i.e., invariance of M, in the dimer adsorption-desorption model in controlling kinetics and spatial correlations, it is instructive to recall that the dimer adsorption-desorption model can be mapped onto a particle diffusion model.27 For such a pure diffusion model, one expects diffusive decay of spatial correlations,31 and perhaps also in site occupancy kinetics.

For the special case where pO2 = dO2, the particle diffusion model corresponds to a randomly hopping gas of particles subject to exclusion of multiple site occupancy. Development of rate equations for the populations of + and − sites for this model shows that these exactly decouple from pair probabilities to form a closed set. See Appendix  B. This constitutes a special case of a more general analysis of Kutner for this problem.32 Solution of these equations reveals exponential decay of coverage to the steady-state, rather than the diffusive decay suggested in Ref. 27. However, evolution of two-point spatial pair probabilities or correlations is governed by an infinite coupled set of closed equations for different separations of the pair sites. These have the structure of diffusion equations in the space of site separations. Thus, evolution corresponds to diffusion of correlations from shorter to longer separations, and asymptotic behavior will reflect this diffusion process. In a time interval t, diffusion will spread correlations over a typical range ∼t1/2, corresponding to a volume ∼td/2 of separation space for dimension d. Thus, the magnitude of the nearest-neighbor and other short-range correlations should decay like t−d/2. This behavior also follows from a detailed spectral analysis of the linear evolution equations describing diffusion of correlations. See Ref. 31 and Appendix  B.

For pO2≠dO2, including the case of most relevance where ε = dO2/pO2 ≪ 1, the evolution of single-site quantities couples to that of the spatial pair correlations, which in turn couple to multi-point spatial correlations. See Appendix  B. While exact analysis of the pair correlations is likely not viable, one anticipates the same asymptotic diffusion behavior as in the above special case. (Diffusive behavior of the time-autocorrelation function for this 1D model was proved for pO2≠dO2.27) Thus, we claim that short-range spatial correlations should again decay like t−d/2 for a d-dimensional model. Furthermore, since the evolution of spatial correlations is coupled to that of coverages, we propose that the asymptotic approach to the steady-state of the coverage in dimer adsorption-desorption models should also exhibit diffusive decay of the form t−d/2.

For the case dO2 = 10−4 shown in Fig. 4, the true asymptotic approach to the steady-state is so slow that precise simulation analysis is not viable. Thus, to assess the above proposed behavior, we perform simulations with larger dO2 = 10−2. Results shown in Fig. 5 are consistent with the proposal that [E] − [E]ss ∼ t−d/2, as t → ∞.

FIG. 5.

KMC results for long-time dimer adsorption-desorption kinetics for pO2 = 1 and dO2 = 10−2. |[E] − [E]ss| versus t for (a) d = 1 (216 site lattice). Dashed line has slope −1/2. (b) d = 2 (512 × 512 site lattice). Dashed (dotted) line has slope −1 (−1/2). Solid (open) symbols correspond to an initially O-covered (empty) surface. [E] approaches [E]ss non-monotonically for an initially empty surface.

FIG. 5.

KMC results for long-time dimer adsorption-desorption kinetics for pO2 = 1 and dO2 = 10−2. |[E] − [E]ss| versus t for (a) d = 1 (216 site lattice). Dashed line has slope −1/2. (b) d = 2 (512 × 512 site lattice). Dashed (dotted) line has slope −1 (−1/2). Solid (open) symbols correspond to an initially O-covered (empty) surface. [E] approaches [E]ss non-monotonically for an initially empty surface.

Close modal

Starting with the exact evolution Eq. (1b), the MF site approximation to the kinetics makes the replacement [E ⋅ E] ≈ [E]2 and [O ⋅ O] ≈ [O]2 yielding

(7)

which corresponds to exponential decay to the steady-state with a time-scale τ MF = ε 1 / 2 / ( 4 p O 2 ) = 1 4 ( p O 2 d O 2 ) 1 / 2 . However, this treatment of kinetics fails to capture the quasi-steady-states and produces far too rapid evolution to the random steady-state.

In the pair approximation, we analyze evolution equations for the pair quantities [O⋅O], [O⋅E], and [E⋅E], noting that these equations include triplet probabilities. To close the equations, we employ the pair approximation [X ⋅ Y ⋅ Z] = [X ⋅ Y][Y ⋅ Z]/[Y], where X, Y, and Z = O, CO, or E. Then, one recovers [O] and [E] from conservation of probability relations. See Appendix  C. Like the MF approximation, starting with an O-covered surface, the pair approximation fails to recover the quasi-steady-state, but rather evolves directly to the true random steady-state. In Appendix  D, we trace this failure to the inability of the pair approximation to describe key spatial correlations. On the other hand, starting with an empty surface, the pair approximation does capture the quasi-steady-state which roughly corresponds to jammed state for RSA of dimers and includes only isolated empty sites (with no empty pairs). In fact, the pair approximation is known to be exact for RSA of dimers in 1D recovering [E]qss ≈ e−2 ≈ 0.135 34.28,30 For the 2D problem, the pair approximation predicts that [E]qss = 1/9, quite close to the exact value of 0.090 32.29,30 The dotted-dashed curves in Fig. 4 show pair-approximation predictions.

To assess asymptotic kinetics in the pair approximation, one linearizes the evolution equations about the steady state to reveal exponential decay with a characteristic time τpair = O(ε−1/pO2) = O(dO2−1) which is fundamentally slower than in the site approximation, but still does not reflect actual model behavior. See Appendix  C.

In addition to traditional MF and pair approximations, it is instructive to develop heuristic rate equations tailored to the regime ε ≪ 1 to describe the evolution of key populations related to empty sites starting with an O-covered surface. In this regime, we anticipate that empty sites are either isolated with probability, [E], or incorporated into empty pairs where [E ⋅ E] is the probability that a specific pair is empty and isolated (i.e., there are no empty triples, etc.). Thus, one has that [E] ≈ [E] + 2d[E ⋅ E]. Below, we first discuss the development of heuristic rate equations for general dimension d. However, we focus on d = 2 where the theory should be most effective. Possible refinements for d = 1 are briefly discussed.

First, consider the evolution of [E]. Isolated empty sites can be created when (i) a slow O2 desorption event occurs at rate dO2 adjacent to an existing empty pair on the surface thereby creating a connected empty linear or bent quartet and (ii) subsequent (fast) adsorption of O2 on the central pair of sites in the quartet creates a separated pair of isolated vacancies. See Fig. 6(a). Thus, the creation rate has the form Rgain = AdO2 [E ⋅ E]. Isolated empty sites effectively diffuse when (i) a slow O2 desorption event at rate dO2 adjacent to the isolated empty site to create an empty triple and (ii) subsequent fast adsorption on the original empty site and a neighbor creates an isolated empty site at a different location. See Fig. 6(b). Isolated empty sites are destroyed by (i) diffusion-mediated formation of empty pairs from two isolated vacancies and (ii) subsequent (fast) adsorption of O2 to annihilate this pair. Thus, the destruction rate has the form Rloss = BdO2([E])2. We conclude that

(8)

The prefactors, A and B, can be assessed by counting the number of configurations associated with the above processes. For effective diffusion coefficient of isolated vacancies, one should also account for distances between the initial and final vacancies in Fig. 6(b). However, we do not provide such a detailed analysis here, and will also set A ≈ B to recover the correct steady-state behavior.

FIG. 6.

2D examples (d = 2) of (a) creation of isolated empty sites E* from isolated empty pairs, EE; (b) effective diffusion of isolated empty sites, E (with both processes mediated by O2 desorption-adsorption). Analogous processes occur in 1D (d = 1).

FIG. 6.

2D examples (d = 2) of (a) creation of isolated empty sites E* from isolated empty pairs, EE; (b) effective diffusion of isolated empty sites, E (with both processes mediated by O2 desorption-adsorption). Analogous processes occur in 1D (d = 1).

Close modal

Second, we consider dominant contributions to the evolution of [E ⋅ E]. Creation of isolated empty pairs occurs via desorption of O2 at rate dO2/d from most of the surface, and annihilation by adsorption of O2 at rate pO2/d. Thus, one has that

(9)

Using Eqs. (8) and (9) to analyze evolution starting from an O-covered surface with [E] = [E ⋅ E] = 0, it is clear that one quickly achieves a quasi-steady-state with [E ⋅ E] ≈ dO2/pO2 = ε. Then, the solution to (8) with A = B becomes

where

(10)

This heuristic result for the kinetics is shown for the 2D model in Fig. 4(d) (solid curve) with A = 1. We find that the life-time of the quasi-steady-state is correctly described by τrlx = ε−3/2/(ApO2), but not by τMF = ε−1/2/(4pO2) or τpair = O(ε−1/pO2). This result is further confirmed by the analysis of (non-asymptotic) exponential decay towards the steady-state in Fig. 4(f). Again, all of these treatments predict exponential decay to the steady-state rather than the actual true asymptotic diffusive decay. For the 1D model, the above heuristic theory is also reasonable as shown in Fig. 4(a) (solid curve). However, there are some additional limitations and subtleties due to the strong recurrence of random walks in 1D. In fact, a suitably refined rate equation treatment in Appendix  E indicates a longer quasi-steady-state lifetime scaling like τrlx = O(ε−2/pO2). The actual kinetics should reflect a spectrum of lifetimes between τrlx and τrlx.

In Sec. II, we have already described the different regimes of steady-state behavior for our reaction model upon varying pCO. Regime I for high pCO ≈ 1 corresponds to a primarily CO-populated surface. Regime II for lower pCO below about 0.1 corresponds to a near O-covered surface with strong spatial correlations, and Regime III for very low pCO corresponding to a random near-O-covered surface. Here, we provide a more detailed analysis (focusing on the latter two regimes) including crossover behavior.

As noted in Refs. 12–15 and in Sec. II, analysis of Regime I with “high” pCO ≈ 1 and pO2 ≪ 1 is perhaps most trivial as this is dominated by random adsorption and desorption of CO on a predominantly randomly CO-populated surface. Thus, CO adsorption-desorption equilibrium implies that [CO]ss ≈ pCO/(pCO + dCO) ≈ 0.909 and that [E]ss ≈ dCO/(pCO + dCO) ≈ 0.091. With pO2 ≪ 1, rare dissociative adsorption and reaction of oxygen on this surface will not significantly perturb the adlayer statistics. Thus, we calculate the TOF from the rate of reaction of O at a specific site as a product of (i) the rate of adsorption of O on that site which is the product of pO2/d times [E ⋅ E]ss ≈ ([E]ss)2 times the number (2d) of empty pairs overlapping the site and (ii) the probability, Pk1 ≈ (2d − 1)k/[(2d − 1)k + pO2] ≈ 1 (for dO2 ≪ k) that this O reacts with a NN⋅CO rather than undergoing recombinative desorption. Thus, one has that

(11)

which is consistent with a MF rate equation treatment for dO2 ≪ k.

With regard to the behavior for a broader range of pCO including crossover to the mixed adlayer with maximum TOF, and subsequent transition to a near O-covered surface, the MF approximation is generally not adequate. However, the standard pair approximation described in Sec. IV D can reasonably describe this behavior.

For Regime III with very low pCO and an essentially random O-covered steady-state, CO adsorption-desorption and reaction should constitute a negligible perturbation on the adlayer statistics. Furthermore, since the rate of CO adsorption on mainly isolated empty sites with concentration [E] ≈ [E]ss ≈ ε1/2 is pCO, and the probability that such a CO reacts with adjacent O rather than desorbing is Pk2 = 2dk/(dCO + 2dk), it follows that

(12)

i.e., TOF ≈ 0.167ε1/2 pCO in 1D, and TOF ≈ 0.286 ε1/2 pCO in 2D. The same form follows from an analysis of the MF rate equations. The first-order dependence of TOF on pCO is clear in the log-log plots of TOF shown earlier.

Here, we consider Regime II with pCO below about 0.1 where the surface becomes mainly O-covered. A key observation by Matera, Meskine, and K. Reuter (MMR)25 for the analogous regime in their model of CO-oxidation on RuO2(110) is that strong correlations in the adlayer are associated with the feature that empty sites are mainly created in NN pairs. The dominance of vacancy pairs, which implies that [E|E] = 1/(2d), was attributed to their creation by dimer desorption and CO + O reaction.

A refined rationalization of this behavior derives from our analysis of the pure dimer adsorption-desorption process in Sec. III. In that case, strong correlations in the quasi-steady-state associated with vacancy pair creation by dimer desorption are ultimately destroyed by the slow creation of isolated empty sites. However, introduction in the reaction model of CO adsorption upon these isolated empty sites can lead to reaction with an adjacent O creating an empty pair. This empty pair is then promptly destroyed by O2 adsorption. This mechanism for “reactive removal of isolated empty sites” should at least partly sustain a strong correlation associated with a high population of empty pairs (a feature which is ultimately destroyed in the pure dimer adsorption-desorption model for pCO = 0). Indeed, our heuristic rate equation treatment in Sec. IV F based upon this picture does imply crossover with increasing pCO from the random near-O-covered Regime III to a regime where empty sites are predominantly in pairs and [E|E] = 1/(2d).

However, actual behavior is more complicated. In Fig. 7, we show KMC simulation results for dO2 = 10−4 of the conditional probability [E|E] starting with an O-covered surface where the simulations are run for various time periods. Significantly, simulations run for shorter time periods do not reach the true steady state, yielding artificially high values of [E|E] closer to the value of 1/(2d). One needs simulation times of at least t ∼ 106 MCS (105 MCS) for d = 1 (d = 2) to access true steady-state behavior. (This difficulty in assessing the true steady-state is reminiscent of the very slow diffusive dynamics seen in the pure dimer adsorption-desorption model.) The actual behavior obtained from precise KMC analysis reveals steady-state values of [E|E] which are O(1) for a broad range of pCO, and only approach the MF value of [E] for very low pCO corresponding to Regime III. However, the maximum value is about 0.3 (0.1) for d = 1 (d = 2) contrasting the value of 0.5 (0.25) associated with only isolated empty site pairs.

FIG. 7.

Conditional probability [E|E] for (a) d = 1 with dO2 = 10−4 and (b) d = 2 with dO2 = 2 × 10−4. Behavior is shown for various simulation times t, with true steady-state behavior achieved only for long t. Solid (dashed) lines indicate MF (pair) approximation predictions. The horizontal dotted lines represent the prediction of Matera et al.25 of [E|E] = 0.5 in 1D for (a), and its extension to 2D of [E|E] = 0.25 for (b).

FIG. 7.

Conditional probability [E|E] for (a) d = 1 with dO2 = 10−4 and (b) d = 2 with dO2 = 2 × 10−4. Behavior is shown for various simulation times t, with true steady-state behavior achieved only for long t. Solid (dashed) lines indicate MF (pair) approximation predictions. The horizontal dotted lines represent the prediction of Matera et al.25 of [E|E] = 0.5 in 1D for (a), and its extension to 2D of [E|E] = 0.25 for (b).

Close modal

Here, we compare precise KMC results for steady-state behavior with predictions of standard mean-field (site) and pair approximations to the hierarchical rate equations and with those from a variety of refined approximations. The standard MF approximation simply factorizes the pair population terms appearing in (1) as a product of single-site populations or concentrations, e.g., [E ⋅ E] ≈ [E]2. For the standard PA, we start with exact rate equations for pair quantities. For example, for d = 1, one has

(13)

where the gain terms correspond to O2 adsorption populating both of the pair of sites under consideration or just one of two where the other is already populated by O. Loss terms correspond either to desorption of O2 or reaction removing one or both of the O in the pair of sites under consideration. One can similarly develop equations for [O⋅CO], [CO⋅CO], [O⋅E], [CO⋅E], and [E⋅E] pairs, where conservation of probability means that one of these quantities is redundant, and also single-site quantities can be obtained from, e.g., [O] = [O ⋅ O] + [O ⋅ CO] + [O ⋅ E]. Closure of these rate equations follows from the pair approximation [X ⋅ Y ⋅ Z] ≈ [X ⋅ Y][Y ⋅ Z]/[Y] for a linear triple as already noted in Sec. III D. This factorization approximation can also be concisely formulated in terms of conditional probabilities as [X|Y ⋅ Z] ≈ [X|Y], where this formulation is more amenable to systematic extension to higher-order approximations.29,30 Results shown for d = 1(d = 2) in Fig. 8 (Fig. 9) for dO2 = 10−4 (and also previously in Figs. 1 and 2) indicate that the pair approximation significantly improves on the MF approximation, although it can still be far from the recovering exact behavior.

FIG. 8.

Comparison with KMC simulation results (solid circles) for [E] (top frame) and TOF (bottom frame) in 1D (d = 1) with dO2 = 10−4 with various analytic approximations: mean-field (MF); standard pair approximation (PA); modified pair factorization (PA1) in Sec. IV D; approximation of MMR;25 refined pair approximation (PA2) of Sec. IV E; and HRE of Sec. IV F.

FIG. 8.

Comparison with KMC simulation results (solid circles) for [E] (top frame) and TOF (bottom frame) in 1D (d = 1) with dO2 = 10−4 with various analytic approximations: mean-field (MF); standard pair approximation (PA); modified pair factorization (PA1) in Sec. IV D; approximation of MMR;25 refined pair approximation (PA2) of Sec. IV E; and HRE of Sec. IV F.

Close modal
FIG. 9.

Comparison with KMC simulation results (solid circles) for [E] (top frame) and TOF (bottom frame) in 2D (d = 2) with dO2 = 2 × 10−4 with various analytic approximations: MF; PA; PA1; PA2; MMR; and HRE. See Fig. 8 for definitions.

FIG. 9.

Comparison with KMC simulation results (solid circles) for [E] (top frame) and TOF (bottom frame) in 2D (d = 2) with dO2 = 2 × 10−4 with various analytic approximations: MF; PA; PA1; PA2; MMR; and HRE. See Fig. 8 for definitions.

Close modal

Some simple refinements to the PA might be motivated by the feature that in Regime II, there is a predominance of vacancy pairs and relatively few isolated empty sites, empty triples, and larger clusters of empty sites. Note that in this case, the standard PA can greatly overestimate the population of vacancy triples as [E ⋅ E ⋅ E] ≈ [E ⋅ E]2/[E] which is the same order of magnitude as [E ⋅ E]. This suggests instead implementing the simpler factorization [E ⋅ E ⋅ E] ≈ [E][E ⋅ E] leading to a much reduced population of empty triplets. Specifically, we implement a slightly modified pair approximation, denoted by PA1, utilizing multiple modified factorizations

and

(14)

and using the standard PA for other triplet populations. Note that this approximation enforces the exact condition that [E ⋅ E ⋅ E] + [O ⋅ E ⋅ E] + [CO ⋅ E ⋅ E] = [E ⋅ E]. PA1 deviates significantly from the standard PA and overall is significantly more accurate for very low pCO. See Figs. 8 and 9.

Finally, we also describe results from the approximation implemented by MMR25 for 1D and the analogous extension to 2D based on the proposal vacancies occur exclusively in NN pairs. This assumption implies that [EE] = [E]/2 for d = 1, which we extend to [E ⋅ E] = [E]/(2d) for general d. In terms of conditional probabilities, the MMR approximation becomes [E|E] = 1/(2d). MMR also proposed that most CO will have adsorbed within a vacancy pair, and thus have exactly one O-neighbor so that [ O CO ] = 1 2 [ CO ] . We naturally extend this relation to claim that [ O CO ] = 3 4 [ CO ] for d = 2. If we also use the MF relation [O ⋅ O] ≈ [O]2, then we can close Eq. (1) for single-site concentrations. Results in Figs. 8 and 9 reveal that the MMR approximation significantly improves upon the MF, PA, and PA1 approximations for d = 1, except for very small pO2. The improvement is not so significant for d = 2, although we note that the effectiveness of the MMR generally improves for lower dO2.

It is clear from our above analysis of various truncation approximations, and from the earlier work of Matera et al.,25 that an accurate treatment of the empty site statistics is the key to a reliable description of reactivity. Thus, we are motivated to improve on the descriptions in Sec. IV D. To this end, we will consider the conditional probability, [E|E ⋅ E], to find an empty site adjacent to a specified empty pair, where both linear and bent configurations are considered for d = 2. The MF approximation sets equal [E|E ⋅ E] ≈ [E], the standard PA sets [E|E ⋅ E] = [E|E], and the MMR treatment sets [E|E ⋅ E] = 0.

For a more precise description, we propose and identify two distinct contributions to [E|E ⋅ E]. The first contribution is not associated with reaction and constitutes a refinement to the MMR picture. Rather than consider essentially all empty sites as being incorporated into isolated empty pairs, one can consider empty sites as being described by a “random dimer-vacancy gas” (DVG). Within this picture, one expects a contribution, [E|E ⋅ E]DV G ≈ [E] to the conditional probability since there is a finite probability for two dimer-vacancies to be adjacent to each other. The key observation is that this relation applies more generally than for random distributions of empty sites. The second contribution is associated with reaction. Specifically, we consider the scenario where a CO adsorbs on one site of an empty pair and reacts with a specific adjacent O rather than reacting with another O or desorbing to create an empty triple. The former occurs at rate pCO, and the latter occurs with probability Pk4 = k/(k + dCO) for d = 1, and Pk4 = k/(3k + dCO) for d = 2. Also, since the lifetime of an empty pair is given by τ2E = d/pO2, the fraction of time for which the additional site is empty is given by [E|E ⋅ E]|RXN = pCOPk4τ2E. See Appendix  F. Thus, we have that

(15)

Fig. 10 shows the high accuracy of this formulation compared with the MF and pair approximations for dO2 = 10−4 and d = 1. For an analysis of steady-state behavior in the reaction model, we implement an approximation PA2 which uses (15), the factorizations [O ⋅ E ⋅ E] = [O] [E ⋅ E], [CO ⋅ E ⋅ E] = [CO][E ⋅ E], and also a standard pair approximation for other triples.33 The results of PA2 shown in Figs. 8 and 9 significantly improve upon all other approximations described previously (MF, PA, PA1, and MMR).

FIG. 10.

Conditional probability [E|E ⋅ E] for d = 1 with dO2 = 10−4 compared with the MF and PA1 approximations, [E|E ⋅ E] = [E] (dotted curve), the PA, [E|E ⋅ E] = [E|E] (dashed curve), and our refined PA2 treatment of Sec. IV E (solid curve).

FIG. 10.

Conditional probability [E|E ⋅ E] for d = 1 with dO2 = 10−4 compared with the MF and PA1 approximations, [E|E ⋅ E] = [E] (dotted curve), the PA, [E|E ⋅ E] = [E|E] (dashed curve), and our refined PA2 treatment of Sec. IV E (solid curve).

Close modal

An alternative strategy to describe steady-state behavior, and in particular to elucidate the broad crossover from Regime III to Regime II, is to develop a HRE treatment in the spirit of Sec. III E. In this treatment, we augment Eq. (8) to account for pathways for creation and destruction of isolated empty sites, E, associated with the CO + O reaction process. Isolated empty sites can be destroyed by CO-adsorption, followed by reaction with adjacent O to create an empty pair, and subsequent adsorption of O2 to fill this empty pair. The rate of CO-adsorption onto an isolated empty site is pCO, and the probability that this CO reacts with an adjacent O rather than first desorbing is given by Pk2 = 2dk/(dCO + 2dk). Thus, the reactive loss rate has the form Rloss(rxn) = pCOPk2[E]. Isolated vacancies can also be created via reaction from vacancy pairs: a CO adsorbs at rate pCO, reacts before desorbing with probability Pk3 = (2d − 1)k/[dCO + (2d − 1)k] with a NN O to create a vacancy triple, and subsequent O2 adsorption produces an isolated vacancy. Thus, the reactive gain rate has the form Rgain(rxn) = 2pCOPk3[EE]. Choosing prefactors A = B and augmenting (8) by these terms yield a heuristic rate equation for [E] having the form

(16)

We supplement this equation using (9) to describe the evolution of [EE], or by simply assuming that [EE] ≈ ε is controlled by dimer adsorption-desorption quasi-equilibrium as in Sec. III E. Then, steady-state analysis of (15) yields

(17)

To determine the corresponding TOF, we use TOF ≈ pCOPk[E]ss where [E]ss = [E] + 2d[E ⋅ E]. Results from this HRE treatment shown in Figs. 8 and 9 are generally as accurate as our most successful refined pair-approximation treatment, PA2.

Finally, we utilize (16) to assess crossover from Regime III to Regime II. The lower boundary, pCO = pCO, of the crossover regime will correspond to the onset of significant deviation of [E]ss from ε1/2 when (pCOPk2)/(2ApO2ε) = O(ε1/2). The upper boundary, pCO = pCO+, of the crossover regime is less precisely determined. One criterion to determine pCO+ is to require that [E]ss be reduced from ε1/2 to O(ε) implying that [E]ss ≈ (ApO2ε2)/(pCOPk2) = O(ε). In conclusion, these analyses suggest that the lower and upper boundaries of the crossover region satisfy

and

(18)

These predictions are compatible with the simulation results in Sec. II (most clearly for those in Fig. 3 for dO2 = 10−3).

Finally, since actual reaction systems at high-P will have at least some degree of surface mobility, it is instructive to explore the effect on steady-state behavior of mobility in our reaction model. Specifically, we introduce hopping of both adspecies to NN empty sites at rate h. Naturally, increasing h which induces adlayer mixing and reduces spatial correlations, and consequently behavior approaches that of the MF prediction. See Fig. 11. However, a significant feature is that the rate of approach to MF behavior with increasing h is much slower for d = 1 than for d = 2. This is expected from the analysis of other 1D reaction-diffusion systems where single-file diffusion (hopping of reactants and/or products to NN empty sites without exchange) has limited ability to reduce strong spatial correlations induced by reaction.34 

FIG. 11.

KMC results for [E] versus pCO showing the effect of adspecies mobility with hop rate h on steady-state behavior. (a) 1D model with dO2 = 10−4. (b) 2D model with dO2 = 2 × 10−4. Solid curve: h-independent MF prediction. Dashed curve: h = 0 PA prediction. PA behavior for increasing h > 0 will smoothly have a transition from h = 0 PA to MF.

FIG. 11.

KMC results for [E] versus pCO showing the effect of adspecies mobility with hop rate h on steady-state behavior. (a) 1D model with dO2 = 10−4. (b) 2D model with dO2 = 2 × 10−4. Solid curve: h-independent MF prediction. Dashed curve: h = 0 PA prediction. PA behavior for increasing h > 0 will smoothly have a transition from h = 0 PA to MF.

Close modal

The simple model considered here for CO-oxidation under high-P conditions was motivated by the modeling of Reuter and coworkers of this reaction of RuO2(110).12,13,23–25 The model exhibits rich behavior. Not only is there a transition from a random CO-populated steady-state for high CO partial pressure, pCO, to a strongly correlated near-O-covered steady-state for low pCO as can be anticipated from the previous similar modeling.24,25 In addition, there is a further transition to a random near-O-covered steady-state for very low pCO. The subtle and slow diffusive kinetics for the pure oxygen adsorption-desorption regime at pCO = 0 impacts behavior for low pCO. This slow kinetics can make it difficult to reliably access the true steady-state via KMC simulation.

Particularly, challenging is a characterization of strong spatial correlations in the near-O-covered steady-state for low pCO as these are kinetically rather than thermodynamically induced. The novel treatment of Matera et al.25 based on the idea that most empty sites occur in pairs greatly improves on MF and PA treatments at least for the 1D model. However, significantly more accurate treatment follows from our refined pair-approximation PA2 describing subtle correlations controlling the population of triples of empty sites. Similar success was achieved with a quite different HRE treatment.

We have focused on the behavior of our simple CO-oxidation model for the specific rate choice pCO + pO2 = 1, dCO = 0.1, k = 0.01 (the reaction-limited regime with low rate of reaction relative to adsorption), with no surface diffusion, often for very low dO2 of 10−4 or below. We have briefly discussed the effect of introducing surface diffusion which is as expected, i.e., randomization of the adlayer and recovery of MF type behavior.

However, it is also appropriate to comment on anticipated behavior of our model in other rate regimes, particularly for lower dCO which could be regarded as corresponding to lower temperature. From this perspective, it is appropriate to recall the classic Ziff-Gulari-Barshad (ZGB) model for d = 235 which includes the same ingredients as our reaction model except that dCO = dO2 = 0 and k = ∞ (which corresponds to the adsorption-limited regime with instantaneous reaction). This model includes a continuous O-poisoning transition which is immediately removed by selecting dO2 > 0, and thus is not germane to our interests. Of more relevance is the feature that the ZGB model for d = 2 exhibits a discontinuous CO-poisoning transition which is preserved for finite k36 and also for dCO > 0 below a finite (small) critical value.37 Thus, we expect a discontinuous CO-poisoning transition in our current reaction model for d = 2 with very low or zero dO2 provided that one sufficiently reduces dCO.38 Since we operate in the reaction-limited rather than adsorption-limited regime (with k = 0.01) with negligible desorption, reaction creates rare empty site pairs which are filled by CO with rate 2pCO or with O2 at rate pO2/2. Thus, adsorption balance requires that the transition occurs at around p CO / p O 2 = 1 4 . See Ref. 36 for a more detailed discussion noting that behavior in this regime is similar to that of the Voter model.39 

Finally, we remark that the general strategies developed here to characterize kinetically induced spatial correlations should extend beyond CO-oxidation on RuO2(110) to other catalyst surfaces and other reactions in the high-pressure regime. Examples where KMC modeling already exists might include CO-oxidation on a Pd oxide surface,14,15 but also selective oxidation of ammonia,40,41 and oxidation of HCl42 on RuO2(110). It should, however, be noted that effective analytic descriptions of the spatial correlations in other reaction systems will likely need to be tailored to the specific reaction mechanism rather than using some generic formulation.

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 and CTC programs. The work was performed at Ames Laboratory which is operated for the USDOE by Iowa State University under Contract No. DE-AC02-07CH11358.

Figure 12 shows the structure of the perfect RuO2(110) surface. This state consists of alternating rows or columns of so-called coordinatively unsaturated (cus) sites on top of Ru atoms, and bridge (br) site between pairs of Ru atoms. The br sites are populated by O in this state, but these can be removed by reaction with CO. As noted in Sec. I, CO adsorption-desorption occurs at single sites of either type, and oxygen adsorption-desorption and CO + O reaction occur at adjacent pairs of sites of any type. In an oxygen-rich environment, both cus and br sites are populated by O. In CO-rich environments, the cus sites are populated by CO, and the O on the br sites is removed and replaced by CO.

FIG. 12.

Arrangement of cus and br sites on the RuO2(110) surface. Reprinted with permission from K. Reuter and M. Scheffler, Phys. Rev. B 73, 045443 (2006). Copyright 2006 American Physical Society.

FIG. 12.

Arrangement of cus and br sites on the RuO2(110) surface. Reprinted with permission from K. Reuter and M. Scheffler, Phys. Rev. B 73, 045443 (2006). Copyright 2006 American Physical Society.

Close modal

As noted in Sec. III B, dimer adsorption-desorption is equivalent to a particle diffusion model wherein particles hop from the + to NN − sublattice sites at rate pO2/d, and from − to NN + sublattice sites at rate dO2/d. Let [X+] ([X]) denote the particle population on the + (−) sublattice. Also, [E+X] denotes the probability of an adjacent filled + site and empty − site, etc. Then, for either the 1D or 2D model, one has

(B1)

and similarly for [X]. From these equations, it is immediately clear that [X+] + [X] is constant. For the special case where pO2 = dO2, there is an exact decoupling of the evolution of the single-site quantities from the pair quantities. The resulting equations reveal simple exponential decay to the steady-state in both 1D and 2D.

For the 1D model, let [X+ − X+]n denote the probability of two filled + sublattice sites separated by an even number, n ≥ 2, of lattice constants. Other pair probabilities, [X − X]n for even n ≥ 2, and [X+ − X]n = [X − X+]n for odd n ≥ 1 (so [X+X]1 = [X+X]), are similarly defined. In addition, we define [XE+ − X+]n = [X+ − E+X]n as the probability for an XE+ pair with E+ separated by odd n ≥ 1 lattice constants from an X+, etc. Then, for even n ≥ 2, one has that

(B2)

Other pair quantities satisfy similar equations for n > 1. To close these equations, an additional equation is needed for adjacent pairs of filled sites corresponding to n = 1.

For the special case where pO2 = dO2, there is an exact decoupling of the evolution of the pair probabilities from that of other quantities. However, characterization of the solutions of this semi-infinite coupled linear set of equations requires spectral analysis of the corresponding linear evolution operator (a semi-infinite matrix). To this end, one looks for special solutions to the above equations for n > 1 of the form31,43

and

(B3)

so that μ corresponds to the eigenvalue of the linear evolution operator. Substitution into the evolution equations for n > 1 with pO2 = dO2 = 1 yields the condition

(B4)

Of the resulting branches of eigenvalues, μj(q), with j = 0-2, one branch j = 0 satisfies μ0 = 0 when q = 0, and more generally μ0(q) ∼ − Aq2 for small q with c > 0. This branch of eigenfunctions with μ0(q) → 0, as q → 0 is analogous to acoustic or Goldstone modes in different applications.26,27 Construction of special solutions for the complete set of evolution equations for n ≥ 1 requires taking linear combinations of the above solutions for e±iqn which can also be regarded as determining the “phase shift” associated with satisfying the boundary condition (B3).31,43 However, the above form of the eigenvalues μj(q) is unchanged.

The complete solution of the initial value problem, we introduce eigenvectors, ϕ ¯ j ( q ) , components of which give the different pair probabilities for various n ≥ 0. Then, the general solution has the form j = 0 - 2 π < q < + π dq c j ( q ) exp [ μ j ( q ) t ] ϕ ¯ j ( q ) . The j = 0 term for large t can be approximated as

(B5)

revealing diffusive decay.

For the 2D model, let [X+ ∼ ∼ X+]n,m denote the probability of two filled + sublattice sites separated by n (m) lattice constants in the x-direction (y-direction), where n + m is even. Other pair probabilities, [X ∼ ∼ X]n,m for even n + m, and [X+ ∼ ∼ X]n,m = [X ∼ ∼ X+]n,m for odd n + m, are similarly defined. The equations for these pair probabilities have a generic form for n + m > 1 and are closed by an additional equation for n + m = 1 for the probability of adjacent pairs of filled site [X+X]n,m with n + m = 1. For a spectral analysis of these equations or the associated linear operator, one looks for special solutions of the generic equations for n + m > 1 of the form

(B6)

where μ is the eigenvalue of the linear evolution operator. Substitution into the evolution equations for n + m > 1 with pO2 = dO2 = 1 yields the condition

(B7)

Of the 3 branches of eigenvalues, μj(qx, qy), with j = 0, 1, 2, one Goldstone branch j = 0 satisfies μ0(q) ∼ − B[(qx)2 + (qy)2] for small q with c > 0. Special solutions for the complete set of evolution equations for n + m ≥ 1 involve linear combinations of the above solutions. Analogous to the above 1D analysis, long-time nature of the solutions comes from

(B8)

In the pair approximation for dimer adsorption-desorption, we can choose to retain two of [O⋅O], [O⋅E], and [E⋅E] as independent variables. Then, one recovers the third pair quantity and also [O] and [E] from conservation of probability relations

and

(C1)

The rate equations for pair quantities for the 1D problem have the form

(C2a)
(C2b)

and

(C2c)

where [E⋅E⋅E] represents the probability of an adjacent triple of empty sites, etc. The gain (loss) terms correspond to different possibilities for creating (destroying) the configurations on the left by adsorption or desorption. To close the equations, we employ the pair approximation [E ⋅ E ⋅ E] = [E ⋅ E]2/[E], [O ⋅ E ⋅ E] = [O ⋅ E][E ⋅ E]/[E], etc.

For the 2D problem, the equations have the same form except that now both linear and bent triplet probabilities appear although these become equal in the pair approximation. Also, various terms have different weights. The pair terms on the right-hand-side (RHS) of the d/dt [O ⋅ O] equation have a weight of 1 2 and the triple terms a weight of 3. All terms on the RHS of the d/dt [O ⋅ E] equation have a weight of 3/2.

To assess asymptotic kinetics in the pair approximation for dimer adsorption-desorption, one linearizes the rate Eq. (C2) about the steady state via [O ⋅ O] = [O ⋅ O]ss + δ[O ⋅ O] and [O ⋅ E] = [O ⋅ E]ss + δ[O ⋅ E], where [O ⋅ O]ss = ([O]ss)2 and [O ⋅ E]ss = [O]ss[E]ss. Forming a 2 × 1 vector δ ¯ with entries δ1 = δ[O ⋅ O] and δ2 = δ[O ⋅ E], the linearized evolution equations have the form

(C3)

The entries of J = for the 1D problem are

(C4a)
(C4b)

An analysis of the eigenvalues λi of J = reveals that λ1 = − 5pO2 + O(pO2ε1/2) and negative λ2 = O(pO2ε). Thus, the key observation is that terms O(ε1/2) are absent in λ2 due to cancellation. Thus, asymptotic decay to the steady-state which is controlled by λ2 occurs fundamentally slower than in the MF approximation with a characteristic time τpair = O(ε−1/pO2) = O(dO2−1).

Exactly, the same scenario applies for the pair approximation analysis for the 2D problem with τpair = O(ε−1/pO2) = O(dO2−1), but now

(C5a)
(C5b)

As noted in Sec. III D, the standard pair approximation fails to capture the quasi-steady-state for dimer adsorption-desorption obtained starting with an O-covered surface. To elucidate this failure, we present a simplified version of the exact rate equations and of the pair approximation tailored to a near O-covered surface. For the 1D problem, using [O] ≈ [O ⋅ O] ≈ [O ⋅ O ⋅ O] ≈ 1, the pair equations can be simplified to

(D1)

The standard pair approximation applies [E ⋅ E ⋅ E] ≈ [E ⋅ E]2/[E] and [OEE] ≈ [OE][EE]/[E] with [E] = [OE] + [EE]. Then, from the structure of the [E⋅E] equation, it is clear that integration starting from an O-covered surface where [E] = [EE] = [OE] = 0 leads to saturation of [EE] at a value of order ε as is appropriate for the quasi-steady-state. However, the structure of the [O⋅E] equation does not lead to any saturation of [O⋅E] at a value of order ε. Instead, [O⋅E] and thus [E] quickly increases through this value to the higher true steady-state value of around ε1/2.

The pair approximation factorization fails to capture the actual behavior [E ⋅ E ⋅ E] ≪ [E ⋅ E] and [O ⋅ E ⋅ E] ≈ [O ⋅ E] in the quasi-steady state. From these relations, the above equations immediately show that [E ⋅ E]qss ≈ [O ⋅ E]qss ≈ ε, so that [E]qss ≈ 2ε for the 1D system, consistent with the analysis in Sec. III A. For the 2D system, basically the same scenario applies. The actual behavior in the quasi-steady-state is described by [E ⋅ E ⋅ E] ≪ [E ⋅ E] for bent and linear triples of empty sites, and [O ⋅ E ⋅ E] ≈ [O ⋅ E]/3 for bent and linear triples. Substituting these results into the rate equations for [O⋅E] and [E⋅E] yields [EE]qss ≈ ε and [OE]qss ≈ 3ε, so that [E]qss ≈ 4ε consistent with Sec. III A.

With regard to heuristic rate equation treatments for dimer adsorption-desorption, we note in Sec. III E some subtleties for the 1D model due to the strong recurrence of random walks. The simplest formulation simply also adopts (8) for [E] and (9) for [E ⋅ E]. Results shown in Fig. 4(a) (solid curve) are reasonable but those in Fig. 4(c) suggest an overestimate of the rate of approach to the true steady-state in the regime. Indeed, (8) fails to account for the feature that diffusion-mediated annihilation of isolated vacancies is intrinsically slower for 1D diffusion (versus 2D diffusion), and also just-created nearby isolated vacancies are much more likely to quickly recombine for the 1D model.

Thus, we develop a refined treatment of diffusion-mediated annihilation and of creation of isolated vacancies. This formulation is more flexible in capturing the differences in behavior for 1D versus 2D diffusion. For diffusion-mediated annihilation of isolated vacancies by reaching sites adjacent to other vacancies followed by dimer adsorption, we write the rate as R(loss) = [E]/τ, where τ is the lifetime of such vacancies before annihilation.44 τ is estimated as the lifetime of a random walker amongst a random distribution of trap sites with density [E]. Thus, one obtains44,45

but

(E1)

For creation of isolated vacancies by dimer desorption next to an existing vacancy pair followed by readsorption, we write the creation rate as R(gain) ∼ dO2[EE] Pgain where the extra factor of Pgain denotes the probability that the just-created nearby isolated vacancies do not immediately recombine. If we define such survival as reaching a separation comparable to the typical distance between isolated vacancies, L ∼ ([E])−1/d, then it follows that46,47

(E2)

and

These results generate the evolution equations d/dt [E] = R(gain) − R(loss), which can be combined with Eq. (9) for [EE], or one can just set [EE] ≈ ε. Note that either the logarithmic correction for d = 2, or the strong correction (introducing an extra factor of [E] in both creation and annihilation terms), does not affect the steady state, so one still recovers the exact result [E]ss ≈ ε1/2. For d = 2, this treatment recovers the evolution Eq. (8) of Sec. III E apart from logarithmic corrections (which are usually ignored). However, for d = 1, the form of the equation

(E3)

is fundamentally different and produces fundamentally slower evolution, and a dimensional analysis reveals a decay time of the form τrlx ∼ ε−2/(ApO2). We do not propose that this treatment describes initial deviation from the quasi-steady state, but rather longer-time behavior. One could anticipate applying this refined treatment to the reaction model for d = 1. However, it is more tuned to capturing long-time behavior in the pure dimer adsorption-desorption model than reactive steady-state behavior.

We consider the contribution [E|E ⋅ E]|RXN to the conditional probability [E|E ⋅ E] associated with adsorption of CO on an empty site in the specified empty pair and subsequent reaction with a NN O to create an adjacent empty site. Consider the evolution of an isolated empty pair of sites from time t = 0 accounting for (i) annihilation at rate pO2/d by dimer adsorption; (ii) creation of an empty triple by CO adsorption and reaction as described above with effective rate 2 pCO Pk4, as described in Sec. IV E; and (iii) destruction of this empty triple by dimer adsorption at rate 2pO2/d. We let P2(t) denote the probability that the original vacancy dimer survives at time t, and P3(t) that the vacancy trimer has been created and survives at time t. Then, one has that

and

(F1)

We identify the lifetime of the vacancy dimer and trimer as T2 = ∫0<t<∞dt P2(t) and T3 = ∫0<t<∞dt P3(t), respectively. Integrating (F1) over time, one immediately obtains

and

(F2)

We identify the contribution [E|E ⋅ E]|RXN to [E|E ⋅ E] as

(F3)
1.
R. A.
van Santen
and
J. W.
Niemantsverdriet
,
Chemical Kinetics and Catalysis
(
Plenum
,
New York
,
1995
).
2.
K.
Reuter
, “
First-principles kinetic Monte Carlo simulation for heterogeneous catalysis
,” in
Modelling Heterogeneous Catalytic Reactions: From the Molecular Process to the Technical System
, edited by
O.
Deutschmann
(
Wiley-VCH
,
Weinheim
,
2011
).
3.
D.-J.
Liu
and
J. W.
Evans
, “
Realistic multi-site lattice–gas modeling and KMC simulation of catalytic surface reactions
,”
Prog. Surf. Sci.
88
,
393
521
(
2013
).
4.
D. P.
Landau
and
K.
Binder
,
A Guide to Monte Carlo Simulations in Statistical Physics
(
Cambridge University Press
,
Cambridge
,
2000
).
5.
R. H.
Fowler
and
E. A.
Guggenheim
,
Statistical Thermodynamics
(
Cambridge University Press
,
Cambridge
,
1956
).
6.
C.
Preston
,
Random Fields
,
Springer Lecture Notes in Mathematics
Vol.
534
(
Springer
,
Berlin
,
1976
).
7.
D.-J.
Liu
and
J. W.
Evans
,
Phys. Rev. B
70
,
193408
(
2004
).
8.
D.-J.
Liu
and
J. W.
Evans
,
J. Chem. Phys.
124
,
154705
(
2006
).
9.
D.-J.
Liu
,
J. Phys. Chem. C
111
,
14698
(
2007
).
10.
M.
Nagasaka
,
H.
Kondoh
,
I.
Nakai
, and
T.
Ohta
,
J. Chem. Phys.
126
,
044704
(
2007
).
11.
S.
Piccinin
and
M.
Stamatakis
,
ACS Catal.
4
,
2143
(
2014
).
12.
K.
Reuter
,
D.
Frenkel
, and
M.
Scheffler
,
Phys. Rev. Lett.
93
,
116105
(
2004
).
13.
K.
Reuter
and
M.
Scheffler
,
Phys. Rev. B
73
,
045443
(
2006
).
14.
R.
Rogal
,
K.
Reuter
, and
M.
Scheffler
,
Phys. Rev. Lett.
98
,
046101
(
2007
).
15.
J.
Rogal
,
K.
Reuter
, and
M.
Scheffler
,
Phys. Rev. B
77
,
155410
(
2008
).
16.
F.
Hess
,
A.
Farkas
,
A. P.
Seitsonen
, and
H.
Over
,
J. Comput. Chem.
33
,
757
(
2012
).
17.
A.
Surda
and
I.
Karasova
,
Surf. Sci.
109
,
605
(
1981
).
18.
19.
M.
Dumont
,
P.
Dufour
,
B.
Sente
, and
R. J.
Dagonnier
,
J. Catal.
122
,
95
(
1990
).
20.
J. W.
Evans
and
M. S.
Miesch
,
Surf. Sci.
245
,
401
(
1991
).
21.
M.
Tammaro
,
M.
Sabella
, and
J. W.
Evans
,
J. Chem. Phys.
103
,
10277
(
1995
).
22.
E. W.
James
,
C.
Song
, and
J. W.
Evans
,
J. Chem. Phys.
111
,
6579
(
1999
).
23.
H.
Meskine
,
S.
Matera
,
M.
Scheffler
,
K.
Reuter
, and
H.
Metiu
,
Surf. Sci.
603
,
1724
(
2009
).
24.
B.
Temel
,
H.
Meskine
,
K.
Reuter
,
M.
Scheffler
, and
H.
Metiu
,
J. Chem. Phys.
126
,
204711
(
2007
).
25.
S.
Matera
,
H.
Meskine
, and
K.
Reuter
,
J. Chem. Phys.
134
,
064713
(
2011
).
26.
M.
Barma
,
M. D.
Grynberg
, and
R. B.
Stinchcombe
,
Phys. Rev. Lett.
70
,
1033
(
1993
).
27.
M. D.
Grynberg
,
R. B.
Stinchcombe
, and
M.
Barma
,
Phys. Rev. E
47
,
4018
(
1993
).
28.
P. J.
Flory
,
J. Am. Chem. Soc.
61
,
1518
(
1939
).
29.
J. W.
Evans
,
Rev. Mod. Phys.
65
,
1281
(
1993
).
30.
R. S.
Nord
and
J. W.
Evans
,
J. Chem. Phys.
82
,
2795
(
1985
).
31.
J. W.
Evans
and
D. K.
Hoffman
,
Phys. Rev. B
30
,
2704
(
1984
).
33.

To satisfy the exact relation [E|E E] + [O|E E] + [CO|E E] = 1, one can instead set [CO|E E] = [CO] and [O|E E] = 1-[CO]-[E|E E]. This modification produces results very similar to PA2.

34.
D. M.
Ackerman
,
J.
Wang
, and
J. W.
Evans
,
Phys. Rev. Lett.
108
,
228301
(
2012
).
35.
R. M.
Ziff
,
E.
Gulari
, and
Y.
Barshad
,
Phys. Rev. Lett.
56
,
2553
(
1986
).
36.
J. W.
Evans
and
M. S.
Miesch
,
Phys. Rev. Lett.
66
,
833
(
1991
).
37.
B. J.
Brosilow
and
R. M.
Ziff
,
Phys. Rev. A
46
,
4534
(
1992
).
38.
D.-J.
Liu
,
A.
Garcia
,
J.
Wang
,
D. M.
Ackerman
,
C.-J.
Wang
, and
J. W.
Evans
, “
Kinetic Monte Carlo simulation of statistical mechanical models and coarse-grained mesoscale descriptions of catalytic reaction-diffusion processes
,”
Chem. Rev.
(in press).
39.
I.
Dornic
,
H.
Chate
,
J.
Chave
, and
H.
Hinrichsen
,
Phys. Rev. Lett.
87
,
045701
(
2001
).
40.
S.
Hong
,
A.
Karim
,
T. S.
Rahman
,
K.
Jacobi
, and
G.
Ertl
,
J. Catal.
276
,
371
(
2010
).
41.
S. I.
Shah
,
S.
Hong
, and
T. S.
Rahman
,
J. Phys. Chem. C
118
,
5226
(
2014
).
42.
F.
Hess
,
P. P. T.
Krause
,
S. F.
Rohrlack
,
J. P.
Hofmann
,
A.
Farkas
, and
H.
Over
,
Surf. Sci.
606
,
L69
(
2012
).
43.
N. G.
Van Kampen
,
Stochastic Processes in Physics and Chemistry
(
North-Holland
,
Amsterdam
,
1981
), Chap. 8.
44.
M. C.
Bartelt
and
J. W.
Evans
,
Phys. Rev. B
46
,
12675
(
1992
).
45.
B. D
Hughes
,
Random Walks in Random Environments
(
Clarendon
,
Oxford
,
1995
), Vol.
1
.
46.
D.-J.
Liu
and
J. W.
Evans
,
Phys. Rev. B
66
,
165407
(
2002
).
47.
S.
Redner
,
A Guide to First-Passage Processes
(
Academic
,
New York
,
2007
).