CO oxidation on Pd(100) at technologically relevant pressure conditions: A first-principles kinetic Monte Carlo study

The possible importance of oxide formation for the catalytic activity of transition metals in heterogenous oxidation catalysis has evoked a lively discussion over the recent years. On the more noble transition metals (like Pd, Pt or Ag) the low stability of the common bulk oxides suggests primarily sub-nanometer thin oxide films, so-called surface oxides, as potential candidates that may be stabilized under gas phase conditions representative of technological oxidation catalysis. We address this issue for the Pd(100) model catalyst surface with first-principles kinetic Monte Carlo (kMC) simulations that assess the stability of the well-characterized (sqrt{5} x sqrt{5})R27 surface oxide during steady-state CO oxidation. Our results show that at ambient pressure conditions the surface oxide is stabilized at the surface up to CO:O2 partial pressure ratios just around the catalytically most relevant stoichiometric feeds (p(CO):p(O2) = 2:1). The precise value depends sensitively on temperature, so that both local pressure and temperature fluctuations may induce a continuous formation and decomposition of oxidic phases during steady-state operation under ambient stoichiometric conditions.


I. INTRODUCTION
The large difference in catalytic activity of "ruthenium" catalysts towards CO oxidation under ultra-high vacuum and ambient pressure conditions 2 has recently been controversially discussed in terms of a possible oxide formation at the surface in technologically relevant environments 3,4 . The well known much lower thermodynamic stability of the bulk oxides of Pd or Pt 5 suggests at first glance that this discussion has no direct relevance for these more noble metals that are equally, if not more often used in oxidation catalysis. Using a constrained firstprinciples thermodynamics approach 6 we could show in a preceding study 7 , however, that the situation is not that clear-cut when considering for the possibility of sub-nanometer thin oxide films, so-called surface oxides. Mapping out a wide range of (T, p O2 , p CO )-conditions we specifically found for the Pd(100) surface that the stability range of the well-characterized ( √ 5 × √ 5)R27 • ( √ 5 in brevity) surface oxide structure 8,9,10,11,12 does indeed extend up to pressure and temperature conditions relevant for technological CO oxidation catalysis (T = 300 − 600 K, p O2 = p CO ∼ 1 atm).
This insight motivates to center a detailed kinetic modeling of the surface structure and composition during steady-state CO oxidation at technologically relevant pressure conditions on the √ 5 surface oxide. More precisely, the question of a possible oxide formation at the surface can be narrowed down to addressing the stability, formation and decomposition of the √ 5 surface oxide structure in the reactive environment. Theoretically, one is then still faced with the challenge to evaluate the kinetics of the surface elementary processes over time scales of the order of seconds or longer, while accurately accounting for the evolving atomistic structure of the surface.
In this study, we meet the time-scale challenge with a first-principles statistical mechanics approach, combining density-functional theory (DFT) for an accurate description of the kinetic parameters with kinetic Monte Carlo (kMC) simulations for the evaluation of the nonequilibrium statistical mechanics problem. At present, such an approach is preferentially carried out on the basis of a lattice model for the surface in order to keep the number of required DFT-based parameters tractable. Due to the different periodicity and Pd-density of the √ 5 overlayer 7,13,14 and the Pd(100) substrate, a combined lattice model embracing both structures will be highly involved and will require numerous kinetic parameters for a manifold of conceivable atomistic pathways for the complex structural transition between the oxidized and pristine Pd surface. Before embarking on such an explicit modeling of the formation and decomposition of the surface oxide structure in the simulations it is thus appropriate to first assess whether the stability range of the √ 5 surface oxide does indeed extend to the catalytically relevant environments, when the kinetics of the ongoing CO 2 formation is fully accounted for. In this study we determine this stability range by focusing on the onset of the decomposition of the surface oxide structure when the partial pressure ratio p CO : p O2 exceeds a critical value. We identify this onset once the average O concentration in the √ 5 structure is reduced below the value corresponding to the intact surface oxide, which allows us to restrict our kMC modeling to the √ 5 lattice.
The central result of our kMC simulations, which we had briefly highlighted before 15 , is that at ambient pressures the surface oxide is stabilized at the surface up to CO:O 2 partial pressure ratios just around the catalytically most relevant stoichiometric feeds (p CO : p O2 = 2 : 1). This suggests an interpretation of the morphologi-cal changes observed in recent reactor scanning-tunneling microscopy (STM) measurements 16 under O-rich feeds in terms of a formation of the √ 5 structure at the surface. The precise value for the critical partial pressure ratio we compute depends furthermore sensitively on temperature, so that both local pressure and temperature fluctuations may induce a continuous formation and decomposition of the oxidic phase during steady-state operation under ambient stoichiometric conditions. We therefore expect that a full understanding of the catalytic activity of the Pd(100) surface in such technologically relevant environments can only be obtained from a modeling that explicitly accounts for the corresponding on-going structural changes, which will then also be able to address the actual role played by the surface oxide for the catalytic activity itself. Nevertheless the rather high peak turnover frequencies obtained in our simulations even at the intact surface oxide show already that the formation of this structure in the reactive environment can not simplistically be equated with a "deactivation" of the model catalyst.

A. First-principles kinetic Monte Carlo simulations
Our target is to describe the steady-state catalytic CO oxidation over the √ 5 surface oxide by explicitly accounting for the kinetics of the different elementary processes taking place in the system. Properly evaluating the time evolution of such a system requires simulation cells that are large enough to capture the effects of correlation and spatial distribution of the chemicals at the surface. Furthermore, most of the considered processes are highly activated and occur on time scales that are orders of magnitude longer than e.g. a typical vibration (∼ 10 −12 s). Due to these so-called rare events we need to evaluate the statistical interplay between the different elementary processes over an extended time scale that can reach up to seconds or more. Since this is unfeasible for molecular dynamics simulations, we employ kinetic Monte Carlo simulations, where the time evolution is efficiently coarse-grained to the rare-event dynamics. In a kMC simulation, each kMC-step corresponds to the execution of a rare event, while appropriately accounting for the short-time dynamics of the system in between. In the generated sequence of system states, each state is assumed to be independent of all previous states, i.e. the kMC algorithm provides an efficient numerical solution to the Master equation of a continuous time Markov process 17,18,19,20,21 .
In order to keep the number of considered elementary processes in a tractable range the kMC simulations are performed on a lattice. Obviously, the chosen lattice model has to properly reflect the important structural features of the real system. For this model a list of all relevant elementary processes p on the lattice is set up with the corresponding rate constants k p that characterize the probability to escape from the present system state by the process p. Starting in a given initial system configuration the sum over the rate constants of all possible processes in this current configuration is evaluated, k tot = p k p . One process i is then randomly chosen by with ρ 1 ∈]0, 1] being a random number. Since the probability of selecting a process i is weighted by its rate constant k i a process with a large rate constant is more likely to be chosen than a process with a small rate constant. The selected process is executed and the system configuration is updated. This way the kMC algorithm simulates a sequence of Poisson processes and an explicit relationship between kMC time and real time can be established 21 . This evolution of real time after each kMC step is given by where ρ 2 ∈]0, 1] is a second random number 21 .
The crucial ingredients to a meaningful kMC simulation are therefore the identification of all relevant elementary processes p and the determination of the corresponding rate constants k p . For the latter we rely on the modern first-principles kinetic Monte Carlo approach 22,23,24,25,26,27,28,29 and calculate the rate constants using DFT and transition-state theory (TST). Concerning the elementary processes during the CO oxidation over the √ 5 surface oxide we consider adsorption, desorption, diffusion and reaction of the chemicals, and follow the recipe of how to determine the rate constants for these processes from DFT and TST given in Ref. 29. In the next Section the resulting expressions for the rate constants are briefly reviewed. The first-principles data required to evaluate the rate constants with these expressions is then given in Section II E.
B. The rate constants

Adsorption
The rate constant of adsorption of a gas phase species i on a surface site of type st depends on the impingement rate and the local sticking coefficient, The only system-specific parameters entering the impingement rate are the mass m i of the impinging gas phase molecules and the surface area A. An explicit, accurate evaluation of the local sticking coefficientS i,st is highly involved and requires a determination of the full, high-dimensional potential energy surface (PES) of the adsorption process, on which molecular dynamics (MD) simulations have to be performed for a statistically relevant number of trajectories 30 . In the present study the local sticking coefficient is instead roughly estimated as 29 where ∆E ads i,st is the highest barrier along the minimum energy pathway (MEP) and A i,st is the so-called active area. Within the concept of an active area it is assumed that only particles of a species i with an initial lateral position within a certain area A i,st around the adsorption site st can actually stick to this site. The rate constant of adsorption in Eq. (3) will then actually become independent of the total impingement area A and only the active area A i,st has to be known. f ads i,st is a factor that reduces the number of impinging gas phase particles by the fraction that is not traveling along the MEP and might thus be reflected at some higher barrier along a different pathway.

Desorption
Since the desorption process is the time-reversed process of adsorption the desorption rate constants are connected to the adsorption rate constants by the detailed balance criterion and can thus be obtained by applying 29 where ∆G i,st (T, p i ) is the change in the Gibbs free energy between a particle in the gas phase and in the adsorbed state. ∆G i,st (T, p i ) is approximated by the difference between the chemical potential of the particle in the gas phase, µ i,gas (T, p i ) = E tot i,gas + ∆µ i,gas (T, p i ) and the free energy of the particle in the adsorbed state To evaluate the desorption rate constant using Eq. (5) we thus need to determine, in addition to the adsorption rate constant, the vibrational partition function of the particle in the adsorbed state, z vib i,st , the temperature and pressure dependent part of the gas phase chemical potential, ∆µ i,gas (T, p i ), as well as the binding energy, E bind i,st = E tot i,gas − E tot i,st of the particle at its adsorption site st.

Diffusion
In a diffusion process a particle i moves from a site st to a site st ′ . If an appropriate saddle point along the diffusion pathway can be identified, harmonic TST can be applied to obtain the corresponding rate constant of diffusion 29 and The diffusion barrier ∆E diff i,st→st ′ denotes the maximum barrier along the MEP of the diffusion process and is given by the energy difference of the transition state (E tot i,st→st ′ ,TS ) and the initial state (E tot i,st ) at site st. The reverse process of the diffusion st → st ′ is simply the backward diffusion st ′ → st and therefore the rate constants of these two processes have to fulfill the detailed balance criterion. To calculate the diffusion rate constants the factor f diff,TST i,st→st ′ (given by the ratio of the vibrational partition functions in the transition state (TS), z vib i,st→st ′ ,TS , and the initial state z vib i,st ) as well as the diffusion barriers have to be known.

Reaction
If a suitable reaction coordinate can be identified that contains a saddle point, again TST can be used to determine the reaction rate constant. In a general form the reaction rate constant can then be written as 29 with i denoting the initial state and f the final state of the reaction process. As for the diffusion events the factor f reac,TST i→f can be obtained from the ratio of the partition functions in the transition and initial state. In addition the reaction barrier, ∆E reac i→f , defined as the energy difference in the transition and initial state, has to be determined.

C. Computational Setup
The most crucial input parameters to determine the rate constants as discussed in the previous Section are the energy barriers for adsorption, diffusion and reaction, as well as the binding energies at a certain adsorption site. We calculate the total energies that are needed to evaluate these quantities using DFT and within the full-potential (linearized) augmented plane wave + local orbital (L)APW+lo method 31,32 as implemented in the WIEN2k code 33 . The exchange-correlation (xc) energy is treated within the generalized gradient approximation (GGA) using the PBE 34 xc-functional. All surfaces are simulated using the supercell approach with inversion symmetric slabs consisting of five Pd(100) layers with the reconstructed surface oxide layer plus additional O/CO on both sides. We verified that further increasing the number of Pd(100) layers to seven changes the binding energies of both adsorbates by an insignificant amount (≤ 2 meV per adsorbate). The vacuum between consecutive slabs is at least 14Å. The outermost adsorbate layers, the surface oxide and topmost palladium substrate layer are fully relaxed. All calculations of gasphase molecules (O, CO, CO 2 ) are done in rectangular supercells with side lengths of (13 × 14 × 15) bohr (for O), (13×14×18) bohr (for O 2 and CO) and (13×14×20) bohr (for CO 2 ).
The muffin-tin radii are set to R Pd MT = 2.0 bohr for palladium, R O MT = 1.0 bohr for oxygen and R C MT = 1.0 bohr for carbon. Inside the muffin-tins the wave functions are expanded up to l wf max = 12 and the potential up to l pot max = 6. For the √ 5 surface unit cell a [4 × 4 × 1] Monkhorst-Pack (MP) grid (8 irreducible k-points) has been used to integrate the Brillouin zone (BZ), and for larger surface unit cells the MP-grid has been reduced accordingly to assure an equivalent sampling of the BZ. For the gas-phase molecule calculations in the large rectangular supercells Γ-point sampling was employed. The energy cutoff for the expansion of the wave function in the interstitial is E wf max = 20 Ry and for the potential E pot max = 196 Ry. The chosen computational setup is thus exactly the same as in the preceding thermodynamic study on this system 7 . As detailed there 7 using these basis set parameters the reported binding energies per oxygen atom or per CO molecule are converged to within 50 meV, which is fully sufficient for the arguments and conclusions presented here.

D. The kMC model
Having established the theoretical framework for the kMC simulations the next crucial step is to build a suitable lattice model that represents the investigated system. For this, we use the insight gained in our preceding constrained thermodynamics study 7 that identified the √ 5 as most relevant oxidic structure under catalytically interesting gas phase conditions (p i ∼ 1 bar, T ∼ 300 − 600 K). As initially discussed in Section I we thus concentrate the kMC simulations on this particular structure and develop a first-principles kMC lattice model that allows for a proper mapping of the √ 5 surface oxide structure. √ 5 surface unit cell, as well as the (1 × 1) surface unit cell of Pd(100) are indicated. As illustrated the √ 5 surface unit-cell consists of two halves that differ only with respect to the underlying Pd(100) substrate. We find the difference in binding energies at the equivalent adsorption sites in the two halves to be very small, and correspondingly use the reduced unit-cell in the kMC simulations. In the lower panel these prominent adsorption sites and the final kMC lattice are illustrated. The three high-symmetry adsorption sites (bridge, top2f, top4f) are indicated by the yellow (light) crosses. The hollow adsorption site of the upper O atoms (red (dark) cross) is likewise included as a site in the kMC lattice, whereas the reconstructed Pd layer as well as the lower O atoms are fixed. In the final kMC lattice only the bridge and hollow sites are considered as shown on the lower left (see text).

The kMC lattice
In the upper panel of Fig. 1 a schematic illustration of the surface oxide on Pd(100) is shown. The √ 5 surface oxide structure characterized in UHV corresponds essentially to a (O-Pd-O)-trilayer of PdO(101) on top of the Pd(100) substrate 13 . Each √ 5 surface unit-cell contains four Pd atoms in the reconstructed oxide layer. In this stoichiometric termination of the surface oxide two of the Pd atoms are fourfold coordinated and two are twofold coordinated by oxygen. There are also two kinds of oxygen atoms, two of them sit on-top of the reconstructed Pd layer (upper O atoms) and two at the interface to the Pd(100) substrate (lower O atoms) forming the previously mentioned trilayer structure.
At first sight the surface oxide structure exhibits six high-symmetry adsorption sites within the √ 5 surface unit cell, where either oxygen or CO can be adsorbed. 7 These are two bridge sites (br), two twofold by oxygen coordinated top sites (top2f) and two fourfold by oxygen coordinated top sites (top4f). The difference between each of these doubly existing adsorption sites arises only from the arrangement of the underlying Pd(100) lattice. These small structural deviations seem to have little consequences on the binding at these different prominent adsorption sites. The computed binding energies of O and CO in the two bridge sites differ by less than 0.1 eV. Similarly for the two top2f and top4f sites the difference is even smaller, ∆E bind < 50 meV (a more detailed discussion of the binding energetics can be found in Ref. 7). In our kMC model, we therefore neglect these small differences and build the lattice from the half sized √ 5 surface unit-cell as shown in Fig. 1, each containing three different adsorption sites, bridge, twofold top and fourfold top.
Inspecting the computed energetics at these three different adsorption sites in more detail the kMC lattice model can even be further simplified. We find that in the top2f and top4f sites the CO is bound much weaker, i.e. by 0.3 − 0.8 eV, compared to the bridge site. In addition, the two top sites correspond only to very shallow PES minima, with diffusion barriers of less than 50 meV to a neighboring bridge site. Considering further that due to the computed strongly repulsive interactions 7 it is not possible to adsorb CO at a top site next to an already occupied bridge site, we realize that the metastable top sites will effectively never be occupied. For oxygen the situation is even more pronounced, since here the additional adsorption of an O atom at the √ 5 surface oxide is only possible in the bridge site, whereas the two top sites are not even metastable. The two different top sites are therefore disregarded in setting up the kMC lattice representation of the surface oxide.
On the other hand, considering only the adsorption into prominent high-symmetry sites at the perfect √ 5 structure is not sufficient to address the CO-induced decomposition of this lattice. For this, we assume that the surface oxide layer starts to destabilize by the removal of the upper oxygen atoms from the trilayer structure. We therefore also include the corresponding hollow sites that are completely occupied by these surface oxygen atoms in the perfectly stoichiometric √ 5 structure as explicit sites in the kMC model, cf. Fig. 1. The basis for the final kMC lattice model is thus formed by the fixed reconstructed Pd layer and the lower O atom in the PdO trilayer. On this lattice there are two different active site types (br and hol, cf. left part of the lower panel in Fig. 1). Each site can either be vacant, or filled with one oxygen atom, or one CO molecule, whereas multiple adsorption at a site is not possible.
Considering all non-concerted processes possible on this lattice we obtain a total number of 26 possible processes. This includes five adsorption processes, the unimolecular adsorption of CO into br or hol sites and the dissociative adsorption of O 2 into adjacent br-br, br-hol or hol-hol sites, as well as the corresponding five desorption processes. Furthermore there are eight site and element specific diffusion processes describing the diffusion of O/CO from br→br, hol→hol, br→hol and hol→br sites, as well as four reaction processes including the reaction of O br + CO br , O hol + CO hol , O br + CO hol and O hol + CO br . Since the reaction processes are modeled as associative desorption there are an additional four processes for the dissociative adsorption of CO 2 . This list of non-concerted processes is exhaustive with respect to the chosen kMC lattice setup. However, it can not be excluded that there might be also other, more complex, concerted processes which could have an influence on the here discussed system.

Simulation setup
All kMC simulations are performed on the final lattice model containing (10 × 10) half √ 5 surface unit cells, i.e. 100 br and 100 hol sites, with periodic boundary conditions. To ensure that the size of the chosen lattice is sufficient simulations were also performed on larger lattices containing up to 1250 br and hollow sites, respectively, without observing significant differences concerning the investigated quantities. These quantities of interest for the present study are the average surface occupations θ i,st of a particle i on a site st under steady-state conditions, which is evaluated as where ∆t n is the time step of the nth Monte Carlo move. In all simulations we checked carefully using different starting configurations and varying simulation lengths that the system actually reached steady-state, before the average surface occupations were evaluated.

E. Evaluation of the rate constants
For all 26 processes that are possible in the kMC model, the rate constants have to be determined. The respective expressions for the four different process types, i.e. adsorption, desorption, diffusion and reaction, have been described in Section II B. Here we determine the corresponding parameters entering the equations for the different rate constants, following the procedure described in more detail in Ref. 29.

Adsorption
To calculate the rate constants of adsorption in Eq. (3) we specify the active area A i,st for the impingement of each species on each site as simply half the surface area of the surface unit cell in the kMC lattice. As shown in Fig. 1 this area corresponds to half the surface unit cell of the To obtain an estimate for the local sticking coeffi-cientS i,st in Eq. (4) we perform DFT calculations to determine first of all if there are any significant barriers, ∆E ads i,st , along the MEP for any of the five different adsorption processes. For the CO molecule this was done by vertically lifting the CO in an upright geometry from a bridge resp. hollow site. For different heights of the CO molecule above its adsorption site the energy of the whole system was then minimized with respect to all other degrees of freedom . For neither of the two sites a barrier was found along this desorption resp. adsorption pathway, i.e. ∆E ads CO,br = ∆E ads CO,hol = 0.0 eV. For the dissociative adsorption of O 2 molecules a PES was mapped out in which the energy is given as a function of the bond length between the two O atoms and the height above the surface of their center of mass (socalled elbow-plot). The lateral position of the center of mass within the surface unit cell, as well as the orientation of the O 2 molecule were kept fixed. Since there are only two hollow sites within the √ 5 surface unit cell a lifting of two oxygen atoms would correspond to a complete removal of all upper oxygen atoms in the surface oxide structure due to the periodic supercell setup of the DFT calculations. Therefore, for these calculations the size of the √ 5 surface unit cell was doubled in the direction of the involved sites. We find that the dissociation of an O 2 molecule in this particular orientation appears to be non-activated. A similar result is also found for the O 2 adsorption in neighboring bridge-hollow sites, i.e. ∆E ads O2,hol−hol = ∆E ads O2,br−hol = 0.0 eV. For the adsorption in two neighboring bridge sites, though, the PES reveals a high barrier of ∆E ads O2,br−br = 1.9 eV. A second parameter entering the sticking coefficient in Eq. (4) is the factor f ads i,st representing the reduction in sticking due to the fraction of impinging molecules not following the MEP. Since the kMC simulations are performed in a temperature range of T = 300 − 600 K, where the impinging gas phase molecules can still be considered as rather slow, it is assumed that the molecules are quite efficiently steered along the MEP and thus the value of f ads i,st is approximated by f ads i,st ≈ 1 for all five adsorption processes. With these parameters the sticking coefficients are unity for the four barrier-free adsorption processes. For the activated adsorption of O 2 in two neighboring bridge sites the sticking coefficient is as small as ∼ 10 −16 even for a temperature of T = 600 K. Even if the crude way of extracting the approximate acti-vation barrier from the restricted elbow plot would thus introduce an error of many orders of magnitude in the sticking coefficient, this adsorption process will therefore still be insignificant compared to the adsorption of O 2 in two neighboring hollow or neighboring bridge and hollow sites. The results of the kMC simulations reported below where indeed found to be unaffected, even when lowering the activation barrier for this adsorption process by up to 0.5 eV.

Desorption
Having determined the five different adsorption rate constants the respective desorption rate constants can be calculated using detailed balance via Eq. (5). To a first approximation the vibrational partition function in the adsorbed state entering Eq. (5) is set to unity, z vib i,st ≈ 1, for all five desorption processes. In the gas phase this approximation would certainly be reasonable over the here investigated temperature range (T = 300 − 600 K), i.e. only the first vibrational state of the O 2 and CO molecules is excited. In the adsorbed state, though, this might be changed due to the lower frequency modes resulting from the adsorbate-substrate bond, which could then lead to a somewhat larger partition function, z vib i,st > 1, and therefore slightly smaller desorption rate constants. Since this enters only linearly into the desorption rate constant we nevertheless found that considering such increased partition functions did not significantly affect the conclusions drawn from the simulations below.
The temperature and pressure dependent part of the gas phase chemical potential ∆µ i,gas (T, p i ) is calculated by approximating the O 2 and CO gas phase with ideal gas phase reservoirs and evaluating the translational, vibrational, rotational, electronic and nuclear partition functions using statistical thermodynamics 35 . This approximation by an ideal gas only introduces a negligible error in the temperature and pressure range considered here 35 .
A further parameter that is needed to evaluate Eq. (5) is the binding energy of the particles in their respective adsorption sites, E bind i,st . Since the binding energies of O and CO in bridge and hollow sites also depend on lateral interactions between the adsorbates, we expand them in a lattice gas Hamiltonian (LGH) 36,37,38,39,40 . Aiming only at a first approximation of the effect of lateral interactions we employ a crude LGH that is restricted to the nearest neighbor pair interactions between adsorbates in different and alike sites, yielding Accordingly, the binding energy of an oxygen atom is calculated. During the associative desorption of an O 2 molecule, respectively, two neighboring sites are depleted, which can be expressed likewise using the LGH. As shown in Fig. 2 where E √ 5−2O corresponds to the total energy of the surface oxide structure without the two upper oxygen atoms. The total energy of an oxygen molecule, E tot O2 , is calculated using the relationship E tot with a highly converged value of the oxygen gas phase binding energy of E bind O2 = −6.22 eV 41 . Table I lists the determined values for the four on-site energies and 10 interaction parameters. With these parameters the binding energy of every adsorbate for any random configuration on the lattice can be evaluated using Eq. (12). With this binding energy the respective desorption rate constants can then be calculated via Eq. (5).

Diffusion
Within the kMC lattice model for the √ 5 surface oxide there are eight different diffusion processes, four between like sites (O/CO hol→hol and br→br) and four between unlike sites (O/CO hol→br and br→hol). To obtain the corresponding rate constants of diffusion the vibrational partition functions of the adsorbates in the initial state and in the transition state are needed, as well as the diffusion barriers, cf. Eqs. (6) and (7).
With a similar justification as for the desorption rate constants, we assume the partition functions in initial and transition state to be comparable, i.e. z vib i,st ≈ z vib i,st→st ′ ,TS yielding f diff,TST i,st→st ′ ≈ 1 for all eight diffusion processes. To calculate the diffusion barriers, ∆E diff i,st→st ′ , we determine the transition state by defining a reaction coordinate connecting the initial and final state and calculate the energy of O/CO for different positions along this coordinate. Since we found the surface oxide trilayer to be rather mobile with respect to a registry shift parallel to the Pd(100) substrate 7 , the lateral positions of the reconstructed Pd atoms are additionally kept fixed in these calculations. The resulting barriers for these approximate diffusion paths are listed in Table II Fig. 3) or -for the diffusion between unlike sites or like sites with a different nearest neighbor coordination -the energies of the initial and final state wells differ (right part of Fig. 3). In the latter case the values listed in Table II represent the minimum value of the diffusion barrier, i.e. if the final state has the same or a lower energy than the initial state the diffusion barrier is equal to the tabulated value. If on the other hand the final state is higher in energy the difference in binding energies between the initial and final state is added to the listed barrier to ensure that the detailed balance criterion is fulfilled for the two time-reversed diffusion processes. Since this difference in binding energies depends on the nearest neighbor coordination of the two lattice sites involved in the diffusion process, lateral interaction effects are this way also considered in the diffusion process rate constants.
The described approximations in determining the exact transition states lead to uncertainties in the stated barrier values which in turn influence the corresponding diffusion rate constants. We checked on the importance of the various diffusion processes on the simulation results and found them to play only a minor role. Essentially, even completely switching off all diffusion processes led only to insignificant changes of the results presented below. We are thus confident that the rather coarse procedure to determine the diffusion barriers is sufficient with respect to the current purpose of the simulations.

Reaction
To calculate the rate constants for the four possible reaction processes (O br +CO br , O hol +CO hol , O br +CO hol and O hol +CO br ) the factor f reac,TST i→f ≈ 1 in Eq. (9) is approximated by unity. Again it is thus assumed that to a first approximation the partition functions of the initial and transition state are in the same order of magnitude, and the ratio of the two is thus close to one.
To determine the reaction barriers, ∆E reac i→f , the corresponding transition states are located by dragging the two reactants towards each other and mapping out the corresponding 2D-PES (O br +CO hol and O hol +CO br ). Hereby only the coordinates of the reactants along the dragging direction are kept fixed whereas all other coordinates are fully relaxed. Similar to the problem encountered when determining the diffusion barriers the reconstructed palladium layer showed a rather high mobility with respect to a lateral shift along the substrate following the dragging direction. To circumvent this problem a two-step approach was chosen to nevertheless identify appropriate transition states. In the first step the lateral coordinates of the reconstructed palladium layer were additionally fixed and the energy is only minimized with respect to the remaining coordinates. In a second step the identified transition state geometry is relaxed with respect to the lateral coordinates of the topmost Pd layer while fixing the coordinates of the reactants O and CO, to release some of the lateral stress that was induced by fixing the lateral positions of the Pd atoms in the first step. For the reaction of O br +CO br and O hol +CO hol this two-step approach was not necessary. Similar to the dissociative adsorption of O 2 an elbow plot was calculated giving the energy as a function of the O-CO distance and the height of the corresponding center of mass above the surface. Due to the symmetry of the calculated geometries no shift of the substrate layer could occur in this case. All correspondingly calculated reaction barriers are summarized in Table III. For the reaction processes the energies of the initial, transition and final states are also influenced by the nearest neighbor interactions. Inspecting the geometries of the transition states they appear to be more related to the initial states and we therefore assume that the transition states are similarly affected by these interactions as the initial states. Thus, the energy difference between the initial and transition states is not influenced by the nearest neighbor interactions and the reaction barriers are constant with respect to different configurations on the surface.
Since the reaction processes are modeled as associative desorption the corresponding time reversed process is the dissociative adsorption of CO 2 into the respective sites. The adsorption barrier for such a process can be obtained from the reaction barrier and the corresponding binding energies of O and CO on the surface with respect to the free CO 2 molecule. The lowest of the four resulting barriers for the dissociative adsorption of CO 2 (leading to CO in bridge and O in a neighboring hollow site) is still very large with ∆E ads CO2,br−hol = 1.4 eV. At a temperature of T = 600 K the local sticking coefficientS CO2,br−hol is even in this case then only of the order of 10 −12 . Even if the CO 2 generated at the surface is not readily transported away and a noticeable CO 2 pressure would build up close to the surface, the re-adsorption of CO 2 on the surface is therefore still negligible due to the very low sticking coefficient and is thus not explicitly considered in the kMC simulations.

III. RESULTS
Having established the setup of the kMC lattice model and the respective rate constants for the 26 different processes based on DFT energetics and harmonic TST we proceed to the results of the kMC simulations. Although the established kMC model has a wider range of applicability, we focus here on the stability of the thin surface oxide layer during steady-state catalytic CO oxidation at the surface.

A. Reproducing the thermodynamic "phase" diagram
In order to make contact to our preceding constrained thermodynamics work 7 , we perform the kMC simulations in a first step only considering adsorption, desorption and diffusion processes, whereas the reaction processes are excluded. The occupation of the different sites at the surface under steady-state conditions reflects then the same situation as a surface in a constrained thermodynamic equilibrium with an oxygen and CO gas-phase, but including the effect of configurational entropy 6 . The results of the constrained thermodynamics approach are summarized in Fig. 4 (for a detailed discussion see Ref. 7). In the shown "phase" diagram the most stable surface structures for any given chemical potential of the O 2 and CO gas-phase are presented (the respective pressure scales for T = 300 K and 600 K are given in the upper x-axes and The differently shaded areas mark the stability regions of various surface structures for a given chemical potential of the O2 and CO gas-phase. Blue/solid areas include O/CO adsorption structures on Pd(100), red/white-hatched areas comprise surface oxide structures and the grey/crosshatched area represents the stability range of PdO bulk oxide. For convenience the dependence of the chemical potential of the two gas-phases is translated into pressure scales for T = 300 K and T = 600 K (upper x axes and right y axes). For the stability region of the surface oxide structures additionally schematic figures of the different surface configurations are shown. Grey, large spheres represent Pd atoms, red (dark) small spheres oxygen atoms and yellow (light) small spheres CO molecules. The white arrows indicate the pressure conditions of the kMC simulations for a temperature of T = 300 K, T = 400 K and T = 600 K. The oxygen pressure is always fixed at pO 2 = 1 atm, while the CO pressure is varied between 10 −5 ≤ pCO ≤ 10 5 atm.
right y-axes). The stability regions of the different surface structures are denoted by differently shaded areas. For the here discussed kMC simulations it is particularly important that the different "phases" can be subdivided into three groups. All areas colored in blue/solid represent structures with O and CO adsorbed on the pristine Pd(100) surface, all red/white-hatched areas comprise structures involving the reconstructed √ 5 surface oxide, and the grey/crosshatched area marks the stability region of the bulk oxide. As discussed in Section II D the here chosen kMC lattice model is based on the lattice structure of the √ 5 surface oxide and does therefore not allow for changes between this structure and the Pd(100) surface or the bulk oxide. Consequently, the constrained kMC simulations can only reproduce the part of the "phase" diagram where structures based on the surface oxide are most stable, i.e. the red/white-hatched region in Fig. 4.
Within the constrained kMC simulations we now evaluate the steady-state average occupation of bridge and hollow sites and compare the obtained structures to the structures predicted within the thermodynamic surface "phase" diagram. As mentioned above the two ap-proaches should yield equivalent results, except at the "phase" boundaries where configurational entropy is expected to become important 6 . The simulations are performed for (T, p)-conditions that cover the entire stability region of the surface oxide (red/white-hatched region in Fig. 4). Specifically, this means for T = 300 K pressures of p O2 = 10 −10 − 10 10 atm and p CO = 10 −10 − 1 atm, and for T = 600 K pressures of p O2 = 1 − 10 10 atm and p CO = 1 − 10 6 atm, respectively. We find that the kMC simulations reproduce accurately all four different surface oxide "phases" that appear in this part of the constrained "phase" diagram. In all four "phases" all hollow sites are always filled with oxygen. The bridge sites are either empty, half filled with O or CO, or completely filled with CO depending on the respective gas-phase conditions and as detailed by the schematic representations in Fig. 4. As expected from the common underlying DFT energetics, we therefore find the results of the constrained kMC simulations to be fully consistent with those obtained previously within the constrained thermodynamics approach. The more extended sampling of configurations due to the LGH expansion within the kMC simulations does not lead to the appearance of new, stable ordered structures in the constrained "phase" diagram and thus confirms that the relevant configurations were included in the preceding thermodynamic approach. The effect of configurational entropy in the constrained kMC simulations is thus verified to "only" smear out the transitions between the different ordered "phases" at the finite temperatures considered 6 , but does not affect the location of the "phase" boundaries themselves. Due to the hereby confirmed consistency between the two approaches, we can directly compare the results in the next Section, when allowing the reaction events to occur in the kMC simulation. This brings us in the position to directly address the effect of the kinetics of the ongoing chemical reactions on the surface structure and composition.

B. Onset of surface oxide decomposition under reaction conditions
By allowing the four different reaction processes to occur in the kMC simulations, we now proceed to evaluate the stability of the surface oxide structure under reaction conditions. Since the kMC lattice is fixed to the structure of the √ 5 surface oxide, the explicit structural change corresponding to the transition from the surface oxide to a pristine Pd(100) surface as predicted in the thermodynamic "phase" diagram (red/white-hatched areas → blue solid areas in Fig. 4) cannot directly be addressed. Instead, we start our simulations always in a (T, p)-range where the surface oxide is definitely the most stable "phase" (high oxygen and low CO chemical potential), and where our kMC lattice model is certainly valid. Keeping the oxygen pressure and the temperature fixed we then increase the CO pressure as indicated by the white, vertical arrows in Fig. 4. Exceeding a certain pressure ratio of p O2 /p CO , we expect the onset of the decomposition of the √ 5 structure to be characterized by a depletion of the otherwise always essentially fully with oxygen occupied hollow sites. In the kMC simulations we therefore monitor the average occupation of hollow sites by oxygen atoms, θ O hol , during steady-state. If all hollow sites are occupied by oxygen, i.e. θ O hol ≈ 100%, then the overall surface composition clearly represents the intact surface oxide structure (this is the case in all four "phases" of the surface oxide structure in Fig. 4). If on the other hand at increasing p CO oxygen atoms get increasingly depleted beyond the amount of thermal fluctuations that create isolated, uncorrelated vacancies or double vacancies due to oxidation reaction or oxygen desorption events, then the deficiency of O atoms at the surface might destabilize the √ 5 surface structure and lead to a local lifting of the surface oxide reconstruction. Here, we therefore interpret a drop of the average occupation below 90% as the onset of a destabilization of the surface oxide. Since the average occupation is a global measure over the whole simulation cell this quantity would be a bad indicator, if there was an appreciable effective attraction between formed O vacancies. A decrease in the average occupation of hollow sites by oxygen of 10% could then involve a complete, local depletion of 10% of the surface oxide area. In our kMC simulations such a locally concentrated depletion of oxygen atoms is not observed. Instead, we find the formed oxygen vacancies to be distributed over the entire simulation cell and are thus confident that monitoring the average occupation, θ O hol , is a suitable measure to determine the onset of surface oxide decomposition.
The simulations are performed for three different temperatures, T = 300, 400 and 600 K. For each temperature the oxygen pressure is fixed to p O2 = 1 atm and the simulations are run for different CO pressures between p CO = 10 −5 − 10 5 atm covering the gas-phase conditions marked by the white arrows in Fig. 4. In Fig. 5 the corresponding simulation results are summarized. The steady-state average occupation of hollow sites with oxygen, θ O hol , is plotted vs. the CO pressure. The solid (green) lines represent the results obtained by using the calculated reaction barriers as listed in Table III. As has been discussed in Section II, several, partly crude approximations were employed in the determination of the individual rate constants entering the kMC simulations. Before we proceed to analyze the obtained simulation results, we therefore critically assess the propagation of these underlying approximations and uncertainties to the final kMC results. We checked on this error propagation by systematically varying the rate constants over several orders of magnitude (corresponding to variations in the respective energetic barriers of ±0.2 eV) and each time monitoring the effect on the average occupation displayed in Fig. 5. In almost all cases, variations of the rate constants had an insignificant effect, which as already described concerns in particular all diffusion processes. Only if altering the barrier for the reaction process of an oxygen in a hollow site with a CO in a neighboring bridge site a noticeable change in the average surface populations occurs. To get an estimate for the corresponding uncertainties in our presented results, we therefore include in Fig. 5 additionally the results for simulations using a 0.1 eV lower barrier for this reaction process O hol +CO br , i.e. ∆E reac O hol +CO br = 0.8 eV, while all other barriers have been left unchanged (dashed (orange) lines). The differences between these results (dashed (or-ange) line) and the results for the unmodified barrier (solid (green) line) provide then an indication of the uncertainty underlying our approach.
Discussing first the top graph in Fig. 5 it can be seen that for a temperature of T = 300 K and an oxygen pressure of p O2 = 1 atm the surface oxide is clearly stabilized for CO pressures up to p CO < 10 −1 atm, i.e. under such gas-phase conditions θ O hol > 95 %. If the CO pressure is further increased the O population at hollow sites starts to decrease and for CO pressures of p CO > 1 atm the surface oxide is certainly completely destabilized (θ O hol ≈ 0 %). For the slightly lower reaction barrier (dashed (orange) line) a steeper decrease in the oxygen population is found, and the surface oxide is stabilized for CO pressures of p CO < 10 −2 atm. To compare these results to the ones obtained within the constrained thermodynamics approach the stability region of the surface oxide as determined in the "phase" diagram in Fig. 4 is indicated by the vertical black lines in the three graphs of Fig. 5. Following the white arrow in the "phase" diagram in Fig. 4 for T = 300 K and p O2 = 1 atm this boundary between the surface oxide structure and a CO covered Pd(100) surface is reached for a CO pressure of p CO ≈ 2 atm. Within the constrained thermodynamics approach the surface oxide will thus be stable for p CO 2 atm, whereas for p CO 2 atm the CO covered Pd(100) surface is the most stable "phase". Comparing this to the kMC simulations which include the effect of the CO 2 formation kinetics, we see that the depletion of oxygen atoms in hollow sites and thus the destabilization of the surface oxide structure starts already at slightly lower CO pressures compared to the transition from the stability region of the surface oxide to the one of a CO covered Pd(100) surface in the "phase" diagram. Although we find the kMC simulations to be dominated by adsorption/desorption processes of the reactants (as anticipated in the constrained thermodynamic approach assuming equilibrium with the gas-phase reservoirs), the kinetics of the ongoing oxidation reaction thus still leads to a slight decrease in the stability of the surface oxide. In turn the kMC simulations predict that at a temperature of T = 300 K oxygen-rich conditions (p CO /p O2 ≈ 1 : 10) are needed to stabilize the surface oxide structure.
In Fig. 5(b) corresponding results are shown for T = 400 K. Compared to the simulations at T = 300 K much more reaction processes take place in the kMC simulations at this higher temperature. For CO pressures of p CO > 10 atm though almost no reaction processes are detected in the steady-state, since under these conditions the surface is already completely filled with CO. Then again this corresponds to gas-phase conditions where also in the constrained thermodynamics approach a CO covered Pd(100) surface is the most stable structure (cf. Fig. 4). Similar to the results for T = 300 K the surface oxide structure is completely stabilized for CO pressures of p CO < 10 −1 atm. At a pressure ratio of p CO /p O2 = 2 : 1, i.e. here p CO = 2 atm, there are still 94 % of all hollow sites occupied by oxygen (solid (green) line in Fig. 5(b)). This result depends, however, critically on the barrier for the reaction of O hol +CO br . If this barrier is decreased by only 0.1 eV the occupation drops to 70 %, and at an even lower barrier of ∆E reac O hol +CO br = 0.7 eV θ O hol decreases even further to ∼ 30 % under these conditions. Whether or not the surface oxide prevails at this temperature for the catalytically most relevant stoichiometric partial pressure ratio can thus not be safely concluded in light of the uncertainties of our approach. For CO pressures lower than p CO < 10 −1 atm, however, the surface oxide is definitely stabilized, just as much as it is definitely destabilized at CO pressures of p CO > 10 atm. We can therefore determine the onset of the surface oxide decomposition to fall into the pressure range 0.1 atm < p CO < 10 atm. Even considering the uncertainties in the reaction rate constants and the limitation to a single kMC lattice, these results agree thus rather nicely with the reactor STM experiments by Hendriksen et al. 16 which have been performed under similar temperature and pressure conditions (T = 408 K, p tot = p CO +p O2 = 1.23 atm). Depending on the pressure ratio of oxygen and CO Hendriksen et al. observed a change in the morphology of the surface, which was assigned to a change from an oxidic to a reduced surface structure and oxygen-rich feeds were required to stabilize the oxidic structure.
The simulation results for the highest considered temperature of T = 600 K are shown in Fig. 5(c). We find that for this elevated temperature the surface oxide is now actually stable up a rather sizable CO pressures. Only for CO pressures as high as p CO = 10 atm the surface oxide starts to decompose (θ O hol < 95 %), which already coincides almost with the border of its stability as determined in the constrained thermodynamics approach (black vertical line). In the "phase" diagram in Fig. 4 it can be seen that for a temperature of T = 600 K and an oxygen pressure of p O2 = 1 atm CO is not stabilized in large amounts as adsorbate on the surface oxide. The CO is thus not readily available as adsorbed reaction partner and the stability of the surface oxide is hardly affected by the ongoing catalytic CO 2 formation. Up to the stoichiometric pressure ratio of p CO /p O2 = 2 : 1 no decomposition of the surface oxide structure occurs, and as apparent from Fig. 5(c) this result holds even in light of the estimated uncertainty of the underlying DFT energetics. Comparing the critical p CO /p O2 ratio determined for the decomposition onset in the temperature range T = 300 − 600 K, we can thus clearly identify an increasing stability of the surface oxide with increasing temperature, which at the highest temperatures studied reaches well up to the catalytically most relevant feeds. Furthermore we find that under these feeds the simulated turnover frequencies for the intact surface oxide alone are already of similar order of magnitude as those reported by Szanyi and Goodman 42 for the Pd(100) surface under comparable gas-phase conditions. While a quantitative comparison is outside the scope of the present work, we note that contrary to the prevalent general preconception at least this particular surface oxide is thus clearly not "inactive" with respect to the oxidation of CO.

IV. SUMMARY
We have employed first-principles kMC simulations to investigate the stability of the √ 5 surface oxide structure at Pd(100) against CO-induced decomposition under steady-state operation conditions of the CO oxidation reaction. The focus has been particularly set on the thin surface oxide, since this structure was identified as the most stable oxidic "phase" under catalytically relevant gas-phase conditions in a preceding constrained first-principles thermodynamics study 7 . The developed kMC model describes the CO oxidation on the √ 5 surface oxide, accounting for all non-concerted processes that are possible on the chosen kMC lattice. This includes 26 site and element specific adsorption, desorption, diffusion and reaction processes. The respective kMC rate constants have been obtained from DFT total energy calculations together with harmonic transition state theory.
Monitoring the onset of the surface oxide decomposition at increasing CO pressures, the central result of the kMC simulations is that much in contrast to thicker bulk-like oxide films the stability of the √ 5 surface oxide can well extend to gas-phase conditions that are relevant for technological CO oxidation catalysis. While the accuracy of present-day DFT functionals does not permit us to conclude on precise values for the critical CO:O 2 partial pressure ratio above which decomposition sets in, our results show undoubtedly that at ambient pressures this ratio is close to the most interesting stoichiometric feed conditions (p CO : p O2 = 2 : 1) and will furthermore shift decisively towards more CO-rich conditions with increasing temperature.
On the basis of these results we conclude that the surface oxide structure can actually be present under catalytic reaction conditions at ambient pressures and our simulated turnover frequencies indicate that it is then also catalytically active towards the oxidation of CO. For the catalytically most relevant feeds the system appears furthermore to be rather close to a transition between the √ 5 surface oxide structure and a reactant covered Pd(100) surface. Particularly in stoichiometric to slightly O-rich feeds even local fluctuations in both temperature or the reactant partial pressures might then lead to significant local or global changes in the structure and composition of the surface. Under steady-state operation conditions oscillations between these two states are thus conceivable, so that in its active state the catalyst surface might exhibit parts that correspond to Pd(100) covered by the reactants, as well as by patches of surface oxide, and potentially it is precisely the on-going formation and decomposition of the oxidic structure that is key to understand the catalytic function of this surface under such technologically relevant pressures. For a full understanding of the catalytic CO oxidation at the Pd(100) surface in corresponding environments it thus appears not to be sufficient to concentrate on either the reduced Pd(100) surface or on an oxidic structure, but both surface states as well as the morphological transition between the two will have to be included in an appropriate model.

V. ACKNOWLEDGMENT
The EU is acknowledged for financial support under contract NMP3-CT-2003-505670 (NANO 2 ), and the DFG for support within the priority program SPP1091.