Enhanced C-H bond activation by tuning the local environment of surface lattice oxygen of MoO3.

The lattice oxygen on transition metal oxides serves as a critical active site in the dehydrogenation of alkanes, whose activity is determined by electronic properties and environmental structures. Hydrogen affinity has been used as a universal descriptor to predict C-H bond activation, while the understanding of the environmental structure is ambiguous due to its complexity. This paper describes a combined theoretical and experimental study to reveal the activity of lattice oxygen species with different local structures, taking Mo-based oxides and C-H bond activation of low-carbon alkanes as model catalytic systems. Our theoretical work suggests that oxygen species with convex curvature are more active than those with concave curvature. Theoretically, we propose an interpretative descriptor, the activation deformation energy, to quantify the surface reconstruction induced by adsorbates with various environmental structures. Experimentally, a Mo-based polyoxometalate with the convex curvature structure shows nearly five times the initial activity than single-crystal molybdenum oxide with the concave one. This work provides theoretical guidance for designing metal oxide catalysts with high activity. This journal is © The Royal Society of Chemistry.

Introduction

Developing highly effective alkane dehydrogenation catalysts is a critical technical challenge for producing olefins from shale gas, a rich source of light hydrocarbons (CH4, C2H6, and C3H8).[1] The oxidation of alkanes begins with C–H bond activation, which is generally considered the rate-determining step for the whole reaction.[2] Due to the inert C–H bond of low-carbon alkanes, effective and selective activation of the C–H bond of low-carbon alkanes has always been challenging in designing such catalysts.[1] Despite non-oxidative dehydrogenation being the industrial technology, the process is limited to using costly platinum or toxic chromium catalysts. Among the non-noble catalysts explored, some reducible transition metal oxides are promising ones.[3-5] The surface lattice oxygen is believed to be a key active site for alkane dehydrogenation. The activation of alkane molecules occurs through the Mars–van Krevelen mechanism, which entails oxidation of hydrocarbons with reduction of the oxide surface through losing surface lattice oxygen species.[6] Generally, the activity of lattice oxygen is determined by the electronic structure and the surrounding environment. However, because a single surface lattice oxygen can activate the C–H bond, many descriptors are proposed to predict the activity of oxygen species solely, such as hydrogen affinity (EH),[7] oxygen vacancy formation energy,[8] and O-2p band center.[9] In contrast, recent studies clearly showed the importance of the local environment near the active site, especially for metal–organic framework (MOF-based) materials[10] and zeolite catalysts.[11] The catalytic microenvironment could form a unique spatial structure accompanied by various electronic properties around the active catalytic sites, which plays an essential role in regulating the catalytic performance.[12,13] However, there is a lack of clear and general comprehension of how the local environment affects the activity. Therefore, an insightful understanding of how the local environment around the lattice oxygen is critical for designing transition metal oxide catalysts with effective C–H bond activation is essential.

Molybdenum trioxide (MoO3) is a reducible oxide used in crucial partial oxidation reactions, such as methanol production[14],[15] alkene metathesis,[16] and hydrocracking.[16] In particular, it has been used in many studies of C–H bond activation of alkanes and alcohols. For example, Choksi et al. studied the partial oxidation of methanol on MoO3 by DFT and microkinetic simulation and elaborated the universal single crystal model reaction for catalytic kinetic modeling.[17] Amakawa et al. used molybdenum oxides as a model system to study the strain effect on the activation of propane oxidation.[18] Yao et al. investigated the kinetics of oxidative and non-oxidative dehydrogenation of ethane on molybdenum oxide and developed a generalized mechanistic framework.[19] These studies have shown that molybdenum oxide is a common and suitable model system for studying C–H bond activation.

This paper describes the influence of the environmental structure of different oxygen species on the reactivity of propane dehydrogenation employing α-MoO3 as a model of transition metal oxide catalysts. α-MoO3 with an orthorhombic cell and pbmn space group is the most common crystal structure in molybdenum oxides. Previous studies have shown that the (010) crystal plane is the most stable and dominant.[20] The α-MoO3 (010) crystal plane provides three distinct oxygen atoms: terminal (Ot), asymmetrical bridging (Oa), and symmetrical bridging (Os) oxygen atoms (Fig. 1a), establishing a great platform with various environmental structures to study the effect. However, molybdenum oxide's low C–H bond activation ability restricts its industrial application. Our theoretical study shows that the bridging oxygen with higher C–H breaking ability is limited by the unfavorable steric structure, which is the main reason for its poor surface activity. Furthermore, through theoretical calculations and experiments, we demonstrated that optimizing such a spatial structure near the surface lattice oxygen can greatly improve its C–H bond activation ability by using a molybdenum-based Keggin-polyoxometalate (POM), one kind of cluster-type molybdenum oxide, as the model catalyst.

Fig. 1

Identification of different oxygen species activity and propane dehydrogenation mechanism. (a) Optimized structures of the α-MoO3 (010) pristine surface and three distinctive lattice oxygen species. (b) C–H bond activation transition states on three lattice oxygen species. (c) Potential energy diagram and transition state structures of the first step, C–H bond activation, over three oxygen species.

Results and discussion

Surface lattice oxygen activity on the α-MoO3 model

Firstly, it is imperative to investigate the activity differences between different oxygen species. Many previous microkinetic and experimental studies on propane dehydrogenation have demonstrated that the first step of C–H bond activation on the methylene of propane is the rate-determining step (RDS) of the overall dehydrogenation reaction.[21-24] Therefore, methylene C–H activation is considered the most crucial elementary reaction. We calculated propane dehydrogenation's activation energies (Ea) at three different oxygen sites (Fig. 1). The terminal oxygen (MoOt) was calculated to be the most active site (Ea = 1.15 eV), while asymmetric bridging oxygen (Mo–Oa–Mo) (Ea = 1.39 eV) was less active. The dehydrogenation of symmetric bridging oxygen (Mo–Os–Mo) was hindered by too high an energy barrier (Ea = 1.74 eV). Consistently, the oxygen vacancy formation energies of Ot and Oa (Table S2†) are similar, while the oxygen vacancy formation energy of Os is much higher. The above results indicate that Os is the most stable and inactive oxygen species. Therefore, Mo–Os–Mo will not be discussed in the following text unless stated. It is noteworthy that through differential charge density analysis, the metal–oxygen species participate in hydrogen abstraction as Lewis acid–base pairs (Fig. S4†), which means that the lattice oxygen (O2−) acts as a Lewis base site for accepting protons and eventually the metal cations as Lewis acid sites accept electrons to be reduced during the activation of the C–H bond.

For a simpler prediction of oxygen activity, Nørskov and co-workers proposed that hydrogen affinity (EH) is one universal descriptor of reactivity for radical-like hydrocarbon activation.[7] They found that EH is proportional to the C–H activation energy of alkanes in various materials. We also calculated EH at three oxygen sites and tried to correlate them with propane activation energy. We found that the EH sequence of the three oxygen species was: Mo–Oa–Mo (EH = −0.48 eV) < MoOt (EH = −0.19 eV) < Mo–Os–Mo (EH = −0.06 eV). However, when comparing the EH and C–H activation energies of Mo–Oa–Mo and MoOt, there was an anomaly (C–H activation energy: MoOt < Mo–Oa–Mo). To search for potential and common patterns, we made the same calculation and comparison for methane and ethane (the homologs of propane) at three active sites (Fig. 2a). Anomalies exist in all of them. Such an anomaly indicates that other factors besides the lattice oxygen greatly affect the activation of propane when compared with terminal oxygen and bridging oxygen. We speculated that this is due to the environmental effect around different oxygen species.

Fig. 2

Effects of different local environmental geometries on oxygen species activity. (a) Comparison of the first step C–H bond activation energy and EH of C1–C3 alkanes on different oxygen species. (b) Non-covalent interaction analysis for PDH at MoOt and Mo–Oa–Mo via the independent gradient model (IGM). In the transition state distance: unreconstructed surface and inactivated C3H8 (left); adsorbate-induced reconstructed surface and activated C3H8 (right). (c) Analysis of the first step, dehydrogenation of propane, in MoOt and Mo–Oa–Mo via the crossing potential model.

Environmental effects on the activity of lattice oxygen species

To support our argument, we used the independent gradient model (IGM) to visually analyze the non-covalent interactions between molecular fragments and active sites to qualitatively show the different environmental effects before and after propane molecules' activation by bridging oxygen (Mo–Oa–Mo) and terminal oxygen (MoOt) (Fig. 2b). The results showed that when propane was close to MoOt, strong attraction (blue) and weak van der Waals force (green) were dominant, while strong repulsion (red) appeared at Mo–Oa–Mo. Furthermore, when comparing the transition states, MoOt underwent no evident reconstruction. However, there was obvious reconstruction in Mo–Oa–Mo, and the previous strong repulsion disappears completely in the reconstructed transition state structure. Strong attraction indicated that the C–H bond was about to break, and meanwhile the O–H bond was about to form, which, together with van der Waals force, constitutes stabilization.

Moreover, we used the crossing potential model to deeply understand how the environmental structure affected the C–H activation energy (Fig. 2c). The model showed that the potential energy varies with the distance from the adsorbate to the surface. It is generally believed that the activation of reactants would undergo a transition from physisorption to chemisorption.[25] The crossing point of the two potential functions is exactly the transition state. The results showed that the C–H bond activation on the two active sites had almost the same O–H bond distance (∼0.98 Å). However, the energy of OaH was lower than that of OtH, which was the direct embodiment of the stronger hydrogen affinity of Oa than that of Ot. Mo–Oa–Mo is more active than MoOt intrinsically. As the hydrogen affinity becomes exothermic, the chemisorption potential moves down as a whole, and meanwhile, the crossing point (transition state) also moves down. Therefore, there is a linear relationship between adsorption and activation energy (i.e., Brønsted–Evans–Polanyi relationship).[26]

However, from the stable physisorption to the transition state, propane continues to approach the surface, causing some overlap of the electron clouds on the surface of propane and the catalyst. It will lead to repulsive interaction (i.e., Born repulsion), which induces surface reconstruction by increasing the potential energy. Comparing the energy changes of both active sites, we found that the energy rise of Mo–Oa–Mo is more obvious than that of MoOt. Combined with the results of IGM analysis, we can know that this stage is mainly due to the strong surface reconstruction in the process of C–H bond activation at the bridging oxygen (Mo–Oa–Mo).

Optimizing the oxygen activity by tuning the environmental structure

Such strong surface reconstruction revealed that the unfavorable environmental structure near α-MoO3 (010) Mo–Oa–Mo originates from the concave local environment of Mo–O–Mo (Fig. 3a). Therefore, it is likely to optimize the steric spatial impact and enhance the catalytic activity to a certain extent. We noted that the molybdenum-based Keggin-polyoxometalate (POM) is an ideal model catalyst to verify our idea. Mo-based POMs (HXMo12O40, X is the central heteroatom, such as P, Si, and Al) are kinds of metal oxide clusters with well-defined chemical composition and structure. POM has excellent thermal stability and redox properties, and there are many remarkable experimental and theoretical studies on the oxidative dehydrogenation of methanol[42,43] and low-carbon alkanes.[44] More importantly, POM has an optimized local environment (convex) of surface lattice oxygen compared to MoO3 and other comparable structural features (Fig. 3b and Table S3†). Phosphomolybdic acid hydrate (H3PMo12O40) is the most common Mo-based POM. Therefore, it is very suitable as a model comparison catalyst. Similar to α-MoO3, there are two kinds of oxygen species (MoO and Mo–O–Mo) on POM. We calculated the EH of different oxygen species according to their symmetries (Fig. S5†). The results suggested that all the bridging oxygens (Mo–O–Mo) are stronger than terminal oxygens (MoO) in hydrogen affinity like α-MoO3. We selected the terminal oxygen and bridging oxygen closest to the average EH as the representative of the average reactivity for the following investigation. Similarly, we compared the C–H activation energy for the first step dehydrogenation of methane, ethane, and propane on different oxygen species on α-MoO3 (Fig. 3c) and H3PMo12O40 (Fig. 3d). The results showed that in contrast to α-MoO3, all the C–H activation energies of Mo–O–Mo are lower than those of MoO. The results showed that the C–H bond activation activity of the bridging oxygen with convex curvature was higher than that of the bridging oxygen with concave curvature.

Fig. 3

Design of more active catalysts by optimizing the microenvironment around the oxygen species. (a and b) Optimized structures of α-MoO3 (010) pristine surface and POM (H3PMo12O40) with concave and convex curvature structure characteristics, respectively. H atoms are hidden; (c and d) activation energy of CH4, C2H6, and C3H8 at MoO and Mo–O–Mo active sites with concave and convex curvature environmental structures. (e and f) Correlation between activation deformation energy and C–H activation energy of alkanes at the active sites of oxygen species with different structures. The red labels indicate the key active sites.

Moreover, to quantitatively describe the surface reconstruction induced by adsorbates mentioned above, we propose a new descriptor, the activation deformation energy (ADE), defined using the following equation:ADE = where TS and IS are short for the transition state and initial state, respectively, ADE is single-point energy with clear physical meaning, the energy increase caused by the surface reconstruction induced by the activated adsorbates. We correlated the C–H activation energies of the terminal oxygen and bridging oxygen on the surface of α-MoO3 and H3PMo12O40 with ADE (Fig. 3e and f) and found that there was a good linear relationship.

For α-MoO3, ADE (0.9–1.2 eV) on the Mo–O–Mo active site is significantly larger than that on MoO (0.4–0.7 eV), indicating that there is a strong reconstruction on the bridging oxygen, but not on the terminal oxygen; For POM, the surface reconstruction is relatively weak (0.48–0.75 eV) for both terminal oxygen and bridging oxygen. It is worth noting that the difference between ADEs of the same reactant at different active sites is almost the same. For example, it is 0.50 eV for α-MoO3 and 0.11 eV for H3PMo12O40. When comparing EH and the constant difference, we found that for α-MoO3, EH could not cover the difference, but for H3PMo12O40, it could. The surface reconstruction induced by H atom adsorption is almost negligible due to its tiny volume. However, the real reactants may induce strong surface reconstruction in the transition state, and such information is not included in the EH. Therefore, for some potential active sites (e.g., Mo–O–Mo on α-MoO3), the obvious surface reconstruction induced by adsorbates is the fundamental reason why EH cannot give an accurate prediction.

Experimental verification

A thermogravimetric test (Fig. S6†) was done on the H3PMo12O40/Al2O3 catalyst to ensure consistency between the theoretical model and experimental structure. Then 723 K was determined to be the appropriate reaction temperature.

To investigate the structural stability of the theoretical model at the reaction temperature and the changes in the microenvironmental features of Mo–O–Mo, we did ab initio molecular dynamics simulation (AIMD) for 30 000 fs at 723 K for H3PMo12O40 (Fig. S7†). The results showed that the overall structure and energy are stable. Besides, the statistical average bond angle about the active site (Mo–O–Mo) is 127.5°, and its key convex curvature structural features remained unchanged.

To verify the theoretical results, we tested the initial performance of two Mo-based oxide model catalysts for propane dehydrogenation at 723 K (Fig. 4a). The propane conversion on H3PMo12O40/Al2O3 was nearly five times higher than on α-MoO3/Al2O3 without significantly decreasing propylene selectivity. The carbon balance was more than 98%, indicating that there was little coking. The results showed that the yield of propane dehydrogenation catalyzed by supported H3PMo12O40 was significantly better than that of α-MoO3 under the same test conditions.

Fig. 4

Experimental verification. (a) The instantaneous performance at the 5th min of 20 wt% α-MoO3 physically dispersed on the Al2O3 support and 20 wt% POM(H3PMo12O40)/Al2O3. Reaction conditions: T = 723 K, 1.4 atm (C3H8 : N2 = 0.125) and 0.5 g sample. (b) DFT calculated rate constant of RDS at active sites (MoO for MoO3; Mo–O–Mo for H3PMo12O40) and experimentally measured productivity. (c) TEM and XRD for α-MoO3 (d) TEM for 20 wt% H3PMo12O40/Al2O3. (e) SAED for α-MoO3 sample. (f) EDS elemental mapping for 20 wt% H3PMo12O40/Al2O3.

Furthermore, to better explain the difference in the intrinsic activity of the active sites in the catalyst experiments, we calculated the activation free energy as well as the rate constants of the rate-determining step (RDS) for the active sites on the two model catalysts (Tables 1 and S4†). We found that the ratio of the rate constants of the RDS at the active sites of the two model catalysts and the experimentally observed productivity ratios are consistent in the order of magnitude (Fig. 4b).

Comparison of Gibbs's activation free energy (ΔG≠) and rate constant (k) at different active sites

Catalyst	Site	ΔG≠a/eV	k(RDS)b/(s−1)
α-MoO3	MoO	2.08	2.35 × 10−2
H3PMo12O40	Mo–O–Mo	2.00	8.72 × 10−2

Temperature: 723 K.

RDS: C3H8(g) → C3H7(g) + H*.

TEM and XRD (X-ray diffraction) (Fig. 4c), and SAED (selected area electron diffraction) (Fig. 4e) confirmed the rationality of our theoretical model, that is, α-MoO3 was a lamellar structure and the (010) facet was the dominating plane. TEM (Fig. 4d) and energy-dispersive X-ray spectroscopy (EDS, Fig. 4f) for 20 wt% H3PMo12O40/Al2O3 showed that POM was uniformly dispersed on the alumina support.

Conclusions

Our work showed that the C–H bond activation of alkanes could be significantly enhanced by tuning the spatial structure around the surface lattice oxygen of molybdenum oxides. Firstly, we explored the differences in the activity of different surface lattice oxygens and their nature using α-MoO3 (010). The activity of terminal oxygen (MoO) is higher than that of bridging oxygen (Mo–O–Mo) on the α-MoO3 (010) plane. However, the hydrogen affinity (EH) showed that Mo–O–Mo is the potential active site with higher intrinsic activity but limited by its unfavorable environmental structure with concave curvature.

Furthermore, we designed a cluster-type molybdenum oxide (i.e., Keggin-type polyoxometalate, H3PMo12O40) with well-defined structures with concave curvature for Mo–O–Mo. We proved that such Mo–O–Mo would significantly enhance the C–H bond activation of low-carbon alkanes compared to α-MoO3 (010) by density functional theory calculations. Moreover, we found that a universal descriptor of hydrogen affinity (EH) cannot predict the real activity of the oxygen species with concave curvature. Therefore, a quantitative descriptor, activation deformation energy (ADE), is proposed to describe the surface reconstruction induced by adsorbates. ADE provides us with a quantitative method with clear physical meaning to help deeply understand how the influence of adsorbate-induced surface reconstruction affects the activation energy. Also, we found that EH can give the correct qualitative activity prediction about lattice oxygen only when EH is greater than ADE; otherwise, it cannot.

Finally, experimentally, we used alumina-supported H3PMo12O40 and α-MoO3 as contrasting catalysts to test the propane dehydrogenation performance and confirmed the obvious improvement (∼5 times) of propane activity and propylene productivity when using H3PMo12O40/Al2O3.

This work provided more insights into the influence of the environmental structure around lattice oxygen species on metal oxides and some guidance for designing highly active metal oxide catalysts by tuning the surface environment around lattice oxygen.

Methods

The calculations were performed using the Vienna ab initio simulation program[27,28] (VASP), version 5.4.4 using the projector augmented wave (PAW) method.[29-31] Self-consistently calculated total energies were evaluated with the Perdew–Burke–Ernzerhof (PBE)[32] functional. The Kohn–Sham equations were solved using a plane-wave basis set with a cutoff energy of 400 eV. The Gaussian smearing method was adopted with a width of 0.05 eV to determine how each wave function's partial occupancies are set. Dipole correction was considered where necessary. All calculations were spin-polarized unless otherwise stated.

For systems such as transition metal oxides, the electronic bandgap may be underestimated, and even a qualitatively wrong metallic ground state may be predicted. The DFT + U approach[33,34] has been proved to study a large variety of strongly correlated compounds with considerable improvement concerning LDA or GGA results. The Hubbard U value was applied to the molybdenum centers in this work. By fitting the enthalpies of formation and reduction of MoO3, we determined that the appropriate Ueff (where Ueff = U − J) value is 2.0 eV, which is consistent with a previous report.[20] The experimental data of enthalpy of formation and enthalpy of reduction were obtained from the NIST-JANAF database (Fig. S2†).

For bulk optimization, the sampling of the Brillouin zone was performed using a Monkhorst–Pack scheme[35] of (4 × 1 × 4). The convergence of the k-point is confirmed with the H adsorption energy. All atoms in the one-bilayered slabs for 2 × 2 to 5 × 5 supercells were relaxed. Moreover, as a descriptor, oxygen vacancy formation energy was used to test the convergence and confirm that the appropriate model in the XY direction is 4 × 4.

For the Z direction, adsorption energies of some important intermediates on one-bilayered 2 × 2 and two-bilayered 2 × 2 slabs were nearly identical, with differences of less than 0.05 eV. Thus the one-bilayered slab was used to simulate the Z direction. Slabs were separated from their periodic images with 15 Å of vacuum (Fig. S3†).

All calculations are performed on the Mo-based Keggin POM (H3PMo12O40) clusters (diameter ∼1.1 nm) placed in the center of 2 × 2 × 2 nm3 cells to prevent the electron density overlap between clusters in adjacent cells. Gamma point (1 × 1 × 1) mesh was used to sample the first Brillouin zone. Dipole/quadrupole corrections were used for H3PMo12O40 to eliminate long-range electrostatic interactions and thus allow simulation of isolated clusters within periodic boundaries of VASP. The structures of reactants, products and stable intermediates were optimized until the atomic force was less than 0.05 eV Å−1. The other calculation details are consistent with MoO3.

Transition state geometries were located using the climbing image nudged elastic band method (CI-NEB)[36,37] with total 6 images, followed by refinement with the dimer method.[38] Vibrational frequencies, computed from the Hessian matrix under the harmonic approximation, were used to calculate zero-point energy (ZPE) corrections. Calculating the frequency of all transition state structures confirmed that they have only one imaginary frequency.

Independent gradient model (IGM) is a new electron density (ED)-based methodology developed by Lefebvre et al. to identify different types and strength of interactions between chemical fragments, especially weak non-covalent interactions (NCI), such as vdW interaction and hydrogen bonds.[39] The expression of the ED gradient in terms of atomic components furnishes the basis for the Independent Gradient Model (IGM). Using an intra/inter uncoupling scheme, a descriptor (δginter) is then derived that uniquely defines intermolecular interaction regions. For a detailed definition and discussion, please refer to the original ref. 39 All IGM analyses were carried out using Multiwfn software (version 3.8).[40] The visualization of the 3D isosurface (isovalue = 0.01 a.u.) is realized by VMD software.[41]

Data availability

The data that supports the findings of this study is available from the corresponding author upon reasonable request.

Author contributions

J. G. conceived and coordinated the research. C. J. conducted the density functional theory calculations and wrote the draft. X. W. contributed to catalyst synthesis, catalytic experiments and characterization. C. J, X. C, and Z. Z. analysed the data. C. J., X. C., Z. Z. and J. G. wrote the manuscript. All authors discussed, commented on and revised the manuscript.

Conflicts of interest

There are no conflicts to declare.