A More Accurate Kinetic Monte Carlo Approach to a Monodimensional Surface Reaction: The Interaction of Oxygen with the RuO2(110) Surface.

The theoretical study of catalysis would substantialy benefit from the use of atomistic simulations that can provide information beyond mean-field approaches. To date, the nanoscale understanding of surface reactions has been only qualitatively achieved by means of kinetic Monte Carlo coupled to density functional theory, KMC-DFT. Here, we examine a widely employed model for oxygen interaction with the RuO2(110) surface, a highly anisotropic system. Our analysis reveals several covert problems that render as questionable the model's predictions. We suggest an advanced approach that considers all the relevant elementary steps and configurations while smoothing the intrinsic errors in the DFT description of oxygen. Under these conditions, KMC provides quantitative agreement to temperature-programmed desorption experiments. These results illustrate how KMC-based simulations can be pushed forward so that they evolve toward being the standard methodology to study complex chemistry at the nanoscale.

Introduction

Understanding microscopic processes at the nanoscale is one of the crucial developments required to achieve atomistic engineering in catalysis. Methods based on density functional theory, DFT, that determine reaction profiles have been proven useful, but as the systems become more complex, the analysis is less straightforward. Mean-field microkinetic approaches based on DFT allow the integration of time scales and can provide coverages, reaction orders, apparent activation energies, and descriptors,[1,2] but that might not yet be enough. In this perspective, the successful implementation of multiscale methods based on kinetic Monte Carlo (KMC) techniques with parameters obtained through DFT represents the most robust pathway to reach fully consistent ab initio simulations. The main challenge ahead of KMC is to derive quantitative data directly comparable to experiments, an objective that seemed far-fetched just a few years ago.[3]

Compared with other mean-field approaches, KMC properly addresses surface anisotropy and lateral interactions[4,5] but implies a high computational burden because initial, final, and transition states for a wide variety of local configurations need to be located and updated for different local distributions.[6] In addition, the accuracy problems that are intrisic to DFT are amplified when calculating rates,[7] and in many cases, the calculation of the prefactors through statistical mechanics is either overlooked, simplified, or reduced.[8,9] As a consequence, despite the enourmous potential of KMC for the study of catalytic reactions on materials, its application has been confined mostly to qualitative interpretations.[3] In the present manuscript, we show how for a truly monodimensional system the kinetic Monte Carlo approach is accurate enough to address the most stringent tests based on simulation of temperature-programmed desorption data, provided that (i) all the relevant structures are taken into account, (ii) the intrinsic errors of DFT are alleviated, (iii) full statistical thermodynamic contributions in the transition state theory are considered, and (iv) a sufficiently complex set of lateral interactions is taken into account.

To this end, we have chosen a highly anisotropic system, with industrial interest and for which a large amount of experimental data has been collected. RuO2 is a powerful catalyst in oxidation reactions[10,11] that has industrial applications in the Deacon process.[2,12] The most common (110) exposed surface presents a large anisotropy, thus making it a wonderful scenario to test KMC techniques. From the very beginning[13] of active research on RuO2, in-silico simulations[3,7,14,15] were found essential to understand[12,16,17] its structural and catalytic features. Naturally, the interaction of RuO2(110) with oxygen is the keystone in every oxidation process on the surface and, considering the number of publications[9,18−22] that cover various aspects of this interaction, one may expect the subject to be well-understood and generally closed. Below, we show that such a conclusion is premature.

To the best of our knowledge, all existing models[9,18,23] of chemical processes on RuO2(110) treat oxygen adsorption and desorption from the surface as single step elementary reactions. At the same time, it is well known[22] that oxygen may be molecularly adsorbed to the surface, and it may stay there in different possible configurations[20] prior to its dissociation. Is it a must to be accounted for in a reliable model of RuO2(110) surface chemistry and whether it may be safely omitted in some cases has not been discussed.

Model and Methods

To approach the above-stated question, we employ ab initio simulations of RuO2(110) within density functional theory and kinetic Monte Carlo techniques. DFT was performed in a p(1 × 2) supercell (Figure 1). The data obtained were post-treated by cluster expansion to reconstruct the energies of larger pieces of the surface and to reveal the binding energies (Figure 2; the Supporting Information contains additional details) of oxygen molecules and atoms at the cus sites of the surface. Our actual DFT values, shown at the figure by the black profile, are somewhat different from those reported in ref (20) as a result of the use of different DFT functionals (RPBE instead of PW91). RPBE is known to provide more accurate adsorption energies compared to any other GGA functional.[24] Still, available DFT methods are not accurate in the prediction of oxygen dissociation energy because of degeneration of the O2 ground state in terms of both spin and orbit (3∏u and 3P, for O2 and O, respectively). The error is estimated[25] as high as ζ = 0.56 eV in the case of the RPBE functional. The value with our setup was reduced to 0.48 eV, but it is still large enough to significantly impact outcomes of KMC simulations; thus, it must be compensated. From a comparison to previous ways to overcome the problem[26] and by analyzing the sensitivity of the KMC simulations series to this particular parameter when compared with experimental data on temperature-programmed desorption (TPD) of oxygen from RuO2(110), we learned that a shift of the DFT-calculated energy of oxygen in the gas phase by 0.39 eV (70% of ζ) is necessary. Because the molecularly adsorbed oxygen still bears some leftover spin, we have correspondingly shifted the DFT-calculated energies of O2* (mono) and O2** (dihapto) configurations by the values 0.08 and 0.04 eV, proportional to the gas-phase correction and the remaining magnetic moment of the adsorbate. The resulting corrected energy profile is shown in red in Figure 2, and it is referred to in the rest of discussion.

Figure 1

Schematic illustration of a RuO2(110) surface and of its representation as a KMC simulation lattice. The dashed rectangle shows the unit cell’s size. Under normal conditions, all the bridge sites are occupied by oxygen atoms.

Figure 2

Energy levels of O2 on a RuO2(110) surface (center), calculated for various configurations of molecular (O2, O2*, O2**) and atomic (O*) oxygen, illustrated in the bottom. All the values are given in eV. The black profile shows the actual output of DFT, and the red one corresponds to an empirically corrected one (more details in the text). The energies for O2* and 2O* configurations split into several levels (α, β, γ) for different oxygen coverages of the neighboring surface sites are shown on the left of the illustration (the red frames show the targeted adsorption sites). Finally, the one-step profile at the right corresponds to the oxygen adsorption and desorption representation in ref (18).

Let us consider adsorption of an oxygen molecule to a cus surface site, surrounded by oxygen atoms, as shown at the subimages α, β, and γ of Figure 2. Once it has approached the site, the O2 molecule binds by one of its atoms to the Ru atom corresponding to the site (configuration O2*), and it gains 0.78 (case α) or 0.58 eV (cases β, γ), compared with the energy in the gas phase. In the absence of close neighbors in both directions along the cus row, the O2* state is not stable, as reorientation of the molecule parallel to the surface with connection of its second atom to the neighboring Ru (configuration O2**) gains an extra 0.16 (case α) or 0.36 eV (cases β, γ), without any noticeable energy barrier on the way. Further dissociation of the oxygen molecule is controlled by a 0.57 eV activation barrier; once it is crossed, the system finds itself at the energy level 1.43–1.58 eV below the initial O2 state. To sum up, at moderate coverages of O*, further O2 adsorption and desorption turn out to be the two-step process O2 ⇄ O2** ⇄ 2O*, rather than the one-step process, O2 ⇄ 2O*, adopted in refs (9), (18), and (23) and in numerous other works of the same research groups. One has to bear in mind that crowded surfaces (or, at least, large occupation of active sites) are expected under catalytic conditions.[2] We also note that the O2* state permits quite an easy diffusion of a separate oxygen molecule along the cus row via “flip-flops”, O2**⇄ O2*, with only a small energy barrier of 0.12 eV. For convenience, in the rest of article, we refer to the our two-step model of oxygen desorption as M-II, and to the one-step one as M-I.

To gain a deeper insight into the mechanics of oxygen adsorption, we need to calculate the rates of elementary reactions identified above. According to transition state theory,[9,27] the rate of an elementary reaction iswhere k and h are Boltzmann and Planck constants, T is the temperature of the system, q′ and q0 are the partition functions of the system in the transition and initial states of the reaction, and ΔE is the activation energy. In the case of the molecular adsorption step from the gas phase to the surface (O2→ O2** in our model, and O2 → 2O* in the one-step adsorption models), expression 1 must be additionally multiplied[27] by pV/kT, where p and V are the pressure and volume of the gas phase. Once it is done, and the partition functions q′ and q0 are properly expanded, expression 1 turns[27] into the well-known result of kinetic gas theoryfor nonactivated adsorption, where A is the area of adsorption site, and m is the mass of adsorbing molecule.

Fair calculation or experimental-based fitting of partition functions for the molecules adsorbed at the surface is far more challenging than obtaining the adsorbates’ binding energies because it demands the evaluation of all relevant vibrational frequencies of the system in the initial and transitional states of the considered elementary reaction. It is common practice, usually employed in the implementations[9,18,23] of M-I, to assume that q′ and q0′ are equal to unity for the molecules and atoms adsorbed to the surface, or in transition states corresponding to on-surface diffusion. In our case, to reach a better accuracy, we calculate explicitly the vibrational frequencies of the adsorbates in O2** and O* states, as well as in the dissociation transition state of the elementary reaction O2** ⇄ 2O*, and use these data to calculate all partition functions in our M-II model (check the Supporting Information for details).

Results and Discussion

The calculated rates are presented in Figure 3. For the best interpretation, they should be examined alongside the temperature desorption profile (Figure 4) corresponding[22] to the desorption of oxygen from cus sites. Desorption becomes possible around 400 K, when the rate of desorption from the O2** state (red line on Figure 3) overrides the rate of dissociation from that state (green line). However, the net rate of desorption is limited by the formation of oxygen molecules from separate atoms (thin black line on Figure 3), because O2** state itself is very short-lived at temperatures around 400 K. Comparison of this desorption-limiting rate to the oxygen desorption rate according to M-I, as implemented in ref (18) (thick line on Figure 3), reveals why previously reported simulations of the O2 TPD peak looked successful. The two lines (ref (18) and the present model) intersect at T ≈ 470 K, and they are very close to each other in the vicinity of the TPD peak. On the other hand, the difference rapidly grows with small temperature changes. It is not that important for the lower temperatures because there, both lines show that desorption of an oxygen molecule from surface is quite a rare event, but for temperatures above 500 K, the two models predict radically different lifetimes of the O2 molecule on the surface. Thus, for T = 600–800 K the average time spent by an oxygen molecule on the cus-row of the surface is two-3 orders of magnitude longer, according to M-II, than it is predicted by M-I. It should be noticed that the industrially implemented Deacon reaction (HCl oxidation) runs slightly below 600 K.[12] Thus, the lifetime differences between the two models is too large to rely on predictions of M-I when it is used in the simulations of complex reaction networks on RuO2.

Figure 3

Temperature dependences (on top) of the rates of elementary reactions involved in oxygen adsorption to RuO2(110): O2 → O2** (solid blue for p = 1 bar, dashed blue for p = 10–10 bar); O2** → O2 (red); O2** → 2O* (green); 2O* → O2** (oxygen association, thin black). The thick black line shows the oxygen desorption rate according to our M-I model. The bottom plot shows the absolute deviation of the black thick line from the thin solid black one.

Figure 4

O2 TPD profile of oxygen in 350–550 K range of temperatures. Experimental[22] data (dashed) are shown alongside with results obtained by KMC simulations in ref (18) (thick black) and by our own KMC simulations with use of the model described in the present work (thin black line). The red line shows a special postprocessing, explained in the text, of our simulation results, that shows a perfect match with the experimental measurements. All the peaks are normalized for the same area.

At this point, we want to highlight several issues that are not obvious from the discussion so far. First of all, an additional argumentation on the inconsistency of M-I model. The way the rate of 2O* → O2 reaction is calculated[9,18,23] within M-I, with the application of the microreversibility condition between 2O* → O2 and O2→ 2O* reactions and with the assumption of nonactivated adsorption, implies that the transition state of the molecule during adsorption and desorption is similar to the state in the gas phase and that the energy of the molecule in this state is the same as in the atmosphere. Somewhere along the way from this transition state to the surface, the molecule must split into separate atoms. With the strength of the oxygen’s double bond of De ≃ 5.20 eV,[28] this dissociation is possible only due to the catalytic effect of the surface. However, such interaction demands a very special binding of O2 to the surface (namely, it should be in the state similar to the transition state between O2** and 2O* in M-II). Reaching such a transition state will be the bottleneck of the dissociative adsorption and associative desorption because the actual “effective” transition state is very different from the one implicitly assumed in M-I.

We also have to emphasize the errors brought into the rates by the unjustified assumption that partition functions of molecules and atoms adsorbed to the surface are approximately equal to unity. In Figure 5, we compare the rates of the O2** → O2 and 2O* → O2** reactions of M-II as we calculate them (solid lines) with their values calculated with approximations of partition functions usually applied in implementations of M-I. In addition, we compare them with the values obtained by shifting the energy barrier values by ±0.1 eV (dotted lines). In the most interesting range of temperatures for practical applications, T = 200–800 K, the 0.1 eV error in energies leads to subsequent errors in the rates of 1–2 orders of magnitude. The approximation of the partition functions leads to a smaller error for the oxygen association reaction, but again, the error associated with the desorption of O2** → O2 has a similar size, 1–2 orders of magnitude, for the temperatures below 400 K and above 700 K. The errors demonstrated in these two examples can be tolerated in many particular cases, but they obviously also may lead to crucial mistakes under other circumstances. Thus, in general, an exact consideration of partition functions is essential. Under reaction conditions, some kind of compensation[29] between the rates and the coverages can be expected; however, for very complex reactions, it would be impossible to assess the correct behavior of this error cancellation.

Figure 5

Rates of O2** → O2 (red) and 2O* → O2** (black) reactions, calculated according to our implementation of the M-II model (solid) are compared with the approximate values obtained with partition functions of adsorbed molecules and atoms set to unity (dashed) or with energy barrier value modified by ±0.1 eV (dotted). The bottom plot shows the absolute deviations of the approximate lines from the original ones.

Finally, we should notice that the previous implementation[9,18,23] of the M-I model has not taken into account lateral interactions between the adsorbates. According to our approach, for example, the binding energy of O2 shifts (Figure 2) by 0.07–0.20 eV, depending on the coverage of neighboring sites by oxygen atoms. According to the evaluation provided above, such an energy shift results in a significant dependence of the rates of elementary reactions on surface coverage by adsorbates. Indeed, the explicit account of lateral interactions in M-II results in a visible alteration of the calculated O2 TPD peak (Figure 4). Instead of the symmetric peak produced[18] by the M-I model, the M-II model leads to a clearly asymmetric profile, formed by several overlapping peaks of different sizes, corresponding to oxygen desorption from nearby energy levels, determined by the lateral interactions. The resulting peak shape is much closer to the experimentally measured[22] one. We have further improved the match of our calculated profile with the experimental one by the application of smoothingon the simulated oxygen desorption signal f(t′), where t′ is the system’s time, and τ is a smoothing constant. Underlying this formula is the fact that in an experimental setup, desorption spectra are recorded as the amount of desorbed molecules seen by a spectrometer’s detector at each moment of time. It is natural to assume a finite delay between the actual desorption of a molecule out from the sample’s surface and its registration by the detector. Avoiding a complicated account for exact geometry of an experiment, it is natural to assume that the mentioned delays between desorption and registration of molecules are distributed exponentially, effectively leading to smoothing 3 of the spectra. The resulting TPD profile (Figure 4, red line), calculated for τ = 4.0 s, matches the experimentally measured one almost perfectly (dependence of the result on the choice of τ is further explored in Figure S1 in the Supporting Information). We consider that such smoothing efficiently imitates the response curve of a spectrometer used in TPD experiments.

Conclusions

In conclusion, we have investigated some features of a widely applied model (M-I) of oxygen interaction with so-called cus sites of RuO2(110) surface. Although a number of works report successful application of the M-I model (single step O2 adsorption/desorption) for calculations of O2 temperature desorption profiles and for investigation of various chemical processes on the RuO2(110) surface without visible problems, we have found significant hidden differences in the oxygen desorption rates predicted by M-I and our alternative, more complex, model, M-II. These differences emerge from (i) oversimplification of the oxygen adsorption process in M-I; (ii) the lack of lateral interactions of the adsorbates on the surface; (iii) the assumptions on equality of partition functions of adsorbed atoms and molecules to unity; and finally, (iv) intrisic errors associated with the accuracy of DFT. Our results for the very stringent test of O2 TPD, obtained with the use of the M-II model, show a visible improvement in the shape of the simulated desorption peak. We also expect that improvements of M-II compared with M-I are crucial for simulations of complex chemical processes on RuO2(110) surfaces, especially at temperatures above 500 K, where most of the interesting oxidation chemistry appears.[12]

We firmly believe that the analysis presented paves the way for a quantitave use of kinetic Monte Carlo tools for systems showing high anisotropy at the nanoscale.