Beyond mean-field approximations for accurate and computationally efficient models of on-lattice chemical kinetics

Modeling the kinetics of surface catalyzed reactions is essential for the design of reactors and chemical processes. The majority of microkinetic models employ mean-ﬁeld approximations, which lead to an approximate description of catalytic kinetics by assuming spatially uncorrelated adsorbates. On the other hand, kinetic Monte Carlo (KMC) methods provide a discrete-space continuous-time stochastic formulation that enables an accurate treatment of spatial correlations in the adlayer, but at a signiﬁcant computation cost. In this work, we use the so-called cluster mean-ﬁeld approach to develop higher order approximations that systematically increase the accuracy of kinetic models by treating spatial correlations at a progressively higher level of detail. We further demonstrate our approach on a reduced model for NO oxidation incorporating ﬁrst nearest-neighbor lateral interactions and construct a sequence of approximations of increasingly higher accuracy, which we compare with KMC and mean-ﬁeld. The latter is found to perform rather poorly, overestimating the turnover frequency by several orders of magnitude for this system. On the other hand, our approximations, while more computationally intense than the traditional mean-ﬁeld treatment, still achieve tremendous computational savings compared to KMC simulations, thereby opening the way for employing them in multiscale modeling frameworks. ' 2017 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). [http://dx.doi.org/10.1063/1.4991690]


I. INTRODUCTION
Heterogeneously catalyzed reactions are widely used in the chemical industry, but also in everyday applications. Examples of such applications range from petroleum refining to automotive emission control. 1 In this type of catalysis, the reactants adsorb onto the catalyst surface, via the formation of chemical or physical bonds. After surface reaction, the product desorbs from the surface and diffuses away. 2,3 Clearly, the presence of the catalyst provides lower energy pathways or alternative elementary steps to get the desired chemical product from that in its absence. That is why there is a continued interest in the chemical industry and in academia to develop more active, selective, stable, and less expensive catalysts. 4 However, the development of such novel catalysts is not an easy task. A central aspect of this endeavor is to understand the elementary reaction steps and model the dynamics of catalytic processes. To this end, kinetic modeling approaches, such as microkinetic mean-field (MKM) models or the so-called kinetic Monte Carlo (KMC) simulation, are of paramount importance.
MKM models of chemical kinetics as well as KMC simulations have indeed proved very useful in studying elementary processes occurring on reactive surfaces. 5-7 From a fundamental viewpoint, both approaches originate from the same Markovian master equation (MME) for the time evolution of a) m.stamatakis@ucl.ac.uk the catalytic system. 8 The MME is derived from first principles and its transition probabilities or rate constants are computed with quantum chemical methods. 7,9 General MKM models result by applying a system-size expansion, thereby focusing on the behavior of the master equation at the thermodynamic limit of very large lattices. In this limit, stochastic fluctuations become negligible and one can formulate ordinary differential equations (ODEs) that describe the temporal evolution of the average surface coverage. 10 Instead, KMC simulations provide trajectories (stochastic realizations) whose statistics follow the master equation. 11,12 Although traditional MKM models are highly efficient, they may lead to inaccurate predictions because they ignore details about the spatial correlations in the adlayer. Such correlations can arise from slow diffusion in tandem with reaction or from adsorbate-adsorbate lateral interactions. [13][14][15] An example of the former case is CO oxidation on RuO 2 (110). 13 In this situation, a microkinetic model incorporating the correct expressions of the pair probabilities of finding vacancies or adatoms in neighboring sites is enough to yield results in quantitative agreement with KMC. 14,15 A different route is to numerically approximate the solution of the highdimensional MME directly. The tensor train approximation recently reported by Gelβ et al., 16 for solving a MME for the aforementioned CO oxidation process, is an interesting example of this approach. However, as noted by Gelβ et al., it is not-trivial to develop such approximations for other surface reaction models.
Treating lateral interactions is also non-trivial. There have been some attempts to improve the accuracy of microkinetic models to better approximate the master equation. 17,18 These approaches generalize traditional MKM methods by introducing an infinite set of evolution equations for spatial correlations, which is truncated at some level using moment closure techniques. The main trade-off between these moment closures is between accuracy and simplicity. KMC simulations, on the other hand, explicitly treat all these correlations but are computationally expensive. 19 Thus, more efforts to find new general methodologies that can capture kinetics of surface catalyzed reactions in a computationally efficient way without sacrificing accuracy are needed.
In this article, we focus on the development of a hierarchy of approximations, based on the so-called cluster mean-field approximation for kinetic lattice-gas (LG) models, 20,21 that allow us to calculate catalytic rates or turnover frequencies (TOFs) in an accurate and efficient way. Our conceptually simpler approach can capture adlayer inhomogeneity due to lateral interactions between adsorbates together with the coverage dependence of the activation barrier of reaction events. Instead of deriving and truncating an infinite system of equations for the correlations, we couple the dynamics of an explicitly treated cluster of lattice sites with a surrounding mean-field. We further assess the error of our approximations in a reduced model for NO oxidation and NO 2 reduction reaction on Pt(111), a realistic chemistry of practical importance. We adopt an analogue of the minimal kinetic model of Wu et al. 22 in which the NO oxidation/NO 2 reduction process is rapid, the O 2 dissociative adsorption/oxygen associative desorption process is rate-limiting, and oxygen is the most abundant surface intermediate. In agreement with the low barriers for oxygen diffusion on Pt(111), 23 we assume that the oxygen atoms diffuse very fast on the surface.
The paper is organized as follows. In Sec. II, we present the Hamiltonian formulation framework together with the description of our cluster mean-field approach. Subsequently, in the same section, we validate and benchmark our approaches by calculating relevant kinetic quantities in the reduced model for NO oxidation and NO 2 reduction on Pt(111). In Sec. III, we compare the accuracy of our approximations under different operating conditions. Finally, summary and conclusions are presented in Sec. IV.

II. METHODOLOGY
In this section, we present the LG Hamiltonian first, which is an exact formalism for treating correlations of adparticles on a lattice, and then discuss a set of approximations which treat explicitly a single cluster of sites coupled with a mean-field. We subsequently benchmark these approaches by calculating relevant kinetic quantities in a reduced model for NO oxidation and NO 2 reduction on Pt(111).

Lattice-gas Hamiltonian formulation
For the case of fast diffusion of adsorbates, we are allowed to assume quasi-equilibrated adlayer. Hence, we can compute quantities of interest as averages, using standard statistical mechanical approaches. To this end, we need a Hamiltonian, and an appropriate choice is the LG, since we are working on discrete lattices. 24 This approach is a common tool to accurately describe spatial correlations and ordering in chemisorbed layers due to adspecies interactions. Let us define a regular lattice formed by N L sites. We number all of the sites on the surface by a single index i and define an occupancy variable σ i , where σ i = 1, if the site is occupied by an adsorbate, and σ i = 0, if the site is empty. Then, the LG Hamiltonian of any state σ = {σ 1 , σ 2 , . . . , σ N L } is expanded in polynomial "clusters" or "figures" of local variables σ i , where N L is the number of lattice sites, h 0 is the free energy of the surface in the absence of the adsorbate, and h 1 is the free energy of one adsorbate. Interactions between adsorbates are captured in pairwise (J ij ) and higher-order terms. The summations run over all sites of the lattice and the factor that multiplies the third term takes care of counting the same pattern multiple times. Note that for N L sites on the lattice, the total number of possible adsorbate configurations is 2 N L . Assuming that the adsorption and desorption events are Markovian, a function P σ can be introduced, which gives the probability that a given microscopic configuration σ is realized at time t. Then, the temporal evolution of the system is determined by the so-called Markovian master equation or MME, where σ and β are configurations of the adlayer, P σ and P β are their probabilities, and k βσ is the rate constant or average escape rates from configuration σ to configuration β in units of inverse time. 8 Because the adlayer is in quasi-equilibrium, it is clear that with the grand canonical partition function given by where N σ is the number of adsorbates of state σ, µ is the chemical potential, k B is the Boltzmann constant, and T is temperature. The average number of adsorbates on the lattice is given by In kinetics, a relevant quantity is where N β is the number of adsorbates of state β and the sum is over all possible configurations σ and β. The double sum of the right hand side can be broken into terms, each of which pertains to one of the possible elementary events (e.g., dissociative O 2 adsorption) happening at perhaps different "environments," i.e., different configurations of spectators.
When averaging such terms, (N σ −N β ) will always be the same and equal to a stoichiometric coefficient (e.g., 2 for dissociative O 2 adsorption), and the remaining factors (k σ β P β ) will give rise to an average rate constant at which this process evolves. Note that if configuration β cannot give rise to configuration σ via a single elementary event, the corresponding k σ β term is zero.

Cluster mean-field approximation
We now proceed to introduce a type of multi-effectivefield cluster mean-field approximation 20,21 that allows us to calculate accurate surface coverages and TOFs with minimal cost. The main idea is to obtain these quantities on a small cluster Ω on the lattice containing N Ω L < N L sites, handling the neighborhood as a "cloud" of particles that interact with the explicitly treated particles within the cluster. When the degrees of freedom outside of the cluster are eliminated, many effective fields acting on single adsorbates as well as on multiplets of adsorbates appear in the cluster. The Hamiltonian of the small clusters under the "cloud" of particles (H C (σ)) is defined as where with h 0,Ω = h 0 N Ω L /N L being only a reference energy that does not change any of the computed properties (e.g., coverages and correlations). The new terms in Eq. (7) contain the effective fields or corrections acting on single adsorbates and multisite clusters of adsorbates. Such corrections are not applied to the adsorption energy at the central site (single body term) and any pairwise (or higher) interaction containing the central site. Note that from now on σ = {σ 1 , σ 2 , . . . , σ N Ω L }. The constants a i and b ij are the one-body and two-body effective fields, respectively. These effective fields are determined self-consistently by imposing (i) the equality of the coverages of all sites and (ii) the equality of correlation functions for symmetrically equivalent clusters of sites (e.g., pairs and triplets). In formulating the appropriate consistency equations, the symmetry operations (translations and rotations) of the original lattice are applied. The correlations are obtained where S is a product of the state variables of the individual sites whose correlation function is sought, and the summation is taken over all states of the cluster, with P σ being the probability of the corresponding configuration of the cluster. Thus, S = σ 1 σ 2 would yield the 1NN correlation, whereas S = σ 1 would simply give the coverage. In Subsection II C, we analyze a set of hierarchical approximations based on the above formalism.
In the hierarchy of approximations that we develop, the smallest such cluster is a single site (which can be vacant or occupied by an adsorbate), surrounded by the implicitly treated "cloud" of adsorbates (see the left panel of Fig. 1). This is identical to the well-known Bragg-Williams or mean-field (MF) approximation. 25 The second level in the hierarchy of approximations is the so-called Bethe-Peierls (BP) approximation in which the central site and its 1NN sites are treated explicitly (see the middle panel of Fig. 1). 26 This approximation entails keeping track of the state of each of the seven sites of the cluster and solving relevant quantities for the 2 N Ω L configurations of the cluster, with everything beyond the 1NN sites treated as a "cloud" of adsorbates that interact with the cluster. Any approximations on clusters bigger than the 7-site cluster of Fig. 1 will be referred to as Kikuchi approximations. In particular, we treat explicitly sites up to second (K2NN) and third (K3NN) nearest-neighbors of the central site (see the right panel of Fig. 1 for the case of 2NN sites). This last approach is inspired by the extensively reformulated Kikuchi cluster-variation method (CVM). 26 Note that in Eq. (7) subscript C will be replaced by the name of the corresponding approximation.
There are variations of each of these approximations (MF, BP, K2NN, and K3NN), each with a different Hamiltonian [Eq. (7)] and corresponding effective field parameters. After solving for these parameters using a system of equations of the form of Eq. (9), we have a Hamiltonian that can be used to obtain average coverages [Eq. (5)] or correlations [Eq. (9)], but also average TOFs from Eq. (6). We will demonstrate and benchmark these approximations using the NO oxidation/NO 2 reduction system, but first let us briefly introduce this system in Sec. II B.

B. NO oxidation and NO 2 reduction on Pt(111)
NO oxidation to NO 2 is relevant to emission control technologies and platinum is a common catalyst used in these applications. 27,28 NO oxidation rates have been measured on single-crystal Pt surfaces as well as on supported Pt nanoparticles. 29,30 However, the NO oxidation mechanism and its transferability from surfaces to nanoparticles remain unresolved. 31 It is known that the greatest TOFs of this reaction occur over the close-packed Pt(111) facets of supported particles. 29,30 Moreover, it has been recognized that, during reaction conditions, oxygen (denoted as O) dominates the Pt surface, and high O coverages have been observed to account for NO oxidation activity. 32 Due to short-range repulsive interactions, the adsorbed oxygen atoms exhibit superlattice ordering, with their primary location given by threefold face-centered-cubic (fcc) sites of the Pt(111) hexagonal surface. 33,34 Many studies also show that the barrier of O 2 dissociative adsorption on Pt(111) is strongly affected by the presence of co-adsorbed oxygen atoms, 32,35 while conversion of NO to NO 2 has an activation barrier that is rather insensitive to the presence of adsorbed oxygen atoms. 32 Therefore, an accurate description of NO oxidation must consider both realistic surface O coverage and its influence on the energetics and kinetics of the catalytic system.
The first KMC simulations of NO oxidation on Pt(111) were reported some years ago. 36,37 The so-called cluster expansion (CE) Hamiltonian techniques were implemented to improve the description of O repulsive interactions on Pt(111). 38 Assuming that O 2 dissociation is rate-limiting, Wu et al. 22 combined CE Hamiltonians incorporating longrange interactions and many-body terms with equilibrium Monte Carlo simulations for oxygen on Pt(111) and Brønsted-Evans-Polanyi (BEP) relations for O 2 dissociation to estimate catalytic NO oxidation rates. More recently, detailed KMC simulations of the full reaction mechanism based on a lattice-gas (LG) Hamiltonian formulation have also been reported. 19 However, Monte Carlo simulations are computationally expensive, motivating the development of more efficient approaches.
In our analysis, the catalytic reaction proceeds via the following steps: 19,22,31,32 (10) (11) (12) with * and O * denoting a vacant adsorption site and adsorbed oxygen atoms, respectively. The last reaction denotes diffusion of O * between neighboring adsorption sites. We adopt the assumptions by Wu et al. 22,32 in which reaction step (10) is rapid and reaction step (11) is the rate-limiting process. We assume that, in agreement with the low barriers for O * diffusion on Pt(111), 23 the oxygen atoms diffuse very fast on the surface. This ensures oxygen adlayer equilibration and the applicability of the lattice-gas Hamiltonian formulation described above for the calculation of average properties. Moreover, we assume an equilibrated interconversion between NO and NO 2 which allows us to relate the chemical potential of surface oxygen to that of the gas phase species as where G o X is the standard state Gibbs free energies for species X at 1 bar, T is temperature, and k B is the Boltzmann constant. The Gibbs free energies are calculated from the NIST Chemistry WebBook. 39 From all previous assumptions, one can deduce that oxygen coverage on the surface is set by the pressure ratio of NO 2 to NO. It is also straightforward to conclude that because the O 2 dissociative adsorption/O * associative desorption step is the rate limiting process, the NO oxidation rate or TOF is determined by the overall reaction rate of chemical step (11). The assumptions behind Eq. (13) are based on kinetic observations by Getman et al. 32 In particular, their experimental results show that while the rates of O 2 and NO 2 dissociations are comparable at low oxygen coverage, NO 2 dissociation is many orders of magnitude faster at O coverages typical of NO oxidation catalysis. This indicates that O 2 dissociation can actually be considered as a rate-limiting process under these conditions, and the chemical potential of oxygen on the surface can be equated to the difference of the chemical potentials of NO and NO 2 . Such an assumption, combined with a first-principles-based cluster expansion, was subsequently successful to reproduce observed O coverages and recovered NO oxidation rates and reaction orders in good agreement with experiments (see, for example, the work of Wu et al. 22 ). However, we anticipate that Eq. (13) may not necessarily hold at a very low coverage, where it is believed that the oxidation of NO is not fast any more. 40 In Sec. II C, we apply the methodology described in Sec. II A to calculate the total NO oxidation rate from equilibrium configurations of O * on a Pt(111) surface at catalytically relevant ranges of µ O * and temperature T.

C. Cluster mean-field approach to NO oxidation and NO 2 reduction
Let us define a hexagonal lattice formed by the N L fcc sites of a Pt(111) surface. On this lattice, the average stationary overall reaction rate or TOF is the superposition of the corresponding overall reaction rates of the elementary reaction steps. In our case, these elementary reaction steps are O 2 dissociative adsorption and O * associative desorption. Therefore, one can define the average overall reaction rate as where R ads and R des denote the average overall adsorption and desorption rates, respectively. 7,41,42 From Eq. (6), we can get the intensive TOF (normalized per site) as where in the summation k σ β = k ads σ β > 0 if σ can originate from β by the adsorption of two oxygen atoms, otherwise k σ β = 0. Similarly, for oxygen associative desorption, one has that where, as before k σ β = k des σ β > 0 when, in the summation, configuration σ can originate from configuration β by the desorption of two oxygen atoms from the surface.
The local rate constants of O 2 dissociative adsorption in configuration β are necessarily influenced by interactions between or with adsorbed oxygen species, and thus, they are distinguished by the local environment. According to transition state theory, the local rate of O 2 dissociative adsorption is given by where A is the pre-exponential factor described in the supplementary material. Using the Hamiltonian formulation [see Eq.
(1) or (7)], one can calculate the overall energy changes (reaction energies) due to candidate O 2 dissociative adsorption events, from the differences in surface energies of initial and final surface states that differ by the addition of adjacent oxygen atoms into previously vacant sites, 19 where ∆E O 2 is the change in energy of oxygen gas molecules. These dissociative adsorption energies are mapped to the local activation energies for O 2 dissociation via BEP-type relations as where E ‡ a,0 and ∆E rxn,0 are the activation energy and the reaction energy at the zero coverage limit, having the values of 0.020 eV and 2h 1 + J = 2.100 eV, respectively. Density functional theory (DFT) calculations reveal that the proximity factor ω is equal to unity 22 for O 2 dissociation. Note that BEP equation (19) thus parameterised is mathematically identical to the BEP developed by Wu et al. 22,32 and is used consistently for each of the approximate models developed in our work.
Desorption of oxygen atoms from two occupied sites of the lattice represents the reverse process to dissociative oxygen adsorption, and as such, it has to satisfy the detailed balance condition. Thus, the expression for the rate constant of O 2 desorption is with A being a pre-exponential factor that guarantees a detailed balance. This factor is described in the supplementary material of this work. Similar to the adsorption case, we capture the desorption barrier via a BEP as with where E ‡ d (β) and E ‡ d,0 are the activation energies for desorption in configuration β and the surface configuration which contains two neighboring adsorbates but without spectators, respectively.
To capture the oxygen adlayer energetics on Pt(111), Schneider and co-workers used an Ising Hamiltonian formalism together with a density functional theory (DFT) database of 66 Pt(111)/O formation energies and developed four cluster expansion Hamiltonians of increasing accuracy by including 3, 5, 8, and 12 figures or clusters. 38 In a subsequent work, Nielsen et al. adopted a lattice-gas Hamiltonian formulation and obtained the lattice-gas parameters by mapping those of the Ising model. 19 In this work, we use such lattice-gas parameters (more details of the mapping can be found in Ref. 19). However, as a benchmark, we assume only two-body 1NN interactions between adsorbed oxygen atoms (pairwise interactions) and J ij = J ∀ i, j. In other words, a cluster lattice-gas Hamiltonian containing only 3 figures will be implemented. 19,22,38 Then, we have h 0 = 0.027 eV, h 1 = 1.200 eV + corrections (zero-point and vibrational energies), J = 0.300 eV, and E ‡ a,0 = 0.020 eV. 19 The restriction to 1NN interactions is indeed an idealisation. For example, Nielsen et al. 19 compared cluster expansions of increasing fidelity (fitted to the same DFT data) and demonstrated that one needs to include up to 3NN interactions for an accurate prediction of the TOF. However, in this work, we start from a simple case to test the performance of our approximations in the first instance and obtain a fundamental understanding of their behaviour.
In Subsections II C 1-II C 3 and III, we evaluate the performance of the cluster approximations introduced in Sec. II A 2 by calculating coverages and TOFs from the lattice gas Hamiltonian model for NO oxidation and NO 2 reduction on Pt(111). Our theoretical results are compared with Metropolis Monte Carlo (MMC) simulations for oxygen adlayers [43][44][45] and KMC simulations of the full reaction model (more details about the Monte Carlo simulations in the supplementary material). The self-consistency conditions necessary to get the effective fields or corrections imposed to the cluster Hamiltonian are obtained in Matlab using the built-in function fsolve, which solves systems of equations employing the trust-region dogleg algorithm. 46

Mean-field approximation
In this approximation, the treatment of the whole lattice is restricted to a single specified site whose occupancy is denoted as σ 1 (the smallest possible cluster). This site can be occupied or unoccupied by an oxygen atom. We approximately include its interaction with the rest of the lattice into an effective field representing a "cloud" of oxygen atoms. Then, from Eq. (7) with two-body 1NN interactions, we get the cluster Hamiltonian given by with The expectation value σ j of the sites around σ 1 is assumed to be equal to the oxygen coverage θ O on the surface, z is the coordination number of the site, and J is the two-body or pairwise interaction constant. Thus, the self-consistency condition is reduced to The reaction energies for O 2 dissociative adsorption with surface spectators are given by This means that the adsorption and desorption rate constants are only functions of the overall coverage on the surface or θ O . However, they are not configuration-dependent because the notion of locality does not exist in this approximation. Thus, adsorption and desorption rates are given by and  (19) and (21)]. After substituting Eqs. (15) and (16) and carrying out the necessary algebraic manipulations, we obtain a TOF given by 7,10 In Sec. III, we compare this TOF with those obtained from the next two approximations and the corresponding KMC simulations of the full system.

Bethe-Peierls approximation
The MF approximation ignores the fact that the probability of one site being occupied is actually dependent on the probability of an adjacent site being occupied. A natural correction to the MF approximation is the well-known BP approximation. In its original formulation, the idea behind the BP approach is to consider a cluster of a few atoms interacting with its environment in a mean-field sense. To apply this idea to the equilibrated adlayer of chemisorbed oxygen on our hexagonal lattice, we treat explicitly a central site and its z = 6 first nearest-neighbors with everything beyond the cluster treated as a "cloud" of atoms interacting with the explicit atoms of the cluster (a cluster of z + 1 sites rather than just onesite cluster, as in the MF approximation). The main advantage of this approximation lies in the fact that, to evaluate Eqs. (15) and (16), we only need to take into account the 2 7 possible states of the cluster (of course symmetry arguments can be used to reduce this number of states and improve computational efficiency). From Eq. (7), the cluster Hamiltonian with two-body 1NN interactions for this approximation is given by This Hamiltonian contains only one-body corrections mediated by a single effective field (note that we did assume that a i = a 1 ∀ i, where i runs from 2 to 7). As shown in Table I, one can still explore several levels of approximations within a single cluster. In the next level, which we call BPE approximation, TABLE I. Different levels of approximations within the Bethe-Peierls approach. The middle column presents the patterns corresponding to the effective fields of Eq. (7). Arabic numerals correspond to single-body terms on the sites noted (and all equivalent sites). Latin numerals correspond to the patterns shown in Fig. 2. The last column shows the self-consistency equations to be solved for the effective field parameters. Note that an approximation in a subsequent row in the table also employs the consistency corrections of all approximations in previous rows; thus the BPEC approximation uses both self-consistency equations noted.

Approximations
Patterns Consistency conditions lateral interactions between adatoms bound to edge sites are considered. In this case, the cluster Hamiltonian is expressed as Finally, in the BPEC approximation, effective fields to edge atom-atom interactions are added. These effective terms will allow us to correct the 1NN pairwise correlation functions so that we achieve self-consistency. The cluster Hamiltonian, when these corrections are taken into account, is written as where b 1 is the new effective field. [Note also that in Eq. (7) we consider b ij = b 1 ∀ i, j.] For the BP and BPE cases, one only needs to solve a single self-consistency equation, while for the BPEC approximation, an extra equation is needed to obtain b 1 [see Figs. 2(a) and 2(b) for the cluster and the energetic patterns used to solve those self-consistency calculations]. Once these effective fields have been determined, any relevant quantity of the system, for example, the TOF from Eqs. (15) and (16) or the surface coverage from Eq. (5), can be calculated by properly weighted averaging based on these Hamiltonians. The choice of which correlations to correct is somewhat arbitrary. For instance, we could assume that the 1NN σ 1 σ 2 correlation is the "most accurate" one and introduce a correction for σ 2 σ 3 , manifested as an added parameter to the corresponding interaction term (J + b 1 )σ 2 σ 3 (with appropriate additional terms for σ 3 σ 4 , σ 4 σ 5 , etc., due to rotational symmetry).
On the other hand, we could correct σ 1 σ 2 , σ 1 σ 3 , etc., using σ 2 σ 3 as the reference. We choose to keep the correlations of the central site as the reference (the "correct" ones) and impose correction parameters to the 1NN correlations among edge sites. The self consistency equations will then have the form appearing in Table I.

Kikuchi approximations
The original version of the BP approximation does not take into account long range correlations. In this part, in order to increase the accuracy of our calculations, we include clusters of larger sizes. As mentioned above, we call this method the Kikuchi approximation because of its conceptual similarity with the traditional Kikuchi CVM. In contrast to the CVM, which aims at finding an approximate expression of the entropy of the interacting system of particles, 20 the main idea of our approach is to come up with an approximate Hamiltonian that allows us to average any quantity (such as the TOF of a surface reaction). Moreover, given that it is recognized that the use of large clusters alone does not necessarily improve the predictions of the BP approximation (even when there are only 1NN interactions as in our case 26 ), in this work, we include 2NN effective fields through the so-called "ghost interactions." 47 We start with a cluster containing up to 2NN of the central site. We denote this approximation as K2NN. The case with corrections up to 1NN correlations or K2NNC1 approximation leads to a cluster characterized by the Hamiltonian, where, as Table II shows, the single body and two-body effective fields a 1 , a 2 , b 1 , and b 2 are obtained from four selfconsistency conditions [see Figs. 3(c) and 3(d) for the patterns used to obtain those values]. Note that the correlations of patterns (i) and (vi) are assumed to be the "correct" ones, so no such effective fields appear in Table II.
In order to increase the accuracy of our approach, we correct correlations with range longer than 1NN, by considering "ghost interaction" terms up to two-body 2NN, added as extra terms to the Hamiltonian. Those extra terms contain effective fields that affect the correlations between two adatoms with 2NN separation. We denote this case as the K2NNC2 approximation for which the Hamiltonian is given by  (7). Arabic numerals correspond to single-body terms on the sites noted (and all equivalent sites). Latin numerals correspond to the patterns shown in Fig. 3. The last column shows the self-consistency equations to be solved for the effective field parameters. Note that an approximation in a subsequent row in the table also employs the consistency corrections of all approximations in previous rows; thus, the K2NNC2 approximation uses all self-consistency equations noted.

Approximations
Patterns Consistency conditions where the new ghost effective fields p 1 and p 2 are obtained self-consistently from two extra equations [see Table II and Fig. 3(d)]. Now, we consider a cluster containing up to 3NN of the central site or K3NN approximation. For the case of 1NN effective mean-field corrections denoted as K3NNC1, the cluster Hamiltonian is simply  a 1 , a 2 , a 3 , b 1 , b 2 , b 3 , and b 4 . TABLE III. Different levels of approximations within the Kikuchi approach with a cluster up to 3NN of the central site. The middle column presents the patterns corresponding to the effective fields of Eq. (7). Arabic numerals correspond to single-body terms on the sites noted (and all equivalent sites). Latin numerals correspond to the patterns shown in Fig. 3. The last column shows the self-consistency equations to be solved for the effective field parameters. Note that an approximation in a subsequent row in the table also employs the consistency corrections of all approximations in previous rows; thus, the K3NNC2 approximation uses all self-consistency equations noted.

Approximations
Patterns Consistency conditions When the K3NNC2 approximation (2NN effective mean-field corrections) is taken into account, the cluster Hamiltonian is expressed as where the new ghost effective fields p 1 , p 2 , and p 2 are obtained from three new self-consistency equations [see again Table III and Figs. 3(c) and 3(d)].
In Sec. III, the accuracy of these approximations will be assessed by comparing their predictions with Monte Carlo data.

III. COMPARISON BETWEEN APPROXIMATIONS
In this section, we test our approximations against KMC simulations of the full NO oxidation system and MMC simulations of oxygen adsorption/desorption. Figure 4 depicts the average oxygen coverage as a function of the chemical potential as well as the 1NN pair-correlation as In all cases, MMC simulations of oxygen adsorption/desorption for a lattice of 96 × 96 sites at T = 480 K are plotted in red squares. We also plot in all cases results from the MF approximation as solid black lines. The predictions of the approximations become progressively closer to MMC results as the cluster size is increased. The BP approximation is already accurate, and the curves obtained with the K3NN approximations are almost indistinguishable on the plot. a function of the oxygen coverage, at a representative temperature of 480 K. We compare our approximations against MMC simulations (see the supplementary material for simulation details). As a reference, MMC simulation results (red squares, highest accuracy) and predictions from MF approximation (black lines, lowest accuracy) are plotted in all panels. Note that these equilibrium properties can also be obtained with high accuracy and low computational expense by using the matrix transfer method, 26,48 as discussed in the supplementary material.
Before starting to discuss how accurate are our methods in predicting the quantities under consideration, let us first briefly discuss the MMC simulation results. The first feature one can notice is the formation of two well-defined and pronounced plateaus due to repulsive interactions. This observed behavior can be interpreted as follows: 43,44 (i) for very low coverages, the oxygen adatoms basically do not interact; however, for coverages around 0.01 ML, the MF curve deviates appreciably from the MMC curve (see the supplementary material for a magnified view around this region). This is a strong indication of the effect of repulsive interactions, which result in the MMC-predicted correlations being close to zero, as the particles are avoiding each other in this coverage range. Above this range, the adsorption sites are filled until the formation of a ( 3)R30 • ordered phase; (ii) for an oxygen coverage between 1/3 and 2/3, the filling continues up to the formation of a similar ( 3)R30 • structure but with a surface coverage of 2/3.
Finally, at very high chemical potentials, the surface becomes almost totally covered. The flat regions associated with different structural rearrangements of the adsorbed oxygen atoms are followed by characteristic behaviors in the pair-correlation. Each plateau is associated with a change in the slope of the pair-correlation. We now proceed to discuss the capacity of our approximations to predict the adlayer structures just discussed. In general, it is clear that the MF approximation performs rather poorly. In particular, it does not predict any ordered phases or any critical transitions. This observation clearly highlights the shortcomings with this exceedingly popular approach in the catalysis community, as well as the need for more accurate computational schemes. Figures 4(a)  model after counting the NO 2 molecules produced (or consumed) per site per time. They also show the corresponding TOF predictions obtained from Eq. (14) for the different approximations of Subsections II C 2 and II C 3 and from Eq. (29) for a single-site cluster (MF approximation). Figures 5(d)-5(f) present similar plots, but for T = 680 K. In all panels, black vertical dashed lines highlight the condition of chemical equilibrium computed from the thermodynamics of the gas phase reaction for the given oxygen molar fraction (y O 2 = 0.1). 19 It is given by where y NO 2 and y NO are molar fractions. The Gibbs free energy of the overall reaction is From our simulation results, one can notice that the NO oxidation rate decreases with the NO 2 to NO ratio until the equilibrium point is reached. After that, NO 2 reduction starts to dominate. It is also interesting to see that the MF approximation overpredicts the TOF by up to 4 orders of magnitude. As mentioned above, this result emphasizes the need of more accurate methods. On the other hand, Fig. 5(a) shows that predictions from the BP with edge atom-atom interactions and their respective effective corrections are still inadequate, due to the exponential dependence of reaction rate on activation energy which amplifies errors in the activation energy predictions. Figures 5(b) and 5(c) show that the Kikuchi approximations exhibit increasingly higher accuracy as the cluster size increases and more effective mean-field corrections are incorporated in the cluster Hamiltonian. In particular, the Kikuchi approximation with up to 3NN sites of the central site and 2NN correction terms or K3NNC2 approximation perform rather well in reproducing the TOF. In the supplementary material, we also present plots of the oxygen coverage as a function of the NO 2 to NO ratio for the two considered temperatures. These plots not only represent another example of the performance of our approach, but also show that the oxygen coverage we consider is large enough to justify the assumptions of Sec. II B.
The KMC simulations presented in the Figs. 5(d)-5(f) show that, as expected, the TOF increases with temperature (in this case we use T = 680 K). As for the case of lower temperature, the MF approximation not only overpredicts the TOF, but it also fails to predict the order of the reaction (calculated from the slope of the TOF curve). The BP approximations overpredict the TOF, and the K2NNC1 approximation (see Subsection II C 3) does not seem to offer much higher accuracy than the BP. On the other hand, the K2NNC2 approximation (where 2NN corrections are taken into account) performs much better. The K3NNC1 approximation yields lower quality predictions than the K2NNC2, showing that the use of larger clusters alone is not necessarily beneficial. The K3NNC2 approximation, however, is much more accurate. For the case of NO 2 reduction, all approximations beyond K2NNC2 perform well. We have thus seen that it is possible to construct highly accurate approximations, which exhibit quantitative agreement with the KMC results. It is desirable that these approximations are also computationally efficient. We will thus compare the computational times required to obtain the results just shown, with the different approximations versus the KMC simulations.
The KMC simulations are extremely computationally intensive, primarily because of time scale separation: it may take hundreds of NO oxidation, NO 2 reduction, and oxygen diffusion events [see reactions (10) and (12) in Sec. II B] to quasi-equilibrate the adlayer to the correct chemical potential, up to the point that a single O 2 adsorption/desorption event is simulated (rate limiting process). However, the approximations have this quasi-equilibration built into the formalism and only focus on the important (rate-limiting) elementary events. Furthermore, in KMC simulations apart from detecting energetic cluster contributions every time a new possible event arises, it is necessary to deal with the updates of the lattice processes after the execution of the corresponding elementary event. To demonstrate the computational efficiency of our approximations, we present in Table IV the cost (in units of CPU seconds) to get the KMC simulation curve from the method discussed in the supplementary material and of obtaining the different theoretical curves plotted in Figs. 5(a)-5(c), each corresponding to a different approximation. This table shows that our approach offers a tremendous computational gain with negligible loss of accuracy. In particular, the most expensive of our methods, the K3NNC2 approximation offers up to 5 orders of magnitude reduction in computational time and an excellent agreement with KMC results.

IV. SUMMARY AND CONCLUSIONS
It is well-recognized that a KMC scheme incorporating detailed models for the adlayer energetics and coverage dependence of the activation barrier of reaction events is a powerful tool to investigate surface reactions. However, introducing all this complexity into KMC algorithms comes at a significant computational cost. On the other hand, efficient but rather approximate microkinetic mean-field models are another widely used approach. These two approaches have their own advantages and disadvantages. Therefore, it is imperative to develop new methodologies that can capture coverage effects in catalytic kinetics in a computationally efficient way without sacrificing accuracy.
Motivated by this need, we developed a hierarchy of approximations based on the so-called cluster mean-field approximation and demonstrated the methodology to a model of NO oxidation/NO 2 reduction on Pt(111), a reaction of practical and scientific interest. We assumed that O 2 adsorption/desorption is the rate-limiting process and used a latticegas Hamiltonian with a 3-figure cluster expansion to simulate oxygen adlayer on the surface. However, the same approach can be implemented using cluster expansions with longerrange terms. Activation energies were mapped to reaction energies (calculated from the Hamiltonian formalism) by Brønsted-Evans-Polanyi relations. Reaction rates or TOFs, oxygen coverages, and pair-correlations were derived from the cluster mean-field approximation and compared with those obtained from KMC simulations of the full reaction model and MMC simulations of oxygen adsorption/desorption.
In our cluster mean-field approximation, all relevant quantities were calculated by averaging in small clusters of the lattice, handling the neighborhood as a "cloud" of oxygen atoms that interact with the explicitly treated atoms within the aforementioned cluster. These global interactions were incorporated into the Hamiltonian formulation as effective meanfield correction terms that must be obtained self-consistently. The smallest such cluster is a single site, surrounded by the implicitly treated "cloud" of oxygen atoms. This is the wellknown Bragg-Williams or mean-field approximation in which the adsorbates are randomly distributed and uncorrelated. In the second level of approximations, a central site and its six first nearest-neighbors were treated explicitly (BP approximation). This approximation is inspired by Bethe's original work on the statistical theory of superlattices. 26 Finally, inspired also by Kikuchi's work, 26 we considered higher-level approximations in which we treated explicitly sites up to second and third nearest-neighbors of the central site (Kikuchi approximations).
We found that the predictions obtained with the MF approximation were rather poor, as spatial correlations are neglected. The approximation was unable to capture even qualitative features observed in KMC or MMC simulations. As expected, it only agreed with those computational methods at very large or small coverages. The BP approximations, on the other hand, exhibited increasingly higher accuracy as lateral interactions between adsorbates bound to the edge sites and their respective effective mean-field corrections were taken into account. Although the BP approximation performed well for predicting coverages, it failed to reproduce the quantitative behavior of the TOF since small errors in the activation energies computed with this approximation were amplified, due to the exponential dependence of reaction rate on activation energy. On the other hand, we found that the highest-level Kikuchi approximation considered yields results in quantitative agreement with the KMC. To appreciate the computational efficiency of our methodology, we compared the computational cost for different levels of approximations with KMC simulations. The enormous computational savings offered by the BP and Kikuchi methods and the higher accuracy exhibited by the Kikuchi approximation show the promise of this methodology. A detailed analysis of the approximation errors is a topic of our ongoing work.
We would like to remark that our methodology can be extended in a straightforward way to include long-range or many-body lateral interactions. For example, when considering the K3NNC2 approximation, one has to include the corresponding interaction parameters J 1 and J 2 (for 1NN and 2NN, respectively) in the approximate Hamiltonians, but the number of unknowns and the number of evaluations to compute a partition function remain the same. Extending the formalism to multiple adsorbates and site types is certainly feasible as well, though one would have to tackle a combinatorial increase in complexity when calculating partition functions and averages. On the bright side, these calculations are "embarrassingly parallel," making it possible to effectively address this challenge. Importance sampling techniques and algorithmic enhancements [parallelisation, e.g., using Field Programmable Gate Arrays (FPGAs)] would be helpful, and we will be exploring these in future research.
We would like to emphasise that these approximations are valid in the so-called fast diffusion limit. For systems whereby diffusion is not sufficiently fast to equilibrate the adlayer, one would like to use KMC anyway as the approximations do not hold. Our aim here was to establish novel approximations for the kinetic modeling of equilibrated adlayers and compare their behavior against the mean-field approach and KMC. While our study focuses on a simple system, these approximations could make the coupling of highly accurate chemistry models and fluid dynamics computationally efficient, paving the road for first-principles-based reactor design.

SUPPLEMENTARY MATERIAL
See supplementary material for a complete description of the pre-exponential factors accompanying the rate constants of our reduced NO oxidation/NO 2 reduction model, average oxygen coverage as a function of the NO 2 to NO ratio and oxygen chemical potential, the transfer matrix method applied to calculate coverage and correlations in our system, and the MMC and KMC algorithms used in our simulations.