Transitions between strongly correlated and random steady-states for catalytic CO-oxidation on surfaces at high-pressure

Traditional mean-field (MF) rate equations of chemical kinetics¹ 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.³ In contrast, for high-P conditions,² 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.⁶ 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.³ 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 Scheffler^12,13 which was crafted to capture the essential features of CO-oxidation on the RuO₂(110) oxide surface under high-P conditions. The RuO₂(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 O₂ 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 CO₂ 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 O₂ 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.²³ 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.¹⁶ 

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 coworkers^24,25 of a transition from primarily CO-populated surface for higher CO-partial pressure, p_CO, with a random distribution of CO (or of isolated vacancies) to a near O-covered surface for lower p_CO 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)/RuO₂(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 p_CO.

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 p_CO = 0, as its spatial and dynamical features impact the reaction model for low p_CO. 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 p_CO. 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 O₂ 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)/RuO₂(110) system. The adsorption rates are chosen so that p_CO + p_O2 = 1 setting the time scale. Here, p_CO is the adsorption rate of CO per site per unit time. Also, p_O2 (p_O2/2) is the dissociative adsorption rate and d_O2 (d_O2/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 ε = d_O2/p_O2 and we will focus in our simulation study on the range from 10⁻⁶ (a typical value in experiment) up to 10⁻². The desorption rate of CO is set to d_CO = 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

and

The positive gain terms are associated with adsorption, and the negative loss terms are associated with desorption and reaction. CO adsorption occurs at rate p_CO [E], as a single empty site (which occurs with probability [E]) is required, and impingement occurs with rate p_CO. O₂ adsorption occurs at rate 2p_O2 [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 2p_O2. 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 CO₂. 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⁻⁶ is shown in Fig. 1 for both 1D and 2D models. One complication is that even these extended simulations for 10⁶ MCS still cannot access the true steady-state for very low p_CO (below about 10⁻⁵ in 1D and 10⁻⁴ 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 p_CO, 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 d_O2. Thus, Fig. 2 shows steady-state behavior for the 1D model for d_O2 = 10⁻⁴, 10⁻³, and 10⁻² where simulation can more readily access the true steady-state. Trends in behavior for the 2D model with increasing d_O2 (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.

Regime I for high p_CO ≈ 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 p_CO 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 p_CO 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 d_O2 of about 10⁻³ or above where they are compressed and shifted to higher p_CO. See Fig. 2.

For a comprehensive understanding of the reaction model of Sec. II for low p_CO, it is invaluable to develop a detailed characterization of limiting behavior of the pure oxygen adsorption-desorption limit with p_CO = 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 p_CO = 0 (so p_O2 = 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.

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.²⁷ One simply interchanges occupied and empty sites on only the + sublattice. Then, one obtains a pure particle diffusion model where particles hop with rate p_O2/d from + to NN − sites, and with rate d_O2/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 L^d − M on a lattice with L^d 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]_ss² and [OO]_ss = [O]_ss², 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

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 d_O2 = 10⁻⁴ 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.

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 p_CO.²⁵) For dimension d = 1 or 2, this implies that

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

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

and

This result consistent with our simulation analysis for ε = 10⁻⁴ 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 by^28,30

so that

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.²⁷ For such a pure diffusion model, one expects diffusive decay of spatial correlations,³¹ and perhaps also in site occupancy kinetics.

For the special case where p_O2 = d_O2, 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.³² 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 ∼t^1/2, corresponding to a volume ∼t^d/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 p_O2≠d_O2, including the case of most relevance where ε = d_O2/p_O2 ≪ 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 p_O2≠d_O2.²⁷) 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 d_O2 = 10⁻⁴ 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 d_O2 = 10⁻². Results shown in Fig. 5 are consistent with the proposal that [E] − [E]_ss ∼ t^−d/2, as t → ∞.

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

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⁻² ≈ 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(ε⁻¹/p_O2) = O(d_O2⁻¹) 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 O₂ desorption event occurs at rate d_O2 adjacent to an existing empty pair on the surface thereby creating a connected empty linear or bent quartet and (ii) subsequent (fast) adsorption of O₂ 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 R_gain = Ad_O2 [E ⋅ E]^∗. Isolated empty sites effectively diffuse when (i) a slow O₂ desorption event at rate d_O2 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 O₂ to annihilate this pair. Thus, the destruction rate has the form R_loss = Bd_O2([E]^∗)². We conclude that

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.

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

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]^∗ ≈ d_O2/p_O2 = ε. Then, the solution to (8) with A = B becomes

where

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/(Ap_O2), but not by τ_MF = ε^−1/2/(4p_O2) or τ_pair = O(ε⁻¹/p_O2). 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(ε⁻²/p_O2). 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 p_CO. Regime I for high p_CO ≈ 1 corresponds to a primarily CO-populated surface. Regime II for lower p_CO below about 0.1 corresponds to a near O-covered surface with strong spatial correlations, and Regime III for very low p_CO 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” p_CO ≈ 1 and p_O2 ≪ 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 ≈ p_CO/(p_CO + d_CO) ≈ 0.909 and that [E]_ss ≈ d_CO/(p_CO + d_CO) ≈ 0.091. With p_O2 ≪ 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 p_O2/d times [E ⋅ E]_ss ≈ ([E]_ss)² times the number (2d) of empty pairs overlapping the site and (ii) the probability, P_k1 ≈ (2d − 1)k/[(2d − 1)k + p_O2] ≈ 1 (for d_O2 ≪ k) that this O reacts with a NN⋅CO rather than undergoing recombinative desorption. Thus, one has that

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

With regard to the behavior for a broader range of p_CO 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 p_CO 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 p_CO, and the probability that such a CO reacts with adjacent O rather than desorbing is P_k2 = 2dk/(d_CO + 2dk), it follows that

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

Here, we consider Regime II with p_CO below about 0.1 where the surface becomes mainly O-covered. A key observation by Matera, Meskine, and K. Reuter (MMR)²⁵ for the analogous regime in their model of CO-oxidation on RuO₂(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 O₂ 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 p_CO = 0). Indeed, our heuristic rate equation treatment in Sec. IV F based upon this picture does imply crossover with increasing p_CO 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 d_O2 = 10⁻⁴ 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 ∼ 10⁶ MCS (10⁵ 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 p_CO, and only approach the MF value of [E] for very low p_CO 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.

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]². For the standard PA, we start with exact rate equations for pair quantities. For example, for d = 1, one has

where the gain terms correspond to O₂ 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 O₂ 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 d_O2 = 10⁻⁴ (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.

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]²/[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

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 p_CO. See Figs. 8 and 9.

Finally, we also describe results from the approximation implemented by MMR²⁵ 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]², 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 p_O2. The improvement is not so significant for d = 2, although we note that the effectiveness of the MMR generally improves for lower d_O2.

It is clear from our above analysis of various truncation approximations, and from the earlier work of Matera et al.,²⁵ 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 p_CO, and the latter occurs with probability P_k4 = k/(k + d_CO) for d = 1, and P_k4 = k/(3k + d_CO) for d = 2. Also, since the lifetime of an empty pair is given by τ_2E = d/p_O2, the fraction of time for which the additional site is empty is given by [E|E ⋅ E]|_RXN = p_COP_k4τ_2E. See Appendix  F. Thus, we have that

Fig. 10 shows the high accuracy of this formulation compared with the MF and pair approximations for d_O2 = 10⁻⁴ 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.³³ The results of PA2 shown in Figs. 8 and 9 significantly improve upon all other approximations described previously (MF, PA, PA1, and MMR).

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 O₂ to fill this empty pair. The rate of CO-adsorption onto an isolated empty site is p_CO, and the probability that this CO reacts with an adjacent O rather than first desorbing is given by P_k2 = 2dk/(d_CO + 2dk). Thus, the reactive loss rate has the form R_loss(rxn) = p_COP_k2[E]^∗. Isolated vacancies can also be created via reaction from vacancy pairs: a CO adsorbs at rate p_CO, reacts before desorbing with probability P_k3 = (2d − 1)k/[d_CO + (2d − 1)k] with a NN O to create a vacancy triple, and subsequent O₂ adsorption produces an isolated vacancy. Thus, the reactive gain rate has the form R_gain(rxn) = 2p_COP_k3[EE]^∗. Choosing prefactors A = B and augmenting (8) by these terms yield a heuristic rate equation for [E]^∗ having the form

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

To determine the corresponding TOF, we use TOF ≈ p_COP_k[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, p_CO = p_CO⁻, of the crossover regime will correspond to the onset of significant deviation of [E]^∗_ss from ε^1/2 when (p_COP_k2)/(2Ap_O2ε) = O(ε^1/2). The upper boundary, p_CO = p_CO+, of the crossover regime is less precisely determined. One criterion to determine p_CO+ is to require that [E]^∗_ss be reduced from ε^1/2 to O(ε) implying that [E]^∗_ss ≈ (Ap_O2ε²)/(p_COP_k2) = O(ε). In conclusion, these analyses suggest that the lower and upper boundaries of the crossover region satisfy

and

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

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.³⁴ 

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 RuO₂(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, p_CO, to a strongly correlated near-O-covered steady-state for low p_CO 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 p_CO. The subtle and slow diffusive kinetics for the pure oxygen adsorption-desorption regime at p_CO = 0 impacts behavior for low p_CO. 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 p_CO as these are kinetically rather than thermodynamically induced. The novel treatment of Matera et al.²⁵ 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 p_CO + p_O2 = 1, d_CO = 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 d_O2 of 10⁻⁴ 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 d_CO 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 = 2³⁵ which includes the same ingredients as our reaction model except that d_CO = d_O2 = 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 d_O2 > 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 k³⁶ and also for d_CO > 0 below a finite (small) critical value.³⁷ Thus, we expect a discontinuous CO-poisoning transition in our current reaction model for d = 2 with very low or zero d_O2 provided that one sufficiently reduces d_CO.³⁸ 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 2p_CO or with O₂ at rate p_O2/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.³⁹ 

Finally, we remark that the general strategies developed here to characterize kinetically induced spatial correlations should extend beyond CO-oxidation on RuO₂(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 HCl⁴² on RuO₂(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 RuO₂(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.

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 p_O2/d, and from − to NN + sublattice sites at rate d_O2/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

and similarly for [X₋]. From these equations, it is immediately clear that [X₊] + [X₋] is constant. For the special case where p_O2 = d_O2, 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₋]₁ = [X₊X₋]), are similarly defined. In addition, we define [X₋E₊ − X₊]_n = [X₊ − E₊X₋]_n as the probability for an X₋E₊ pair with E₊ separated by odd n ≥ 1 lattice constants from an X₊, etc. Then, for even n ≥ 2, one has that

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 p_O2 = d_O2, 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 form^31,43

and

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

Of the resulting branches of eigenvalues, μ_j(q), with j = 0-2, one branch j = 0 satisfies μ₀ = 0 when q = 0, and more generally μ₀(q) ∼ − Aq² for small q with c > 0. This branch of eigenfunctions with μ₀(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

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

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

Of the 3 branches of eigenvalues, μ_j(q_x, q_y), with j = 0, 1, 2, one Goldstone branch j = 0 satisfies μ₀(q) ∼ − B[(q_x)² + (q_y)²] 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

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

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

and

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]²/[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)² and [O ⋅ E]_ss = [O]_ss[E]_ss. Forming a $2×1 vector δ ¯$ with entries δ₁ = δ[O ⋅ O] and δ₂ = δ[O ⋅ E], the linearized evolution equations have the form

The entries of $J =$ for the 1D problem are

An analysis of the eigenvalues λ_i of $J =$ reveals that λ₁ = − 5p_O2 + O(p_O2ε^1/2) and negative λ₂ = O(p_O2ε). Thus, the key observation is that terms O(ε^1/2) are absent in λ₂ due to cancellation. Thus, asymptotic decay to the steady-state which is controlled by λ₂ occurs fundamentally slower than in the MF approximation with a characteristic time τ_pair = O(ε⁻¹/p_O2) = O(d_O2⁻¹).

Exactly, the same scenario applies for the pair approximation analysis for the 2D problem with τ_pair = O(ε⁻¹/p_O2) = O(d_O2⁻¹), but now

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

The standard pair approximation applies [E ⋅ E ⋅ E] ≈ [E ⋅ E]²/[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.⁴⁴ τ is estimated as the lifetime of a random walker amongst a random distribution of trap sites with density [E]^∗. Thus, one obtains^44,45

but

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) ∼ d_O2[EE]^∗ P_gain where the extra factor of P_gain 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 that^46,47

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

is fundamentally different and produces fundamentally slower evolution, and a dimensional analysis reveals a decay time of the form τ_rlx^∗ ∼ ε⁻²/(Ap_O2). 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 p_O2/d by dimer adsorption; (ii) creation of an empty triple by CO adsorption and reaction as described above with effective rate 2 p_CO P_k4, as described in Sec. IV E; and (iii) destruction of this empty triple by dimer adsorption at rate 2p_O2/d. We let P₂(t) denote the probability that the original vacancy dimer survives at time t, and P₃(t) that the vacancy trimer has been created and survives at time t. Then, one has that

and

We identify the lifetime of the vacancy dimer and trimer as T₂ = ∫_0<t<∞dt P₂(t) and T₃ = ∫_0<t<∞dt P₃(t), respectively. Integrating (F1) over time, one immediately obtains

and

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