Plasma-assisted CH₄ activation on Cu/CeO₂ catalysts: insights into the effect of catalyst surface and vibrational excitation

Abstract

The lack of chemical understanding and efficient catalysts impedes the development of plasma–catalytic CH₄ conversion. In this work, we employ density functional theory calculations to understand the effects of vibrational excitation on the dissociative chemisorption of both CH₄ and CH₃ on surfaces relevant for (plasma-assisted) catalysis, i.e., Cu(111), CeO₂(111), and a single Cu atom supported on CeO₂(111). The single-atom Cu catalyst (Cu₁/CeO₂(111)) shows the lowest energy barrier (0.35 eV) for CH₄ dissociation among the three surfaces. The vibrational mode-specific reactivity of CH₄ and CH₃ is assessed using the sudden vector projection (SVP) model, in which the stretching mode of CH₄ is dominant for CH₄ dissociation on these three surfaces. Additionally, depending on the reaction mechanism of CH₃ chemisorption and dissociation, either the stretching or bending modes are predicted to be more effective at promoting reactivity. Furthermore, vibrational efficacies for dissociative chemisorption of CH₄ on the investigated catalyst surfaces are compared using a simplified model, which also employs SVP calculations, to reveal the importance of mode specificity and the structural dependence of the catalyst, offering valuable insights into catalyst design in heterogeneous and plasma catalysis.

---

Introduction

Methane (CH₄) is a very stable molecule. Dissociative chemisorption of CH₄ on catalyst surfaces is typically considered the rate-controlling step in CH₄ conversion.^1–3 Therefore, understanding and predicting C–H bond activation are key to comparing catalytic performances across various catalysts and for guiding the rational design of effective catalysts.⁴ Plasma, a partly ionized gas, is an alternative method to drive reactions, viz. with electrical energy instead of thermal energy. When coupled to catalytic materials, plasma catalysis offers a means to drive chemical reactions under ambient conditions, without the need for elevated temperatures or pressures.^5–11 Typically, the function of a catalyst is to reduce the barrier height for the rate-controlling reaction or to regulate the product distribution. This clear-cut function could also apply to plasma catalysis.

Specifically, reactants can be activated in the plasma before interaction with the catalyst takes place.⁵ Hence, energy can be converted from different internal degrees of freedom (DOFs) of the reactant into the reaction coordinate, effectively increasing the amount of energy available for the reaction. The vibrational modes of the reactants can be selectively activated by the reduced electric field (i.e., the ratio of the electric field to the gas number density) in non-equilibrium plasmas.^6–8 The energy stored in specific vibrational modes can increase the vibrational temperature (T_vib > 1000 K), which may promote energy-efficient gas conversion in a plasma or on a catalyst.^6–8 For example, plasma-induced excitations can play an important role in CH₄/CO₂ conversion and NH₃ synthesis by plasma catalysis, as shown in both experimental and theoretical studies.^12–17 Nozaki et al. studied the role of vibrational excitation of CH₄ in plasma catalysis by emission spectroscopy and observed that the vibrational temperature of excited species largely increased with the packing bed temperature, when using Ni/SiO₂ as a catalyst.¹² Furthermore, the importance of the bending mode of vibrationally excited CO₂ molecules on a Pd₂Ga/SiO₂ alloy catalyst system was investigated by Kim et al.¹³ These authors found that a lower effective barrier height, due to vibrational excitation, can significantly improve the reaction performance, in which the CO₂ conversion is not only increased more than a factor of two but also breaks the thermodynamic equilibrium limitation compared to thermal conditions.¹³ Mehta et al. combined microkinetic modelling with experiments to investigate plasma–catalytic NH₃ synthesis from N₂ and H₂. They suggested that vibrationally excited N₂ can lower the barrier height in the catalyst–adsorption step and thereby increase the catalyst activity for metals that bind N₂ weakly in thermal catalysis.¹⁴ In a follow-up study, Engelmann et al. specifically investigated the effect of various N₂ vibrational energy distribution functions, including those characteristics for dielectric barrier discharge plasmas commonly used in plasma catalysis. They also compared the role of vibrationally excited N₂ molecules vs. plasma-generated radicals.¹⁵ In another study, Engelmann et al. also investigated the reactions of vibrationally excited CH₄ and hydrocarbon radicals on transition metal catalysts, revealing that vibrational excitations and plasma-generated radicals can impact the reaction rates and product selectivity.¹⁶ Finally, Michiels et al. also revealed by microkinetic modeling that vibrationally excited CO₂, as well as radicals, can increase the turnover frequency of CH₃OH formation in plasma–catalytic CO₂ hydrogenation.¹⁷ It should be noted that in the aforementioned modelling studies, the plasma-induced vibrational excitations for the gas–surfaces reactions were always described by the Fridman–Macheret (F–M) α model, originally formulated for the gas-phase reaction of an atom with a diatomic molecule.¹⁸ Unfortunately, this model seems to be unable to capture the complexity of molecule–metal surface interactions arising from high-dimensional potential energy surfaces (PESs), which cause dynamical and state-specific effects in plasma catalysis.^19–21

Quantum dynamics studies can provide high precision by incorporating all degrees of freedom and quantum effects. A lot of effort has been devoted to constructing a PES to describe the non-statistical nature of the dissociative chemisorption of CH₄.²² In order to avoid the extreme, intractable computational cost of full-dimensional quantum dynamics, some approximate methods have been proposed, e.g., quasiclassical trajectory (QCT), and the reaction path Hamiltonian (RPH) approach.²² In the QCT approach, the quantum mechanical vibrational energy is imparted to each mode in the initial velocity setup by treating classical objects following Newton's laws of motion.²³ QCT is fast and intuitive for large systems but fails to capture quantum tunneling or zero-point energy effects, leading to inaccuracies at low energies. Recently, Gerrits et al. proposed a new approach to more accurately predict the reactivity under catalytically relevant conditions (i.e., low translational and high vibrational energies in the molecule) by extending a ring polymer molecular dynamics (RPMD) approach to include surface atom motion.³ In contrast to the QCT approach, RPMD can include nuclear quantum effects, like tunneling, and remedy the artificial zero-point energy leakage of the molecule into the reaction coordinate for translational energies below the minimum barrier height. This RPMD method offers more accurate predictions of the experimental sticking probabilities (i.e., a measure of reactivity), which could also be employed for non-equilibrium conditions in, e.g., plasma catalysis.³

Besides, Bal et al. developed an indirect approach by implementing a bias potential in molecular dynamics (MD) simulations, in which the selected mode can be excited to higher temperatures, while all others remain at thermal equilibrium.²¹ Furthermore, different Ni surfaces and vibrational modes were investigated using the same approach to understand the impact of vibrational non-equilibrium on the dissociative chemisorption barrier of CH₄.²⁴ The effect of vibrational excitation on the free energy barrier was predicted to be larger on terrace sites than on surface steps. Also, even at a low vibrational temperature, high vibrational efficacies (i.e., the quantitative effect of vibrational energy on the reactivity compared to translational energy) were obtained. The efficacy of the symmetric stretch was greater than that of asymmetric stretches, which, in turn, was higher than that of the bending modes, in agreement with experiments.²⁴

The RPH approach was also proposed by focusing on the minimum energy reaction path rather than the full-dimensional PES to understand the dynamics of vibrational mode-specific chemistry.^25–27 Jackson et al. used this approach to understand the dissociative chemisorption of CH₄ on a Ni(100) surface based on a harmonic expansion of the vibrational modes along the reaction path. Among the vibrational modes, the symmetric stretch (ν₁) exhibited the highest efficacy, as it strongly couples to the reaction coordinate and softens at the transition state (TS).²⁶ Likewise, Roy et al. reported the same trend that the symmetric stretching vibrational mode of CH₄ on the Ni/Pt-bimetallic alloy exhibits the highest reactivity due to significant mode softening near the TS, lowering the effective barrier height.²⁷

Indeed, dynamical simulations can provide a better understanding of the dynamical effects, including energy transfer between the molecule and the catalyst surface, thermal local barrier height modulation, and the bobsled effect.^28–34 Interestingly, Jiang et al. proposed a sudden vector projection (SVP) model to qualitatively predict the vibrational state-specific efficacies in reactions involving polyatomic reactants, without requiring complex dynamical simulations.^35–39 The SVP model assumes the energy transfer as occurring on a timescale too short for full intramolecular vibrational redistribution, which has been demonstrated to be valid in several gas phase reactions and gas–surface interactions.^38–41 Jiang et al. compared the computed efficacies of vibrationally excited CH₄ in the reaction on a rigid Ni(111) surface between MD simulations and the SVP model.³⁵ In the SVP approach, the two stretching modes exhibited higher vibrational efficacy than the bending modes, in agreement with both theoretical and experimental observations.^35,42,43

Recently, Gerrits and Bogaerts proposed an improved model to capture the effect of vibrational excitation on dissociative chemisorption rates in heterogeneous and plasma catalysis.¹⁹ Unlike the widely used F–M α model—which lacks vibrational mode specificity and very poorly matches MD and experimental molecular beam data—the new model combines the forward barrier height, “lateness” of the transition state geometry (ratio of the dissociating and equilibrium bond lengths), and the overlap of vibrational modes with the reaction coordinate using SVP calculations. This approach yields significantly better agreement with MD-derived vibrational efficacies (R² = 0.52 vs. −0.35) at a similar computational cost, making it more suitable for microkinetic modeling of vibrationally excited molecule–metal surface reactions.¹⁹

The key to CH₄ conversion lies in finding efficient catalytic systems with a controllable reaction kinetics process. Single-atom catalysts (SACs) have attracted wide attention as promising candidates for tackling challenges in energy conversion, environmental remediation, and chemical synthesis, owing to their maximum metal dispersion, precise control over catalytic sites, and enhanced reactivity and selectivity.^44,45 Since the typical mean electron energy in a CH₄ plasma is in the range of 1–5 eV, vibrational excitation of CH₄ due to impacting electrons is assumed to be more prevalent than depositing the energy in other channels. The vibrational excitation of CH₄, in turn, can enhance chemical reactivity of the molecule on a catalyst surface, compared to thermal (i.e., heterogeneous) catalysis.^46–48 Thus, in this paper, we use density functional theory (DFT) calculations to investigate the performance of CH₄ activation on a SAC, by comparing three typical surfaces, i.e., Cu(111), CeO₂(111), and a single Cu atom supported on CeO₂(111), denoted as Cu₁/CeO₂(111). The effects of vibrational mode-specificity of CH₄ on the three surfaces are investigated to understand CH₄ activation in plasma catalysis. We note that the F–M α model is currently used as a one-fits-all tool for quantitative prediction of the effect of vibrational excitation in microkinetic models. However, even for the dissociative chemisorption of diatomic molecules, it is inaccurate, where obviously no mode specificity is present.²⁰ In this work, the qualitative comparison between the F–M α model and the alternative η model is also discussed, to indicate the importance of mode-specificity and the relationship between the structural dependence of the catalyst and vibrational efficacy.

Finally, we will also discuss the further dissociation of a CH₃ radical from the plasma on the aforementioned three surfaces, because besides vibrational excitation, microkinetic models revealed that radical chemistry might also be important for plasma catalysis.^11,15–17

Methodology

DFT calculations

Spin-polarized DFT calculations were performed using the Vienna ab initio simulation package (VASP, version 6.2.1).^49–51 The Perdew–Burke–Ernzerhof (PBE) exchange–correlation density functional⁵² was employed and the projected augmented wave method was used to describe the core electrons.^53,54 A cutoff energy of 500 eV was adopted throughout. To adequately describe the electron localization in the Ce 4f orbital, an on-site Coulomb repulsion was applied, as described by the Hubbard U parameter.⁵⁵ The U value was set to 4.5 eV for the Ce 4f orbitals, providing sufficient localization of the electrons left on CeO₂.^56–58 The maximum force criterion for ionic convergence was set to 0.02 eV Å⁻¹. van der Waals dispersion forces between adsorbates and surfaces were accounted for using the zero damping DFT-D3 method of Grimme.⁵⁹ All configurations were optimized using the conjugate gradient algorithm. The climbing image nudged elastic band (CI-NEB) method and minimum-mode following dimer methods were adopted to determine the TS structures of elementary reaction steps as implemented in the VASP Transition State Tools package.^60,61 All TSs were validated by vibrational frequency analysis.

The projected crystal orbital Hamilton population (COHP) curves were calculated using LOBSTER to analyze the bonding and anti-bonding states,^62,63 where the pbeVASPFit2015 basis set was used for the projection of wave functions. Bader charges were calculated for the electron population analysis.⁶⁴

The adsorption energy, E_ads, is defined as:

E_ads = E_{adsorbate+surface} − (E_surface + E_adsorbate) (1)

Here, E_{adsorbate+surface}, E_surface, and E_adsorbate are the total energies of the adsorbate on the slab, the clean slab, and the gaseous adsorbate, respectively.

Surface structures

The close-packed (111) surfaces of Cu and CeO₂ are studied because they represent the most stable terminations observed experimentally. Furthermore, the CeO₂(111) surface in particular has been widely used as a support for various metals in many theoretical studies.^57,65,66 In addition, we also investigate a single Cu atom on CeO₂(111), i.e., Cu₁/CeO₂(111), because single-atom catalysts are gaining attention for their maximal metal dispersion, precise site control, and improved reactivity and selectivity in energy, environmental, and chemical applications.^44,45 These three specific surfaces help us theoretically understand the CH₄ dissociation mechanism and its vibrational excitation due to their distinct TS structures. As shown in Fig. 1A, the Cu(111) surface was represented by a four-layer slab and a (5 × 5) supercell, with a vacuum gap of 15 Å. The bottom two layers were fixed at the ideal bulk positions, while the positions of all other atoms were allowed to relax. The surface model of CeO₂(111) was represented by three O–Ce–O triple-layers with (4 × 4) supercells, again with a vacuum layer of 15 Å (Fig. 1B). The bottom two O–Ce–O layers were fixed, while the top O–Ce–O layers were allowed to relax. For a Cu atom adsorbed on CeO₂(111), the hollow site (Fig. 1C) is the most stable adsorption site.⁶⁷ A (2 × 2 × 1) γ-centered k-point mesh was used to sample the Brillouin zone for the Cu(111) surfaces, with a Methfessel–Paxton smearing width of 0.2 eV. For the surfaces of CeO₂(111) and Cu₁/CeO₂(111), only the Γ-point sampling was used for k-point sampling and a Gaussian smearing width of 0.05 eV. Indeed, employing Γ-point sampling with the VASP_GAM version significantly reduces the computational cost, while the difference in total energy remains below 0.01 eV (Table S1).

Fig. 1 Top view (above) and side view (below) of the optimized structure for (A) Cu(111), (B) CeO₂(111) and (C) Cu₁/CeO₂(111) (colour code: light blue: Cu; gold: Ce; red: O).

CH₄ vibrational excitation analysis

As explained in the introduction, we compare the efficacy of vibrationally exciting CH₄ for promoting dissociation by using two methods. The first is the Fridman–Macheret (F–M) α model, which assesses the influence of plasma-induced vibrational excitations. This model was originally developed for reactions in the gas phase.¹⁸ However, it has also been used in gas–surface reactions to predict the enhancements in the dissociation rates by vibrational excitation, due to the lack of more accurate, feasible approaches available at that time.^14–17 In the F–M α model, the efficacy is determined by the proportionality constant α_FM:^14–17

Here, E^f_b and E^r_b are the energy barriers for the forward and reverse reactions, respectively. If E^f_b is equal to zero, there is no enthalpy barrier, and thus, the reaction is diffusion-limited. In this case, α_FM is equal to zero (i.e., the vibrationally excited levels play no role in enhancing the reaction). Vice versa, if E^r_b is equal to zero, the reaction is enthalpy-limited, and thus, the efficacy of the vibrationally excited levels to lower the reaction energy barrier is at maximum. In this case, α_FM is equal to 1.

The other method is the alternative η model proposed by Gerrits and Bogaerts, which fits vibrational efficacies obtained from MD simulations by incorporating three variables,¹⁹i.e., (1) the forward barrier height (E^f_b), (2) the ratio of the dissociating bond length between the transition state (TS) and the reactant (R_TS/R_gas), and (3) the SVP values, using the semi-empirical parameters of α₁ = 0.008259 mol kJ⁻¹, α₂ = 2.4405, and α₃ = 0.2032.

Specifically, the SVP values can be calculated by comparing the projection of a reactant normal mode onto the reaction coordinate at the transition state (TS). A larger SVP value indicates a stronger coupling between the reaction coordinate and the vibrational mode, thus resulting in more energy being available for the reaction.³⁹

In this study, we will evaluate these two approaches to compare the impact of vibrational excitation of CH₄ by calculating the forward rate constant k_v:

where A is the pre-exponential factor, E_b is the energy barrier of the dissociative chemisorption of CH₄, T is the gas temperature, k_B is the Boltzmann constant, E_v is the vibrational energy, and H(x) is the Heaviside step function (H(x) = 1 when x ≥ 0, and H(x) = 0 when x < 0). Specifically, the vibrational efficacy α can be calculated from eqn (2), or alternatively employ η calculated with eqn (3).

Results and discussion

CH₄ and CH₃ activation on the different catalyst surfaces

The three catalyst surfaces, i.e., Cu(111), CeO₂(111), and Cu₁/CeO₂(111), were investigated to compare their ability for C–H activation. The optimized CH₄ adsorption sites on these three surfaces are shown in Fig. S1–S3. We compare the CH₄ orientation based on one of the H atoms towards the surfaces. On Cu(111) and CeO₂(111), the most stable configuration corresponds to CH₄ adsorbed at a hollow site, where three H atoms lie approximately parallel to the surface plane. The resulting adsorption energies are −0.26 eV and −0.22 eV, respectively. In contrast, on Cu₁/CeO₂(111), CH₄ preferentially adsorbs at a bridge site because a single Cu atom already occupies the hollow site of CeO₂(111), leading to a slightly increased adsorption energy of −0.28 eV. These results (Table 1) indicate that CH₄ is only weakly physisorbed on all three surfaces. Although the orientation of CH₄ varies noticeably among the different adsorption sites, the adsorption energy differences remain small. The dissociative chemisorption (Fig. 2A, upper panel) of CH₄, i.e., CH₄ → CH₃* + H*, is endothermic on Cu(111) and CeO₂(111), with a reaction energy of 0.68 eV and 1.41 eV, respectively. The energy barrier of CH₄ on Cu(111) is 1.41 eV, which is similar to the result obtained using the optB86b-vdW functional.^29,68 It should be mentioned that the barrier height is affected by the choice of density of functional (DF).^29,69 For example, SRP32-vdW-DF1,^28,31,70 which is a DF generally performing well for the dissociative chemisorption of CH₄,^71,72 yields a considerably higher energy barrier of 1.72 eV for CH₄ on Cu(111).²⁹ The variations in barrier heights obtained from DFs will be further discussed below in relation to the mode-specificity analysis.

Table 1 Barrier height and reaction energy of CH₄ and CH₃ (dissociation) on Cu(111), CeO₂(111), and Cu₁/CeO₂(111)

Reaction	Surface	Barrier height (eV)	Reaction energy (eV)
CH4 → CH4*	Cu(111)	—	−0.26
CH4 → CH4*	CeO2(111)	—	−0.22
CH4 → CH4*	Cu1/CeO2(111)	—	−0.28
CH4* → CH3* + H*	Cu(111)	1.41	0.68
CH4* → CH3* + H*	CeO2(111)	1.56	1.49
CH4* → CH3* + H*	Cu1/CeO2(111)	0.35	−0.49
CH3 → CH3*	Cu(111)	—	−1.93
CH3 → CH3*	CeO2(111)	—	−2.21
CH3 → CH3*	Cu1/CeO2(111)	—	−2.16
CH3* → CH2* + H*	Cu(111)	1.43	0.88
CH3* → CH2* + H*	CeO2(111)	1.35	0.84
CH3* → CH2* + H*	Cu1/CeO2(111)	1.26	0.85

---

Fig. 2 (A) Energy diagram for CH₄ and CH₃ dissociation on Cu(111), CeO₂(111), and Cu₁/CeO₂(111). (B)–(D) The TS structures of CH₄ dissociation on Cu(111), CeO₂(111), and Cu₁/CeO₂(111). (E)–(G) The TS structures of CH₃ dissociation on Cu(111), CeO₂(111), and Cu₁/CeO₂(111) (colour code: light blue: Cu; gold: Ce; red: O; grey: C; white: H).

There is no evident difference in the barrier heights for C–H activation between the Cu(111) and CeO₂(111) surfaces (i.e., 1.41 vs. 1.56 eV, respectively), although the geometries of the TSs are quite different. On Cu(111), the CH₄ molecule can undergo dissociation to form Cu–CH₃ as an intermediate, enabling direct cleavage of the C–H bond alongside the coordination of the CH₃ group (Fig. 2B). However, this is not the case on CeO₂(111). Indeed, we find that the vibrational direction of the imaginary frequency (Fig. S4) does not coincide with that of C–H bond dissociation, referred to as a “pseudo-transition state”,⁷³ when we calculate the frequency of the TS structure of the O–CH₃ intermediate formed. Potentially, this can lead to a very different mode-specificity for the different surfaces. Notably, the other distinct TS structure found on CeO₂(111) corresponds to H abstraction from CH₄ to form ˙CH₃ radicals by electrophilic oxygen atoms from the catalyst surface (Fig. 2C).

The emergence of two distinct transition state (TS) geometries on the Cu(111) and CeO₂(111) surfaces stems from different mechanisms: the first TS features a stabilized CH₃ group by interacting with surface atoms, while a radical-like TS is formed instead when the geometric access or energetic favorability of CH₃–surface interaction is hindered.⁴ Notably, the Cu₁/CeO₂(111) surface reduces the barrier height of CH₄ dissociation to 0.35 eV, and the reaction becomes exothermic (reaction energy of −0.49 eV). Indeed, such a low barrier can also be obtained on the SAC Pd₁/CeO₂(111) with a similar exothermic process for dissociative chemisorption of CH₄.⁷⁴ Low-barrier catalysts for CH₄ at a low temperature may avoid unwanted side reactions to selectively convert CH₄ to value-added products.⁷⁵ Additionally, the Cu–CH₃ intermediate (Fig. 2D) can be formed in the TS for Cu₁/CeO₂(111), instead of the ˙CH₃ radical being suspended in the gas phase, as observed for CeO₂(111).

Next to dissociation on the surface, CH₄ can also be already dissociated in the plasma before it interacts with the catalyst, i.e., CH₄(g) → CH₃(g) + H(g), e.g., upon electron impact or upon reaction with other molecules or radicals. The CH₃ radicals produced in the gas phase can then adsorb on the catalyst surface. Therefore, we also compare the C–H bond dissociation of CH₃ on the three surfaces. The adsorption energy of CH₃ (Table 1) follows the order CeO₂(111) (−2.21 eV) > Cu₁/CeO₂(111) (−2.16 eV) > Cu(111) (−1.93 eV), indicating that a metal oxide surface (e.g., CeO₂) adsorbs CH₃ more strongly. The structures of adsorbed CH₃ for the three surfaces are shown in Fig. S5. Note that the adsorption energy of CH₃ is an order of magnitude larger than that of CH₄. However, the type of surface has only a limited effect on the barrier height for the dissociation of CH₃ (Fig. 2A, lower panel), with a maximum barrier height of 1.43 eV on Cu(111), similar to the results obtained with the optB86b-vdW DF,⁶⁸ and only slightly lower on CeO₂(111) (1.35 eV) and on Cu₁/CeO₂(111) (1.26 eV). The TS structures of CH₃ dissociation on the three surfaces are shown in Fig. 2E and F.

Electronic structure analysis

The different catalysts influence the barrier for CH₄ dissociation, resulting in distinct TS structures and energies, while showing very similar behavior towards CH₃ activation. To elucidate the effect of bonding on the dissociation of CH₄, we investigate the projected density of states (pDOS) and the charge density differences of the CH₄ TS on Cu(111), CeO₂(111), and Cu₁/CeO₂(111). The pDOS of a gaseous CH₄ molecule (Fig. S6) shows two clear peaks in which four H_1s orbitals overlap with the C_2s and C_2p orbitals. The stability of CH₄ in its tetrahedral structure arises from the formation of four equivalent C–H bonds, through sp³ hybridization between the C and H atoms.⁷⁶ When CH₄ is adsorbed on Cu(111), the H_1s–C_2p peak is split and a small peak appears nearby. Three peaks of H_1s–C_2p are observed if CH₄ is adsorbed on Cu₁/CeO₂(111), whereas only one peak appears on CeO₂(111). This means that the C_2p orbital of CH₄ can be split into C_2px, C_2py, and C_2pz on Cu₁/CeO₂(111).

In Fig. 3A, the overlaps between the H_1s–C_2s orbitals appear for the TS of CH₄ on all three surfaces. No evident spin polarization is observed on either the Cu(111) or the Cu₁/CeO₂(111) surface. When CH₄ dissociates on CeO₂(111), there is a strong spin polarization because the surface oxygen atom can abstract the H atom, resulting in a radical-like TS structure (Fig. 2C). The pDOS of the H_1s and C_2p orbital overlap shows multiple peaks for the Cu(111) surface. The overlap is the largest on the Cu₁/CeO₂(111) surface, indicative of a strong interaction, and thus likely the origin of the comparatively low dissociation barrier for CH₄.

Fig. 3 (A) Projected density of states (pDOS) of the C and H atoms in the TS of CH₄ on Cu(111), CeO₂(111), and Cu₁/CeO₂(111); the (B) positive and (C) negative of electron density difference isosurfaces (Δρ = ρ(TS) − ρ(surface) − ρ(CH₄*)) for the TS of CH₄ on Cu(111), CeO₂(111), and Cu₁/CeO₂(111), with the calculated values of Bader charges (ΔQ) of C and H atoms in between.

In Fig. 3B, the electron density difference isosurfaces show that the density increases around the C atom in the order Cu₁/CeO₂(111) > Cu(111) > CeO₂(111). The highest calculated Bader charge (i.e., 0.48 |e|, Table 2) is found for the Cu₁/CeO₂(111) surface, indicating the largest electron transfer. On Cu(111) (Fig. 3C), there is a slight electron density increase near the H atom; in this case, the Cu atom can act as an electron donor. However, the surface O atom from CeO₂(111) attracts the H atom in CH₄, which decreases the electron density near the H atom in the CH₄ TS. Similarly, the electron density near H is also reduced on Cu₁/CeO₂(111). Overall, the pDOS analysis and electron density differences show that the CH₄ TS interacts more strongly with Cu₁/CeO₂(111) than with the other surfaces, due to an increased electron transfer to and from the catalyst surface.

Table 2 Structural parameters, atomic Bader charge, and iCOHP of the TS for the C–H bond of CH₄ dissociation on the three surfaces

Surfaces	R(C–H)/Å		Q/e		iCOHP
	Initial state	Transition state	Q C	Q H
Cu(111)	1.10	1.72	0.27	0.12	−1.38
CeO2(111)	1.10	1.58	0.25	−0.48	−1.81
Cu1/CeO2(111)	1.10	1.39	0.48	−0.38	−3.55

---

Furthermore, the crystal orbital Hamiltonian population (COHP) is used to analyze the bonding interaction between atoms. The calculated projected COHPs (pCOHPs) were compared to consider the C–H bond interaction strength in the initial states (Fig. S7) and TSs (Fig. 4) of CH₄ on Cu(111), CeO₂(111), and Cu₁/CeO₂(111). The integral of the COHP up to the Fermi level (termed iCOHP) is used to qualitatively measure the strength of a chemical bond between two atoms. The more negative the iCOHP is, the stronger the bond is. As shown in Table 2, the largest value of iCOHP is found for the CH₄ TS on the Cu₁/CeO₂(111) surface (with a value of −3.55), suggesting a strong bond. Correspondingly, the C–H bond length is the shortest among the three surfaces. The barrier height is lower if an intermediate (i.e., TS) is more stable. Therefore, the pCOHPs and iCOHPs results demonstrate that the TS on the Cu₁/CeO₂(111) surface is the most stable structure among the three surfaces, indicating the lowest barrier height (cf.Table 1).

Fig. 4 The crystal orbital Hamilton population (COHP) for the transition states of CH₄ on (A) Cu(111), (B) CeO₂(111), and (C) Cu₁/CeO₂(111).

Mode specificity analysis by the SVP model

As mentioned in the introduction, vibrational excitation of CH₄ can promote its reactivity and improve product selectivity compared to thermal conversion.^15,46,77 The CH₄ molecule has one C atom tetrahedrally bonded to four H atoms. Therefore, the total number of degrees of freedom is 15, i.e., 9 vibrational modes, 3 translation modes, and 3 rotational modes. Since some of the vibrational modes are degenerate, CH₄ has 4 fundamental vibrational modes: the bending modes (ν2 and ν4), a symmetric stretching mode (ν1), and an asymmetric stretching mode (ν3). Here, the calculated vibrational frequencies of gaseous CH₄ (Table S2) are in good agreement with experiments.⁴⁷

To investigate how these modes might affect CH₄ dissociation on the three surfaces studied here, we use the SVP model to compare the gas phase vibrational modes with the reaction coordinate vector at the saddle point. In this study, CH₄ in the gas phase is taken as the initial structure for the SVP calculation to mimic the process of the excited CH₄ molecule reaching the surface under plasma conditions, rather than CH₄ being physisorbed on the surface, which is directly followed by dissociative chemisorption. Note that activated dissociative chemisorption in general should be modelled as a direct reaction of the gaseous reactant with the surface without thermalization, instead of prior physisorption and thermalization of the adsorbate with the surface. The calculated frequencies for the SVP calculations are listed in Table S3. As shown in Fig. 5A, the SVP analysis shows that the stretching modes (v1, v3) of CH₄ have the dominant contribution in the C–H bond dissociation of CH₄ on the Cu(111) surface. The SVP value of the symmetric stretching mode (v1) on Cu(111) is 0.36, consistent with the value for CH₄ + Ni(111).³⁵ There are three degenerate asymmetric stretching modes (v3) with the highest SVP value of 0.72. However, the v3 mode is three-fold degenerate, thus the average SVP value of v3 over all three vectors is slightly lower than that of v1. This means that symmetric stretching has a larger contribution, in accordance with the experimental and theoretical results on Cu(111) and other metal surfaces.^{28,29,35,42,43}

Fig. 5 SVP values of the vibrational modes of gaseous CH₄ and CH₃ (A)–(F) and adsorbed CH3 (G)–(I) onto the reaction coordinate at the TS on the surfaces of (A), (D) and (G) Cu(111), (B), (E) and (H) CeO₂(111), and (C), (F) and (I) Cu₁/CeO₂(111), with the comparison of α values (horizontal purple dotted line) calculated by the F–M α model (v1: symmetric stretching; v3: asymmetric stretching; v2: twisting; v4: scissoring).

The SVP model (Fig. 5B–C) predicts similar results (v1 > v3) for CeO₂(111) and Cu₁/CeO₂(111). Furthermore, the SVP values for the bending modes (i.e., v2 and v4) on Cu₁/CeO₂(111) and Cu(111) are similar. Interestingly, on CeO₂(111), the v2 mode yields an SVP value of almost zero. This is caused by the distinct geometry structure of the TS on CeO₂(111), in which the nearly perpendicular orientation of the TS structure (Fig. 2C) yields a limited projection of the bending mode onto the reaction coordinate, leading to the SVP value being close to zero.

Moreover, we compare the dependence of the SVP calculations for the choice of DFs on CH₄ dissociation over three surfaces in Table S4. As expected, the choice of DF significantly affects the barrier height, as well as the reaction energy, since it is proportional to the barrier height.²⁹ In particular, on Cu(111), the barrier height increases from 1.41 eV (PBE) to 1.92 eV (BEEF-vdW), indicating strong functional dependence. For Cu(111) and Cu₁/CeO₂(111) surfaces, however, the SVP values remain nearly constant, as the TS geometries and imaginary frequencies remain similar across different functionals, suggesting that the SVP and the shape of the PES are relatively insensitive to the choice of DF.⁷⁸ However, the larger functional sensitivity observed for CeO₂(111) can be attributed to distinct TS structures with ˙CH₃ radical being suspended in the gas phase (Fig. 2C). For example, we employed DFT+U (U = 4.5 eV) to compare the SVP values on CeO₂(111), in which PBE yields an SVP of 0.35, whereas SRP32-vdW-DF1 and BEEF-vdW give higher values of 0.43 and 0.46, respectively, indicating a modestly larger vibrational contribution by changing different treatments of DFs and dispersion interactions. Overall, comparisons among various DFs show that the variation in SVP values remains relatively small.

It should be mentioned that generally the SVP values of translation and rotation motion are less useful, possibly due to the complex nature of dynamical effects arising from these degrees of freedom. For example, the large rotational efficacy observed for HCl + Au(111) is due to the shape of the PES and concomitant dynamics, and, therefore, cannot be captured by the PES surrounding the TS.⁷⁹ Likewise, the bobsled effect in late barrier systems is associated with translational motion and is caused by the curvature of the reaction path, again prior to reaching the TS.³³ Thus, the static SVP calculation is more reliable for the vibrational modes. Additionally, the vibrational temperature in plasma can be much higher than the gas temperature, indicating a high prevalence of vibrationally excited reactants.^80,81 Therefore, we mainly discuss the effects of vibrational modes by SVP in this paper.

Furthermore, the dissociation of the CH₃ radical on the three surfaces is also investigated using the SVP model. Two options are considered: first, the CH₃ radical might dissociate directly from the gas phase, without sufficient time to reorient the planar geometry to the more stable bent adsorption geometry. The chemisorption of CH₃ (Fig. S8A) is highly exothermic, with an adsorption energy of 1.93 eV on Cu(111). The adsorption energy is 0.50 eV larger than the CH₃ dissociation barrier (1.43 eV, Table 1). Similarly, the adsorption energies of CH₃ on CeO₂(111) and Cu₁/CeO₂(111) are also larger than the CH₃ dissociation barrier (Fig. S8B and C). These results indicate that CH₃ dissociation might occur directly during or after CH₃ binds to the surface. Thus, we use the SVP model to understand the mode specificity of CH₃ dissociation on the three surfaces. As shown in Fig. 5D and E, the bending modes (v2 and v4) of gaseous CH₃ have a larger contribution to CH₃ dissociation than the stretching modes (v1 and v3) on both Cu(111) and CeO₂(111) because the bending modes align the planar gaseous CH₃ radical with the bent configuration at the TS. In contrast, the stretching modes (Fig. 5F) still take a leading role in CH₃ dissociation on Cu₁/CeO₂(111) because there is a large projection from the stretching vibration onto the reaction coordinate vector at the saddle point. The SVP values of CH₃ dissociation are listed in Table S5.

Secondly, CH₃ might adsorb without immediate subsequent dissociation, giving CH₃ sufficient time to reorient from planar to bent. Furthermore, in that case, the large chemisorption energy needs to be dissipated. Several important dissipation channels might exist, such as vibrational excitation and phonons. Vibrational excitation could again lead to increased reactivity. Moreover, other plasma–catalytic processes might be able to vibrationally excite adsorbates like CH₃. Energy transfer to the phonons, on the other hand, could reduce the reactivity, since there would be less energy available for the reaction. However, it is also likely that this process is relatively slow. In Fig. 5G–I, the stretching modes of adsorbed CH₃ yield larger SVP values than the bending modes on all three surfaces, which is the opposite of when the gas phase vibrational modes are employed. In other words, there are large differences in vibrational efficacy (i.e., stretching vs. bending) for CH₃ dissociation on Cu(111) and CeO₂(111), depending on whether the TS is reached directly from the gaseous or chemisorbed state. The relevant state and concomitant reaction mechanism might be determined in future work using vibrational state-specific MD calculations and molecular beam experiments to investigate the radical chemistry in plasma catalysis.

Vibrational efficacy calculations

Vibrational efficacies are calculated using the F–M α model and the η model for the three surfaces. As shown in Table 3 and Fig. 5A–C, higher α values can be obtained on the Cu(111) and CeO₂(111) surfaces when CH₄ dissociation is endothermic for the F–M α model, while a lower α value (0.29) can be observed on Cu₁/CeO₂(111) when CH₄ dissociation is exothermic. As for the η model, the efficacy of CH₄ excitation is estimated by focusing on the dominant contribution from various vibrational modes, where symmetric stretching (v1) exhibited the highest SVP value on all three surfaces. Besides, the other two features, i.e., the forward barrier height (E^f_b) and the ratio of the dissociating bond length between the TS and the reactant (R_TS/R_gas), have been incorporated to understand the vibrational state-specific efficacies.¹⁹ Evidently, the calculated vibrational efficacies η are much larger than the values predicted by the F–M α model on the Cu(111) and CeO₂(111) surfaces. However, the difference between η and α are not so significant on Cu₁/CeO₂(111), since the forward barrier height is comparatively low. Generally, the η model deviates more strongly from the F–M α model when the absolute barrier height is large. Nevertheless, vibrational mode specific reactivity has been observed for low barrier heights as well, highlighting the importance of taking into account the coupling of specific vibrational modes to the reaction coordinate.¹⁹

Table 3 Vibrational efficacies of CH₄ dissociation by using the F–M α and η models on the three investigated surfaces

Surfaces	Reaction type	α F–M	E ^fb (eV)	R TS/Rgas	SVP (= v1)	η	k η/kFM (500 K)
Cu(111)	Endothermic	0.66	1.41	1.56	0.36	1.67	5846
CeO2(111)	Endothermic	0.96	1.56	1.44	0.35	1.69	243
Cu1/CeO2(111)	Exothermic	0.29	0.35	1.26	0.37	0.34	1.54

---

The integration of vibrational efficacy into existing kinetic models enables the prediction of gas–surface reaction rates dependent on the specific vibrational distribution and the evaluation of how vibrational non-equilibrium in a plasma affects the overall process at a macroscopic level.^14–17 Here, we compare the ratio (k_η/k_FM) of forward rate constants between the η and F–M α model for the three surfaces at different temperatures (Fig. S9). The ratio disparities are greater at low temperatures and diminish as the temperature increases, indicating that the vibrational efficacy can significantly change rate coefficients at low temperatures. For example, the η model yields a rate coefficient that is 5846 times higher on Cu(111) and 243 times higher on CeO₂(111) compared to those obtained from the F–M α model. In contrast, the difference between α and η is quite small on Cu₁/CeO₂(111), thus hardly affecting the reaction rate.

The main reason for the discrepancy between the F–M α model and the SVP model lies in the fact that CH₄ dissociates differently on the three surfaces with the distinct structures of the TS. As for the F–M α model, the vibrational efficacy is solely computed based on the ratio between the forward and backward reaction barrier heights. This approach presents some critical limitations, including its oversimplified dependence on barrier height ratios without rigorous TS validation, its neglect of mode-specific vibrational effects, and failure to account for coupling between vibrational energy and various DOFs.¹⁹ Additionally, its arbitrary restriction (α ∈ [0, 1]) contradicts experimental observations where vibrational energy surpasses translational efficacy for dissociative chemisorption of CH₄ on catalyst surfaces.¹⁹ Therefore, the prediction of vibrational efficacy by the F–M α model is less convincing for the gas–surface reactions. Vibrational efficacy predictions using the η model can be taken at comparable computational cost to understand the importance of mode specificity and the structural dependence of the catalyst. Notably, at low temperatures, the two models diverge significantly in their predictions of the reaction rate, implying that the role of vibrational excitation may be underestimated using the F–M α model, which merits further investigation through combined experimental and theoretical molecular beam studies. These observations were also made in an extensive investigation of the dissociative chemisorption of N₂ on Ru(0001), by comparing various transition state theory-based models with MD simulations.²⁰

Finally, we emphasize the importance of further validating the SVP results for plasma–catalytic CH₄ activation. Performing dynamical calculations of CH₄ dissociation on a realistic PES can offer more detailed insights into reactive collisions and capture the energy dependence of state-specific reaction probabilities.^28,29 These insights might help to further develop and validate the SVP model and the fitted η approach to account for the effects of vibrational non-equilibrium in plasma catalysis.

Conclusions

We studied the dissociation of both CH₄ and CH₃ on Cu(111), CeO₂(111), and Cu₁/CeO₂(111) with DFT calculations. We find that the TS structures for CH₄ are different for the three surfaces. Compared to Cu(111) and CeO₂(111), the single Cu atom supported on CeO₂(111) can significantly reduce the barrier height (0.35 eV), which is attributed to an increased molecule–surface electron transfer and a more stable TS structure, indicated by electronic structure analysis. However, these effects are considerably less important for dissociating the chemisorbed CH₃, where the barrier height is similar on all three investigated surfaces.

The effect of vibrational excitation of CH₄ and CH₃ was investigated with the Fridman–Macheret α model and a novel, alternative approach that is fitted to vibrational efficacies and DFT results in the literature. Notably, the prediction of vibrational efficacy by the F–M α model seems to be less reliable due to its lack of mode specificity or structure dependence, potentially leading to an underestimation of the role of vibrational excitation at low temperatures. Furthermore, the SVP results indicate that the stretching modes of CH₄ play a primary role in its dissociation on these three surfaces, in qualitative agreement with previous experimental and theoretical studies. The relative vibrational efficacies for CH₃ dissociation of the stretching and bending modes show large differences, depending on the specific reaction dynamics. Future MD and molecular beam studies focusing on the vibrational efficacies can help elucidate the reaction mechanism, namely, whether the CH₃ radical reacts directly from the gas phase or first reorients upon chemisorption. Moreover, the vibrational efficacies are compared between the F–M α model and the η model, which show differences of up to three orders of magnitude in computed reaction rates, in particular at low temperature. Future vibrational state-specific molecular beam studies should be able to validate our predictions regarding the dissociative chemisorption of vibrationally excited CH₄ on catalyst surfaces.

Overall, these results offer valuable insights into catalytic C–H activation and the impact of vibrational excitation and may be of help in developing efficient catalysts for plasma–catalytic CH₄ conversion. We hope that our study will inspire further exploration through high-dimensional (quantum) dynamical calculations and experiments, enriching our comprehension of reaction dynamics in plasma catalysis.

Author contributions

The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript. Shangkun Li: conceptualization, validation, formal analysis, resources, data curation, writing original draft, and writing – review and editing. Santu Luo: conceptualization, validation, formal analysis, resources, data curation, writing original draft, and writing – review and editing. Rui Liu: conceptualization, validation, and writing – review and editing. Zhaolun Cui: conceptualization, validation, and writing – review and editing. Yanhui Yi: conceptualization, validation, and writing – review and editing. Erik C. Neyts: formal analysis, resources, data curation, writing – original draft, writing – review and editing, and supervision. Annemie Bogaerts: formal analysis, resources, data curation, writing – original draft, writing – review and editing, supervision, and funding acquisition. Nick Gerrits: formal analysis, resources, data curation, writing – original draft, and writing – review and editing.

Conflicts of interest

There are no conflicts to declare.

Data availability

The authors confirm that the data supporting the results of this research article are available in the article and its supplementary information (SI). Supplementary information contains convergence and CPU timing benchmarks; depiction and several computational results of various adsorbed and transition states; vibrational mode analysis; dependence of barriers on the density functional; temperature dependence of k_η/k_FM. See DOI: https://doi.org/10.1039/d5cp03569d.

Additional data can be made available, upon reasonable request.

Acknowledgements

The research was supported by the China Scholarship Council and the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant Agreement No. 810182 – SCOPE ERC Synergy project). We acknowledge financial support from the National Natural Science Foundation of China [22472018 and 22272015]. The computational resources and services used in this work were provided by the HPC core facility CalcUA of the Universiteit Antwerpen, and VSC (Flemish Supercomputer Center). We thank Dr Jintao Sun and Roel Michiels for their many valuable discussions.

Notes and references

1. P. Schwach, X. Pan and X. Bao, Chem. Rev., 2017, 117(13), 8497–8520.
2. S. Li, R. Ahmed, Y. Yi and A. Bogaerts, Catalysts, 2021, 11, 590.
3. N. Gerrits, B. Jackson and A. Bogaerts, J. Phys. Chem. Lett., 2024, 15(9), 2566–2572.
4. A. A. Latimer, A. R. Kulkarni, H. Aljama, J. H. Montoya, J. S. Yoo, C. Tsai, F. Abild-Pedersen, F. Studt and J. K. Norskov, Nat. Mater., 2017, 16(2), 225–229.
5. L. Lefferts, Angew. Chem., Int. Ed., 2024, 63(10), e202305322.
6. A. Bogaerts, Nat. Chem. Eng., 2025, 2, 336–340.
7. A. Bogaerts, X. Tu, J. C. Whitehead, G. Centi, L. Lefferts, O. Guaitella, F. Azzolina-Jury, H.-H. Kim, A. B. Murphy, W. F. Schneider, T. Nozaki, J. C. Hicks, A. Rousseau, F. Thevenet, A. Khacef and M. Carreon, J. Phys. D: Appl. Phys., 2020, 53(44), 443001.
8. J. Van Turnhout, K. Rouwenhorst, L. Lefferts and A. Bogaerts, EES Catal., 2025, 3, 669–693.
9. Y. Yi, S. Li, Z. Cui, Y. Hao, Y. Zhang, L. Wang, P. Liu, X. Tu, X. Xu, H. Guo and A. Bogaerts, Appl. Catal., B, 2021, 296, 120384.
10. E. C. Neyts, K. K. Ostrikov, M. K. Sunkara and A. Bogaerts, Chem. Rev., 2015, 115(24), 13408–13446.
11. B. Loenders, R. Michiels and A. Bogaerts, J. Energy Chem., 2023, 85, 501–533.
12. T. Nozaki, N. Muto, S. Kadio and K. Okazaki, Catal. Today, 2004, 89(1–2), 67–74.
13. D. Y. Kim, H. Ham, X. Chen, S. Liu, H. Xu, B. Lu, S. Furukawa, H. H. Kim, S. Takakusagi, K. Sasaki and T. Nozaki, J. Am. Chem. Soc., 2022, 144(31), 14140–14149.
14. P. Mehta, P. Barboun, F. A. Herrera, J. Kim, P. Rumbach, D. B. Go, J. C. Hicks and W. F. Schneider, Nat. Catal., 2018, 1(4), 269–275.
15. Y. Engelmann, K. van’t Veer, Y. Gorbanev, E. C. Neyts, W. F. Schneider and A. Bogaerts, ACS Sustainable Chem. Eng., 2021, 9(39), 13151–13163.
16. Y. Engelmann, P. Mehta, E. C. Neyts, W. F. Schneider and A. Bogaerts, ACS Sustainable Chem. Eng., 2020, 8(15), 6043–6054.
17. R. Michiels, Y. Engelmann and A. Bogaerts, J. Phys. Chem. C, 2020, 124(47), 25859–25872.
18. A. Fridman, Plasma Chemistry, Cambridge Univ. Press, New York, NY, 2008.
19. N. Gerrits and A. Bogaerts, EES Catal., 2025, 3, 733–742.
20. F. van den Bosch, N. Gerrits and J. Meyer, EES Catal., 2025, 3(6), 1257–1271.
21. K. M. Bal and E. C. Neyts, J. Phys. D: Appl. Phys., 2021, 54(39), 394004.
22. S. Roy, N. K. J., N. Tiwari and A. K. Tiwari, Int. Rev. Phys. Chem., 2020, 39(3), 267–318.
23. M. Karplus, R. N. Porter and R. D. Sharma, J. Chem. Phys., 1965, 43(9), 3259–3287.
24. K. M. Bal, A. Bogaerts and E. C. Neyts, J. Phys. Chem. Lett., 2020, 11(2), 401–406.
25. W. H. Miller, N. C. Handy and J. E. Adams, J. Chem. Phys., 1980, 72, 99–112.
26. B. Jackson and S. Nave, J. Chem. Phys., 2011, 135, 114701.
27. S. Roy and A. K. Tiwari, Phys. Chem. Chem. Phys., 2022, 24, 16596–16610.
28. N. Gerrits, H. Chadwick and G. J. Kroes, J. Phys. Chem. C, 2019, 123(39), 24013–24023.
29. N. Gerrits, D. Migliorini and G. J. Kroes, J. Chem. Phys., 2018, 149(22), 224701.
30. N. Gerrits, J. Geweke, E. W. F. Smeets, J. Voss, A. M. Wodtke and G. J. Kroes, J. Phys. Chem. C, 2020, 124(29), 15944–15960.
31. F. Nattino, C. Diaz, B. Jackson and G. J. Kroes, Phys. Rev. Lett., 2012, 108(23), 236104.
32. B. Jackson, F. Nattino and G. J. Kroes, J. Chem. Phys., 2014, 141(5), 054102.
33. X. Zhou, B. Jiang and H. Guo, J. Phys. Chem. C, 2019, 123(34), 20893–20902.
34. R. Moiraghi, A. Lozano, E. Peterson, A. Utz, W. Dong and H. F. Busnengo, J. Phys. Chem. Lett., 2020, 11(6), 2211–2218.
35. B. Jiang, R. Liu, J. Li, D. Xie, M. Yang and H. Guo, Chem. Sci., 2013, 4(8), 3249–3254.
36. B. Jackson and S. Nave, J. Chem. Phys., 2011, 135(11), 114701.
37. P. M. Hundt, M. E. van Reijzen, H. Ueta and R. D. Beck, J. Phys. Chem. Lett., 2014, 5(11), 1963–1967.
38. B. Jiang and H. Guo, J. Chem. Phys., 2013, 138(23), 234104.
39. H. Guo and B. Jiang, Acc. Chem. Res., 2014, 47(12), 3679–3685.
40. F. Nattino, H. Ueta, H. Chadwick, M. E. van Reijzen, R. D. Beck, B. Jackson, M. C. van Hemert and G. J. Kroes, J. Phys. Chem. Lett., 2014, 5(8), 1294–1299.
41. H. Guo, A. Farjamnia and B. Jackson, J. Phys. Chem. Lett., 2016, 7(22), 4576–4584.
42. P. Maroni, D. C. Papageorgopoulos, M. Sacchi, T. T. Dang, R. D. Beck and T. R. Rizzo, Phys. Rev. Lett., 2005, 94(24), 246104.
43. S. Nave and B. Jackson, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81(23), 233408.
44. A. Wang, J. Li and T. Zhang, Nat. Rev. Chem., 2018, 2(6), 65–81.
45. J.-C. Liu, Y. Tang, Y.-G. Wang, T. Zhang and J. Li, Natl. Sci. Rev., 2018, 5, 5.
46. J. Sun, Q. Chen, Y. Guo, Z. Zhou and Y. Song, J. Energy Chem., 2020, 46, 133–143.
47. T. Butterworth, A. van de Steeg, D. van den Bekerom, T. Minea, T. Righart, Q. Ong and G. van Rooij, Plasma Sources Sci. Technol., 2020, 29(9), 095007.
48. S. Li, J. Sun, Y. Gorbanev, K. van’t Veer, B. Loenders, Y. Yi, T. Kenis, Q. Chen and A. Bogaerts, ACS Sustainable Chem. Eng., 2023, 11(42), 15373–15384.
49. G. F. Kresse and J. Hafner, Comput. Mater. Sci., 1996, 6, 15–50.
50. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169.
51. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47(1), 558–561.
52. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
53. P. E. Blochl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50(24), 17953–17979.
54. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
55. S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys and A. P. Sutton, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 57, 1505–1509.
56. Z.-Q. Huang, T. Zhang, C.-R. Chang and J. Li, ACS Catal., 2019, 9(6), 5523–5536.
57. L. Chen, I. A. W. Filot and E. J. M. Hensen, ACS Catal., 2023, 13(23), 15230–15247.
58. Z.-Q. Huang, L.-P. Liu, S. Qi, S. Zhang, Y. Qu and C.-R. Chang, ACS Catal., 2017, 8(1), 546–554.
59. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu, J. Chem. Phys., 2010, 132(15), 154104.
60. G. Henkelman, B. P. Uberuaga and H. Jónsson, J. Chem. Phys., 2000, 113(22), 9901–9904.
61. G. Henkelman and H. Jónsson, J. Chem. Phys., 1999, 111(15), 7010–7022.
62. V. L. Deringer, A. L. Tchougreeff and R. Dronskowski, J. Phys. Chem. A, 2011, 115(21), 5461–5466.
63. S. Maintz, V. L. Deringer, A. L. Tchougreeff and R. Dronskowski, J. Comput. Chem., 2016, 37(11), 1030–1035.
64. W. Tang, E. Sanville and G. Henkelman, J. Phys.: Condens. Matter, 2009, 21(8), 084204.
65. Y. Zhang, S. Zhao, J. Feng, S. Song, W. Shi, D. Wang and H. Zhang, Chem, 2021, 7, 2022–2059.
66. J. L. Chen, X. C. Jiang, L. Feng, J. Zhu, J. W. Zhao, J. X. Liu and W. X. Li, Nat. Commun., 2025, 16, 9144.
67. W. Ji, N. Wang, X. Chen, Q. Li, K. Lin, J. Deng, J. Chen and X. Xing, Inorg. Chem., 2022, 61(26), 10006–10014.
68. M. D. Marcinkowski, M. T. Darby, J. Liu, J. M. Wimble, F. R. Lucci, S. Lee, A. Michaelides, M. Flytzani-Stephanopoulos, M. Stamatakis and E. C. H. Sykes, Nat. Chem., 2018, 10(3), 325–332.
69. T. Tchakoua, N. Gerrits, E. W. F. Smeets and G. J. Kroes, J. Chem. Theory Comput., 2023, 19(1), 245–270.
70. H. Guo and B. Jackson, J. Chem. Phys., 2019, 150(20), 204703.
71. N. Gerrits, Front. Chem., 2024, 12, 1481235.
72. H. Chadwick, D. Migliorini and G. J. Kroes, J. Chem. Phys., 2018, 149(4), 044701.
73. M. D. Krcha, A. D. Mayernick and M. J. Janik, J. Catal., 2012, 293, 103–115.
74. S. Tomar, B. S. Bhadoria, H. Jeong, J. H. Choi, S.-C. Lee and S. Bhattacharjee, J. Phys. Chem. C, 2024, 128, 8580–8589.
75. Z. Liang, T. Li, M. Kim, A. Asthagiri and J. F. Weaver, Science, 2017, 356, 299–303.
76. S. S. Sason, C. H. Phillippe, W. John and H. Sons, Chrom., 2009, 69(11–12), 1493.
77. P. Mehta, P. Barboun, D. B. Go, J. C. Hicks and W. F. Schneider, ACS Energy Lett., 2019, 4(5), 1115–1133.
78. A. D. Powell, N. Gerrits, T. Tchakoua, M. F. Somers, H. F. Busnengo, J. Meyer, G. Kroes and K. Doblhoff-Dier, J. Phys. Chem. Lett., 2024, 15, 307–315.
79. N. Gerrits, J. Geweke, D. J. Auerbach, R. D. Beck and G. J. Kroes, J. Phys. Chem. Lett., 2021, 12(30), 7252–7260.
80. J. C. Whitehead, J. Phys. D: Appl. Phys., 2016, 49(24), 243001.
81. J. J. Qiao, Q. Yang, D. Z. Wang and Q. Xiong, Plasma Sources Sci. Technol., 2023, 32(11), 11LT01.

---

Footnote

† S. Li and S. Luo contributed equally.

---