Single-atom catalysis (SACs) [1] has become the most active new frontier in heterogeneous catalysis field in recent years [2-6]. The most important structure feature of SACs contributing to catalytic activity and selectivity is the isolated atoms dispersed on the support, which could be analogical to homogeneous catalysts [7]. In general, the high ratio of exposed individual atoms is possible to provide uniform active sites and maximize the catalytic efficiency compared to heterogeneous catalysts consisting of various active sites [8]. With the great advances in synthetic methodologies, a large number of single-atom catalysis (SACs) has exhibited excellent performances for various chemical reactions [9-12]. Interestingly, when we gaze at the typical metal oxide surface, it could be found that the surface is alike the combination of a series of single metal atom center coordinated with oxygen anions. Moreover, it is worth noting that for the SACs, the bonding capacity of single atom may limit its application in the complicated reaction.

The catalytic combustion of methane (CH₄) [13-15] is a complex process which should include four steps of dehydrogenation and two C-O coupling elementary steps at least. It is a vital alternative to the conventional thermal combustion processes of coal and oil [16] for energy production. However, as a major component of natural gas, CH₄ is also a greenhouse gas, which has as much as 28 times more potent than carbon dioxide (CO₂) [17]. Accordingly, its highly efficient combustion could relieve the greenhouse effect due to exhaust emission [18-20]. To meet the lower level of methane emission requirement in stricter environmental standards such as the Euro-VI standard for natural gas vehicles (NGV), there are urgent requirements for high activity catalysts to maintain high combustion rates for the abatement of methane emissions in the NGV exhaust [21]. In order to achieve high activity and thermal stability for catalytic combustion of methane, extensive studies have been done and various catalysts have been developed by scientists, such as supported palladium oxides [22] and rhodium oxides [23] catalysts. Cargnello et al. recently reported a modular Pd@CeO₂ nanoparticles on functionalized Al₂O₃ that could achieve both catalytic activity below 400 ℃ and thermal stability up to 850 ℃ [24]. Nevertheless, due to the high cost, sintering and susceptibility to poisoning of noble metal catalysts, cobalt-based spinel oxides, e.g. Co₃O₄, which are treated as alternatives also exhibit excellent catalytic performance for methane combustion [25-28]. Recently, Hu et al. [29] investigated a complete catalytic cycle for methane combustion on the Co₃O₄(110) surface based on first-principles calculations. It is found that the cooperation of multiple sites on Co₃O₄ surface may play an important role in promoting methane combustion over Co₃O₄ surfaces. Yet the effect of multi-site on CH₄ combustion activity just lies on the energetic analysis from density functional theory (DFT) calculations. Still, it is not easy to intuitively comprehend the importance of the cooperation between multiple active sites for the overall catalytic activity compared to SACs-like separated single active site.

Kinetics simulation as an efficient approach linking macroscopic and microscopic properties of catalytic systems [30-33] could be used to analyse complex catalytic system quantitatively. Over the last few years, kinetic Monte Carlo (kMC) simulation based on DFT calculations has been widely used in heterogeneous catalytic reactions and biomolecular reactions [34-42]. As kMC simulations intrinsically take all the location information and diffusion issues into account, it is an effective way to describe the single- or multi-site effect of heterogeneous catalysis quantitatively.

In this paper, to quantitatively understand the kinetic influence of the interplay of multi site for the catalysis, the reaction mechanism through multi active sites is systematically compared with that constrained at the single active site for the CH₄ combustion over Co₃O₄(110) utilizing first-principle kMC simulation. First, we explore a complete CH₄ combustion cycle on an ideal isolated single active site on Co₃O₄(110) surface. Based on the energetic data from DFT calculations, we then perform kMC simulations for both two types of catalytic cycles of CH₄ combustion on Co₃O₄(110) respectively through single and multiple active sites. By the virtue of kinetic calculations, we manage to quantitatively analyze the effect of the interplay of multiple active sites for the activity and selectivity of CH₄ oxidation. Furthermore, the key factors governing the kinetic activities are also disclosed, which may help to improve the catalytic activity of spinel-type cobalt oxide catalysts for CH₄ oxidation.

All the spin-polarization calculations in this work were performed using the Perdew-Burke-Ernzerh (PBE) generalized gradient approximation (GGA) functional under the framework of density functional theory (DFT) in Vienna ab initio Simulation Package (VASP) [43-46] with the employment of a plane wave basis set [47]. The projector-augmented-wave (PAW) pseudopotentials [48] with the cores of 1s 2p, 1s, and 1s for Co, O, and C, respectively, were utilized to describe the valence-core interactions. A plane-wave kinetic energy cutoff of 500 eV was set for the expanded plane wave basis-set. As for the strong correlation between the electrons in partially occupied 3d orbitals of Co, the on-site Coulomb repulsion correction term of U within the Hubbard scheme (PBE+U) was applied in the 3d electrons of Co. The effective U_eff value, namely U−J, was set to 2.0 eV, which has been demonstrated to well simulate the properties of Co₃O₄ such as the band gap, lattice constant, and magnetic moment in previous work [28, 49-52]. Van der Waals interactions were determined by the BJ-damped variants of the Grimme's D3 approach [53, 54]. The Brillouin-zone integration was performed on Monkhorst-Pack grids with a 2 × 3 × 1 mesh where a Gaussian-smearing approach with σ = 0.05 eV is used during the ionic optimization.

The Co₃O₄(110) surface was modeled using a stoichiometric slab with of a p(2×1) unit cell, which was based on the antiferromagnetic bulk Co₃O₄ with normal spinel structure (a = 8.124 Å). Each periodic slab contains eight layers of Co-O plane and separated by a vacuum space of approximately 15 Å. Co₃O₄ surface may have two different terminations [55-57]. We investigated the stabilities of the two terminations of the Co₃O₄(110) surface in our previous work [29]. The (110)-A termination exposes two types of cobalt cations (Co²⁺ and Co³⁺) and one type of oxygen anion (3-fold coordinated oxygen anion, O_3c), whereas the (110)-B termination has only one type of cobalt cation (Co³⁺) and two types of oxygen anions (2-fold coordinated oxygen anion O_2c and 3-fold coordinated oxygen anion O_3c). The (110)-B termination which has a lower Gibbs surface energy under the work conditions was utilized in the work. For the spin configuration of the slab, we used the anti-ferromagnetic configuration from the subsurface to the bottom layer of the slab model. It is consistent with the bulk properties. For the surface layer, the spin of the unsaturated Co_4c was opposite to the spin of Co²⁺ in the subsurface. The same spin configuration has also been confirmed as the most stable configuration for the Co₃O₄(110) in previous work [58]. Top and side views of the Co₃O₄(110) surface structure is displayed in Fig. 1(a). Modeling details can also be found in our previous work [29]. In this work, we treat a O_3c-Co³⁺-O_2c as a single active site on Co₃O₄(110) surface as labeled in Fig. 1(a) by black dashed square.

All the geometry structures of adsorption intermediates were optimized using a force-based conjugate gradient algorithm until the forces on all the relaxed atoms were below 0.05 eV/Å. During the geometry optimization, the bottom four Co-O layers were fixed while the top four Co-O layers and adsorbates could be relaxed. The transition states (TS) in reactions were searched with a constrained optimization scheme [59] and the convergence of forces was set to 0.05 eV/ Å. It was achieved when all the forces on atoms vanished and the total energy was a maximum along the reaction coordination but a minimum with respect to the rest of the degrees of freedom (namely saddle point). Each TS was further verified as a first-order saddle point with only one imaginary vibrational frequency and the corresponding vibrational mode along the reaction coordination on the basis of a numerical vibrational frequency analysis.

For the kinetic Monte Carlo simulation of CH₄ combustion process on Co₃O₄(110), we implemented a simulation code based on the VSSM kMC method [34, 60-62]. To overcome the stiff problem introduced by the fast species diffusion on surface and accurately describe its influence on the kinetics, the fast species redistribution (FSR) algorithm [63] was adapted in our code where the fast diffusion species are redistributed within feasible space under quasi-equilibrium assumption during kMC iteration. For the surface grid modeling of Co₃O₄(110), Co, O_3c, O_2c sites and O vacancy are all explicitly modeled in a 30 × 30 two-dimensional grid with periodic boundary condition enabled. In simulation process, surface adsorbate coverages and turnover frequencies (TOF) were calculated from statistical time averages [39]. The steady state was detected based on exponentially weighted moving average (EWMA) protocol [64-66] where the tolerance range parameter L was set to 6 and the decay parameter λ of EWMA was 0.05. The convergence criterion p_c = 0.05 was choosen to ensure enough kMC sampling in our work. The degrees of rate control (X_RC) for elementary reaction i were calculated directly using the formula in the definition [30, 32, 67] based on numerical differentiation:

where r refers to the reaction rate, the partial derivative is obtained through keeping the rate constant, k_i, for the other steps (j≠i) and each equilibrium constant, K_i as constants. All kinetic Monte Carlo simulations were implemented upon free energy landscape [68, 69]. The free energy was calculated based on the total energy from DFT calculations with the thermodynamic correction. The calculation process is as follows:

where G, E and C_P refer to the chemical potential (partial molar Gibbs free energy), electronic energy and heat capacity, respectively. The entropy term can be expressed as the sum of the translational, rotational, vibrational and electronic contributions as to:

And finally, the intrinsic zero-point energy (ZPE) corrections can be included to finally obtain:

For the case of solids and adsorbates, some approximations can be assumed: (1) S_e ≈ 0 at the fundamental electronic level; (2) Translational and rotational motions can be neglected, therefore, S_t ≈ 0 and S_r ≈ 0. In this sense, all the entropy contributions come from vibrations: S = S_v. Similarly, translational and rotational contributions to the heat capacity are neglected. Therefore, Gibbs free energies for the different states have been calculated as to:

In this paper, all energetic data reported are free energy after correction under the reaction conditions of T = 450 ℃, P(CH₄) = 0.01 bar, P(O₂) = 0.2 bar, P(CO₂) = 0.01 bar, and P(H₂O) = 0.02 bar. The original DFT calculated total energies and details of kinetic parameters calculation methods are shown in the Supporting Information.

The complete catalytic cycle and possible reaction network of CH₄ combustion on the Co₃O₄(110) surface has been explored via multi-site mechanism in the total energy landscape in our previous work [29]. The optimal reaction route of CH₄ complete oxidation following the optimal multi-site mechanism can be summarized as: CH₄ → CH₃(Co) → CH₃O_2c → CH₂O_2c → O_2cCH₂O_2c → O_2cCHO_2c → O_2cCHO(Co) → CO₂. Each adsorbate bonded with Co³⁺ is denoted with (Co) throughout this work. The CH₄ deep oxidation relies on the H and the oxidative dehydrogenation intermediate migration between multi active O_2c sites. Methane undergoes four dehydrogenation and two C-O coupling steps. Gaseous O₂ could be dissociated at lattice oxygen vacancy to fill the oxygen vacancy generated by the CO₂ desorption and to generate O(Co) which could help to form H₂O on the surface to achieve active site regeneration [29]. Herein the calculated free energy profiles of the optimal reaction pathways of the methane catalytic oxidation at multi sites under the reaction condition are illustrated in Fig. 2(a). Based on previous PBE+U calculation results, we further investigated the catalytic cycle assuming it is confined at the isolated single active site (O_3c-Co³⁺-O_2c).

The calculated free energy profiles of the optimal reaction pathways of the methane catalytic oxidation at multi and single site are displayed in Fig. 2(a) and Fig. 2(b), respectively. Similar with the cycles on multi sites, the first C-H activation step on single active site still has a high free energy barrier of 1.80 eV (the TS structure is shown in Fig. 1(b)). With the first C–H bond being dissociated, the CH₃ and H fragments are bonded at Co and O_2c sites. Different from the reaction on multiple active sites, on the isolated active sites (O_3c-Co³⁺-O_2c), the H atom is assumed that it could not be removed via the migration between two O_2c sites to recover the active O_2c site. Thus, the CH₃(Co) on Co³⁺ has to be coupled with the O_3c to form CH₃O_3c with free barrier of 1.40 eV (Fig. 1(d)) since the energy barrier of H₃C-O_2cH coupling is even as high as 2.39 eV, which is much higher than that of C-O_2c coupling (0.82 eV, Fig. 1(c)). Although the C-O_3c coupling barrier is lower than that of first C-H breaking, the endothermic reaction of the first step leads to a higher potential energy of the transition state of C-O_3c coupling. Hence, it is possible that the C-O coupling may significantly reduce the overall activity of CH₄ combustion.

Besides the high activation barrier in the first C-H breaking step, the process of CH₂O(Co) generation from (O_3c)CH₂O(Co) also needs to overcome a high free barrier of 1.30 eV, indicating that this step could also be a key elementary step. On multiple sites, the CH₂ structure could adsorb across two O_2c sites to form more stable O_2cCH₂O_2c followed by the oxidative dehydrogenation by adsorbed O₂. However, when the reaction is confined to the single active site, (O_3c)CH₂O(Co) could not be dehydrogenated directly but be converted to CH₂O(Co) via C-O breaking to break the spatial limitation to enable the further dehydrogenation towards CO₂ over Co₃O₄(110). We also find that CH₂O(Co) seems to be desorbed easily to generate HCHO releasing the energy of –1.05 eV. Comparing the further dehydrogenation of CH₂O(Co) on single site (∆^‡G = 1.21 eV, ∆G = –1.04 eV) and direct desorption of O_2cCH₂O_2c on multi sites (∆G = 0.31 eV), it is more possible that CH₄ would be oxidized selectively to HCHO instead of the complete oxidation to CO₂ over the surface when the reaction cycle is confined on single active site.

Furthermore, similar to multi-site mechanism, single-site mechanism also depends on the dissociation adsorption of O₂ at the oxygen vacancy to generate active O(Co), which facilitates the left H at lattice oxygen after CH₄ oxidation to be removed. All the calculated free energy barriers and free reaction energies of the elementary steps involved in the catalytic cycle of CH₄ combustion on Co₃O₄(110) following single-site and multi-site mechanisms are listed in Table 1 and Table 2, respectively. The structures of transition states and intermediates in CH₄ oxidation cycle on limited single active site are displayed in Fig. 3.

Hence, the optimal reaction route of CH₄ complete oxidation at the Co₃O₄(110) single active site can be summarized as: CH₄ → CH₃(Co) → CH₃O_3c → CH₂O_3c + OOH(Co) → O_3cCH₂O(Co) → CH₂O(Co) → CHO(Co) → O_3cCHO → CO₂. In the reaction pathway, methane undergoes four dehydrogenation reactions, three coupling processes and one C–O bond-breaking reaction. Compared with reactions at multiple active sites [29], the methane oxidation process is required to undergo more elementary reactions.

As illustrated in the energy profile of complete catalytic cycle on single active site (in Fig. 2(b)), we can find that the first C-O coupling step has the highest effective free energy barriers (2.26 eV). This indicates that the first C-O coupling step may be the key step with the greatest impact on the entire cycle. However, when analyzing the elementary reactions separately, we can find that the maximum free energy barrier is 1.80 eV for the first C-H breaking step. Therefore, it is difficult to quantitatively analyze the properties of the whole catalytic cycle from the energy profile alone. Furthermore, so far we have obtained both catalytic pictures about methane combustion cycles on both single and multiple active sites over Co₃O₄(110) based on PBE+U calculations. Nevertheless, it is still not enough for quantitative comparison of the effects of single and multiple sites to catalytic oxidation activities. In the following parts, we will perform kinetic analysis on different catalytic cycles using kinetic Monte Carlo (kMC) simulation for the further analysis.

In this section, we performed kinetic analysis utilizing kMC simulations for both CH₄ combustion catalytic cycles on single active sites and multi active sites of Co₃O₄(110) surface. The calculated free energetic data used for kMC simulations for single-site and multiple-site cases are listed in Table 1 and Table 2, respectively. All kMC simulations were performed at the typical temperature of complete conversion of the methane (T = 450 ℃, P(CH₄) = 0.01 bar, P(O₂) = 0.2 bar, P(CO₂) = 0.01 bar, P(H₂O) = 0.02 bar). Table 3 shows the simulated steady-state coverages of different adsorbates on three different sites (O_3c, Co, and O_2c) for CH₄ combustion cycle on multi active sites. Apart from the unoccupied sites, the surface H(O_2c) is the dominant species that has the maximum coverage (6.50 × 10^–2 ML) on the surface. The coverage of CH₃ species is secondly to H(O_2c) (1.61 × 10^–3 ML) on O_2c site, which is consistent with the observations in experiments [28]. From the surface coverage results, it is clear from Table 3 that the C containing species are mainly located at the Co and O_2c site. Therefore, compared with the O_3c site, the active O_2c site plays a more important role in the CH₄ complete oxidation with the help of cooperation of multiple sites on Co₃O₄(110) surface. The simulated turnover frequency (TOF) towards CO₂ is 5.57 × 10^–3 s^–1 following the optimal mechanism. HCHO is a possible by-product from the surface reaction. It has a lower yield of 7.18 × 10^–4 s^–1. However, since 450 ℃ is above the ignition point of HCHO (300 ℃), the generated HCHO would be rapidly burned to completely form CO₂ as well. In order to identify the rate-determining steps, we calculated event frequencies for each elementary reaction in the CH₄ combustion cycle and the degree of rate control (X_{RC, i}) of elementary step i [30, 32] as displayed in Fig. 4 using kinetic Monte Carlo simulation. It is clear that most of the elementary steps reach the quasi-equilibrium, namely the reversibility of forward and reverse step is close to 1, except for the first C-H bond breaking, the first C-O coupling, the second C-O coupling and HCHO desorption. The further degree of rate control analysis discloses that the first C-H bond activation step possesses the greatest X_RC value, indicating that the first C-H activation is the rate-determining step for the catalytic cycle of CH₄ combustion on multi active sites.

To intuitively comprehend the effects of multi active sites on the catalytic activities, we also calculated kinetic properties of CH₄ combustion cycle confined at the single site of Co₃O₄(110) surface as the reference. The simulated surface coverages of different adsorbates on three different sites are listed in Table 4. In the single active site case, as the active O_2c site is passivated by H species, the major active sites for CH₄ combustion turn to Co and O_3c. As shown in Fig. 5, different from the multi-site case, only HCHO is generated with a low rate of 3.98 × 10^–6 s^–1, indicating that HCHO rather than CO₂ is the main product for CH₄ combustion over the Co₃O₄(110) surface if the reaction cycle is limited on the single active site. Although gaseous HCHO is also rapidly converted to CO₂ under the reaction conditions, the turnover frequency is still much lower than that following the multi-site mechanism. This suggests that the single active site has extremely low activity for oxidation of CH₄. Moreover, interestingly, forward and reverse elementary event frequencies in Fig. 5 show that the first C-H activation step is almost equilibrated following the single-site mechanism. This indicates that it may be not the rate-limiting step for CH₄ combustion cycle anymore if following this mechanism. This is proved by the further degree of rate control results. The X_RC value of the first C-H activation step is 0.11 which is significantly less than 0.97 of the first C-O coupling. This indicates that the first C-O coupling instead of the first C-H breaking step becomes the greatest limitation to the overall catalytic oxidation activity if following single-site mechanism. Hence, it is demonstrated from the reverse side that the cooperation of multi sites could effectively ease the deep dehydrogenation and C-O coupling for methane combustion over Co₃O₄(110) but is unable to facilitate the first C-H bond activation.

By virtue of kMC simulation, the effect of active sites and kinetic behaviors of CH₄ complete oxidation on the Co₃O₄(110) surface can be analyzed quantitatively. By comparing the final overall turnover frequencies of different cycles, it can be found that the CH₄ catalytic combustion rate following the multi-site mechanism would be more than 3 orders of magnitude higher than that confined at the separated single site. It indicates the synergistic effect of multi active sites could promote the catalytic activity substantially. From the statistical results of surface coverage and adsorbate distribution on different sites (O_3c, O_2c, and Co) for two different cases, it is found that O_2c and Co are the major active sites of multi sites catalytic cycle, while O_3c and Co are main reaction sites for single site oxidation cycle. It could be attributed to the fact that the multi sites could assist the rapid transfer of H species on the Co₃O₄(110) surface to prevent the passivation of active two-coordinated O sites (O_2c). However, on the isolated single active site, once the O_2c is occupied by H and blocked, the less active O_3c has to be the reaction site. Consequently, it is the CH₃* and O_3c coupling that becomes the rate-determining step for the methane combustion. This significantly limits the whole reaction activity. It suggests that the fast removal of H on O_2c by the assistance of multi sites could switch the C-O_3c coupling to C-O_2c coupling to enhance the activity of CH₄ oxidation. In addition to species migration, the space limitation broken by multi sites could also provide more possibilities for the surface intermediates to reach more stable structures. With the assistance of two neighboring O_2c sites, the CH₂O_2c could be reconstructed to more stable O_2cCH₂O_2c structure where C atom switches from sp² hybridization to sp³ hybridization. The C-O bonding could be therefore changed from the double bond to single bond, promoting the surface Co-O bonding and consequently promoting the intermediate stability. This kind of intermediate stabilization allows the reaction to proceed towards the complete combustion to CO₂ rather than the selective oxidation to HCHO. Vice versa, the CH₄ confined at the single active site can only be oxidized to HCHO selectively. Therefore, from the quantitative kinetics comparisons between CH₄ complete oxidation cycles on single site and multi sites, we can find that the cooperation of multi active sites could not only avoid rapid passivation of single sites but also facilitate to form more stable intermediate to promote methane complete combustion substantially. In other words, it also indicates that it is possible to efficiently regulate the selectivity of CH₄ oxidations by confining the active sites to form separate sites. Nevertheless, the kMC simulation results also disclose that the conversion of methane is possible to be weakened due to the confined and isolated active sites.

In summary, based on the quantitative comparative results from the first-principles kinetic Monte Carlo simulation on single active sites and multi active sites, we find that the cooperation of multi active sites in Co₃O₄(110) could substantially promote the activity of the methane oxidation compared to single active sites. The selectivity of the surface oxidation would be changed as well. The complete oxidation product of CO₂ could be yielded following the multi-site mechanism while the gaseous HCHO would be produced following the single-site mechanism over the Co₃O₄(110) surface. The synergistic effect of multi active sites mainly originates from: (1) providing space for species migration to avoid rapid passivation of single sites, and (2) providing more possibilities for surface intermediates to reach more stable structures. The confinement of separated single active sites can only selectively catalyze CH₄ oxidation to generate gaseous HCHO if not considering the rapid direct combustion of HCHO to CO₂ in the gas phase under the reaction conditions. Consequently, from the other side, it sheds light on the fact it is of great significance to confine active sites in the catalyst surfaces for the selective oxidation of CH₄.