Density functional theory study on the catalytic dehydrogenation of methane on MoO 3 (0 1 0) surface

.


Introduction
Prices fluctuations, the instability of oil and gas production, and the consequences of the COVID-19 outbreak, have motivated the search for transforming methane into value-added chemicals that have higher economical and strategic importance.Furthermore, the expansions in production, rising demand for hydrogen fuel, and global shortage in aromatic compounds and olefins are additional motivations for methane transformation into chemical commodity products [1][2][3].The shift from conventional oil into shale oil has also affected the catalytic production of some hydrocarbons, such as olefins and small aromatics such as ethylene, benzene, and naphthalene.These chemicals are essential for manufacturing polymers, fertilizers, pharmaceuticals, and many other products.This becomes valuable once we know that the global production and demand of natural gas will be doubled by 2045 [4].More importantly, non-oxidative methane conversion yields hydrogen gas, which is considered as the fuel-of-the-future [2,3,5].
The direct conversion of methane to other chemicals is challenging, mainly due to the strength of the C -H bond (439 kJ/mol), and its low polarizability, which makes it very resistive to chemical attacks [1,6,7].The C-H bond breaking is believed to be the rate-limiting step in the process of methane activation [8].Therefore, both experimental and theoretical studies in this field have focused on activating the C-H bond to improve the selectivity and efficiency of methane conversion.It is widely agreed that methane activation should be catalyzed using materials such as molybdenum/zeolites [7,9,10].Nevertheless, improving the efficiency of the process requires better understanding of the catalyst surface, active sites, their acidity towards the C -H bond, as well as the reaction mechanism and kinetics.
Since the early 1990 ′ s, theoretical quantum calculations have led the progress in methane conversion [3,5,[11][12][13][14].These calculations were facilitated by studying reactions which utilize clusters of zeolites, metal oxides, or a hybrid model.However, it is more appropriate to study methane conversion under periodic boundary conditions as it provides better description of the geometries, lattice relaxation, and provides realistic values for the reaction energies [15].In this work, we focused our attention on methane interactions with the surface of neat molybdenum oxide, MoO 3 , in a fundamental effort to understand the methane activation at the molecular level.The results of this study will serve as a benchmark for future experiments on Mo-modified zeolites.
MoO 3 is a crystalline solid that has an orthorhombic structure.The material has several industrial applications in the catalytic partial oxidation of alkanes [8], biomass upgrading [16], reduction of nitrogen oxides (NO x ) [17], and pesticides destruction [18].It is believed that methane activation is initiated by a hydrogen abstraction on the Lewis acid sites (O atoms) on MoO 3 surface [17].This step was proven to take place on the MoO 3 (0 1 0) but not the (1 0 0) surfaces [19,20].The selectivity of the activation might also be influenced by modifying of the acidity of the catalytic sites through doping with other transition elements, such as Fe and Co [2,3,9].
In this work, we employed periodic density functional theory calculations to study the adsorption of methane on MoO 3 surface.Two different (4 × 4 × 4) supercells of the α-MoO 3 (0 1 0) surface were constructed to study the adsorption.We also investigated the reaction mechanism initiated by a C -H rupture in methane to produce hydrogen gas.While previous works targeted methane dehydrogenation using gaseous clusters [11,21,22], the originality of this work lies in locating the transition states of the reaction under periodic boundary conditions.In addition to hydrogen formation, the mechanism showed the formation of primary radicals (e.g., • CH 3 ) that can further react to form small olefins or aromatics.Ethylene was also a main product of this mechanism.

Computational details
This study was performed on α-MoO 3 , which is the most stable crystalline phase of molybdenum oxide [16,18].The structure is orthorhombic that belongs to the Pbmn space group with experimental parameters of a = 3.96, b = 13.86, and c = 3.69 Å [23,24].It is believed that the crystal structure of MoO 3 forms bilayers that are parallel to its (0 1 0) surface, which is MoO 3 ′ s most thermodynamically stable geometry [17][18][19].The density functional theory calculations of this work were achieved using the BAND engine in Amsterdam Modeling Suite (AMS) 2021 [25][26][27].The Perdew − Burke − Ernzerhof (PBE) exchange correlation functional with Grimme's dispersion (PBE-D) [28,29] was implemented with the double-ζ basis sets (DZP).The method was selected based on a benchmark test (shown in Table 1) as well as its previous performance in similar studies [8,16,18,30].Because of the presence of Mo atom (a heavy metal with Z = 42), scalar relativistic effects were treated with the accurate and efficient ZORA approach [30,31].Also, the smallest frozen core approximation was adopted as advised [26,30].All geometry optimizations were completed with an energy cut-off of 0.0001 Ha and a gradient convergence of 0.005 Ha/Å.A zero-spin polarization (corresponding to singlet multiplicity) was assumed for all closed shell spins.Otherwise, the spin polarization (spin) was calculated as the (2S) component of the spin-orbit interaction in open shell molecules.For instance, the spin of • CH 3 would be 1 (doublet) and O 2 is 2 (triplet).

Table 1
Relative percentage errors associated with optimizing the MoO 3 unit cell using different functionals.Lattice optimization was done using the DZP basis set.To validate the theoretical method used in the current work, the MoO 3 unit cell, obtained from an experimental cif file [23], was optimized at different levels as shown in Table 1.The optimization was performed by selecting the (optimize lattice) keyword in AMS using different set of functionals and the DZP basis set.An optimized structure of the unit cell is shown in Fig. 1.The best results were obtained using the PBE-D, PBE-D3, and the Minnesota M06L functionals.Thus, we decided to adapt the PBE-D functional for the rest of this work, as it includes the dispersion component needed for methane interactions with the MoO 3 surface, in addition to its reasonable computational cost.
To study the adsorption of methane on MoO 3 surface, two supercells were constructed as shown in Fig. 2. The first supercell, SC1, was modeled from the 010 surface and it includes 4 × 4 × 4 layers.This produces a formula of Mo 36 O 108 , which maintains the empirical ratio of the original unit cell.SC1 is constructed so the O atoms are dangling in the 001 (z) direction.The second supercell, SC2, on the other side, was constructed from the same surface but with Mo and O atoms constructing a smooth surface without any dangling O atoms, as shown in Fig. 2-b.Having these two different constructions of the (0 1 0) supercells allowed us to investigate different adsorption behavior on the surface.
The methane adsorption was studied on different sites as shown in Fig. 3-a.In the top position, the CH 4 molecule is placed right above the Mo atom.In the bridge position, CH 4 is situated between two oxygen atoms, and in the hollow position, the methane molecule is situated in the cavity formed between four oxygen atoms.Similar sites are defined for SC2 as illustrated in Fig. 3-b.The adsorption was performed by relaxing the methane molecule and the upper layer while freezing the other core layers.The adsorption energy (E ad ) was thus calculated as described in equation ( 1) [17]: where E adsorbate/surface is the total energy of the optimized slab/CH 4 system, E surface is the total energy of the optimized Mo 36 O 108 slab, and E adsorbate is the total energy of optimized methane.The definition of eq. ( 1) implies that a negative E ad responds to stronger adsorption [16,17,30].It is noteworthy that eq. ( 1) does not account for changes in zero-point energies (ZPE) or heat capacities because of the high possibility of having imaginary frequencies in such periodic calculations [8,24,30].
For the reaction mechanism of CH 4 activation on MoO 3 surface, a simple surface cluster (Mo 4 O 12 ) was used, and the calculations were achieved using the PBE-D functional.The triple-ζ polarized (TZP) basis set was used to obtain trusted values of energy changes.Transition states involved in this work were obtained by exploring the potential energy surface (PES) along the reaction coordinate of interest, followed by optimizing the highest energy point into a transition state [32,33].Each transition state was confirmed to have only one imaginary frequency along the reaction coordinate.Furthermore, the transition states were validated by constructing an intrinsic reaction coordinate (IRC) at the same level of theory and ensured that it connects the desired reactants to their products.
To validate the calculation method used in the reaction mechanism, a benchmark study was performed on simple molecules, including water, ammonia, and methane, to estimate their bond dissociation energies, and comparing the results with reported experimental values.The calculations were done by optimizing the species and fractions of each molecule at the PBE-D/TZP level of theory, followed by frequency calculations as explained earlier.The bond dissociation energy (BDE) is defined as the change of enthalpy at 298 K (ΔH 298 ) for the reaction showed in equation ( 2) [34]: The enthalpy change (ΔH 298 ) for a given species can be computed using ab-initio techniques as: where ΔH 0 is the enthalpy change at 0 K and is defined as: where SPE is the single point energy and ZPE is the zero-point energy.
The enthalpy correction to enthalpy (H corr. ) was obtained from frequency calculations after successful optimization for each specie.The results of the benchmark study are shown in Table S1 (Supplementary Information Section).The level of theory used in work demonstrates good agreement with the BDE experimental values.The total average error from all molecules under investigation is estimated at 2.6%.

Methane adsorption on MoO 3 (0 1 0) surface
The adsorption of methane on the MoO 3 supercell (SC1) was studied by placing the methane molecule next to the surface in the bridge, hollow, and top positions.Then, the calculations started by relaxing both the CH 4 and the upper layer, while freezing the bulk layers of MoO 3 .Fig. 4-a shows an optimized geometry of methane in the bridge position at PBE-D/DZP level of theory.Similar geometries for the hollow and top positions are shown in Fig. 4-b and c, respectively.Table 2 lists some selected bond distances between the methane atoms and the surface after optimization.The final geometries show electrostatic interactions between H and O atoms in MoO 3 at distances between 2.3 and 2.7 Å, which lie in the range of hydrogen bond interactions [35].It was reported that small molecules such as O 2 , H 2 O, and CH 4 form weak Van Der Waals interactions with MoO 3 surface [18,36,37].The calculated E ads values on the SC1 slab range between 168 and 179 kJ/mol (1.75 and 1.86 eV).The positive values indicate unfavorable adsorption as per eq.( 1).The E ads

Table 2
Adsorption energies (E ads ) and selected bond distances for methane adsorption on α-MoO 3 (0 1 0) surface obtained at PBE-D/DZP level of theory.The adsorption of methane on the SC2 supercell is shown in Fig. 5.To the best of our knowledge, this surface has never been considered for methane adsorption in the literature.Under this new geometry, O and Mo atoms form a smooth surface with no dangling O bonds.The adsorption energies on SC2 are found to be − 12.7 and − 10.8 kJ/mol for the bridge and hollow positions, respectively.The reaction is exothermic indicating favorable adsorption as defined in eq. ( 1).As shown in Table 2, the O atoms in the pure SC2 bear charges of − 0.313 and − 0.342.After adsorption, the charges slightly decreased to − 0.293 and − 0.340.The C -H polarization was not very different from that of SC1.But the remarkable difference between the adsorption on SC1 and SC2 is in the C-Mo bond distance.In SC2, the distance is as low as 3.5 Å, which is lower than that in SC1.Also, the charges on O atoms in SC2 are more negative.Thus, we conclude that the O δ--C δ-repulsions are more important than O δ--H δ+ attractions in methane adsorption.In other words, the smooth surface of SC2 can better accommodate the methane molecule by avoiding too many O δ--C δ-repulsions.
The adsorption on the top position in SC2 comes with a high unfavorable E ads (201.2 kJ/mol).To explain this finding, the bond distances and atom charges of the final optimized CH 4 /SC2 system are not very different from the other positions.However, a close look at Fig. 5-c shows that methane molecule is placed right above the Mo atom, resulting in a large electrostatic C δ+ -Mo δ+ repulsion leading to the methane molecule being pushed away from the surface.The repulsion is reflected by the long C-Mo distance in SC2 (4.93 Å), which is almost one angstrom longer than the corresponding bridge and hollow positions.Consequently, the final CH 4 /SC2 system has an additional factor of instability, resulting in large E ads for that position.

Methane dehydrogenation on MoO 3 (0 1 0) surface
The reactions of CH 4 with the α-MoO 3 (0 1 0) surface were studied at the PBE-D/TZP level as explained in Section 2 (computational details).A smaller, simplified slab from the overall structure was used for this part of the study to ease the calculations and maintaining periodic boundary calculations and is shown in Fig. 6-a.It was previously claimed that O atoms coordinated to Mo in MoO 3 are the active sites capable of abstracting H from methane [38].Therefore, the Mo 4 O 12 slab was selected with O atoms dangling away of the surface.Fig. 7 depicts the detailed reaction mechanism initiated by Habstraction by O atom.A transition state for this step was located (TS1) and its optimized geometry is shown in Fig. 6-b.The TS has an imaginary frequency of − 1260i cm − 1 that corresponds to a C -H stretching towards the O atom.This primary step, which is considered the rate limiting step, as mentioned in the introduction, leads to the formation of a radical intermediate (M1) with an O-H terminal, in addition to a reactive • CH 3 radical.An Intrinsic reaction coordinate (IRC) for the TS1 that connects the reactant pair to methyl radical and M1 is shown in Fig. S1 (Supplementary material).The spin polarization (spin) of methyl  radicals is 1 (doublet).Considering the principle of spin conservation [39,40], the total spin of these two separate radicals could be equal to 0 (singlet) or 2 (triplet).Although the ground state of MoO 3 is singlet [41], TS1 was only located with a spin equal to 2. All our efforts to locate a singlet transition state for this step have failed.Thus, the H-abstraction step using our calculation method can only take place at a triplet surface.The debate on abstraction reactions at a triplet surface is long standing.In the early 1990 ′ s, a two-state reactivity was developed by Shaik [39,42], Schwarz [43,44], Schroder [45] and others [41,46,47] to describe reactions catalyzed by organometallic systems.Under this approach, reactions can proceed on different electronic states through spin inversion or conical intersection [47,48].The theory was first applied to metal ion species (MO + ), but it is now extended to describe neutral systems.MoO 3 is one example of transition metals compounds that often possess low-spin excited states that are adjacent to their ground states.Thus, their reactivity may involve two or more electronic states [39].Yin et al., reported that high spin states in metal oxides become more accessible as a result of increased exchange stabilization of nonbonding (d) orbitals [41].Consequently, hydrogen abstraction  pathways can happen with higher multiplicities through thermal or photoexcited Mo-O activated complexes [38,49].
Under the level of theory of this work, the energy difference between the singlet Mo 4 O 12 ground-state and its first triplet excited state is 160 kJ/ mol (1.66 eV).The singlet M2 intermediate (cf.Fig. 7) is only − 52 kJ/ mol (0.54 eV) below its triplet state.These values suggest that a spin inversion can cause the H abstraction, which could take place either thermally or photochemically, and the reaction will be similar to a H abstraction from methane with triplet O 2 [50].However, such abstraction would be difficult to occur under ambient conditions, which requires further studies of modified MoO 3 catalyst, such as zeolites.
The activation barrier for the H-abstraction through TS1 is equal to 66.4 kJ/mol as shown in Fig. 8.This barrier is much lower than the one obtained by Zhang et al., at 217 kJ/mol [8].Also, the product pair, M1 and • CH 3 are 97.8 kJ/mol above the reactants, compared to 68 kJ/mol in Zhang's work.This difference is mostly due to the small basis set that was used.Fu et al., estimated the E a 's for H-abstraction from methane by Mo 2 O 6 clusters to be in the range of 190-270 kJ/mol [38].Although these values were obtained at a singlet surface, the difference from our value is because Fu's work was accomplished for separate clusters in the gas phase and no periodic conditions were applied.The reported C -H bond dissociation energy in methane is 439 kJ/mol [51].This is almost 6-folds higher than the first activation barrier (66.4 kJ/mol) to break the C-H bond as per our calculations, reflecting the catalytic ability for the MoO 3 surface to activate methane.A future investigation is needed to confirm this fact experimentally.
Once M1 is formed, a fast combination between the dangling O atom and • CH 3 can take place as previously suggested for Mo 2 O 6 clusters [38,41].Alternatively, • CH 3 self-recombination will form ethane gas.The competition between these two pathways depends on the availably of the O active sites and the reaction environment [8,11,41,45].The attack of • CH 3 on M1 will lead to the formation of the metastable methoxy intermediate, M2, as shown in Fig. 7.We have not found any transition state for this reaction as free radical recombination is expected to proceed with no barrier [52,53].The reaction can then proceed from M2 through a H 2 elimination step as shown in Fig. 7. First, M2 can undergo electronic relaxation from its triplet state to its singlet ground state with an estimated difference of 8.8 kJ/mol.Then, the hydrogen elimination takes place in one concerted step to form the O--CH 2 intermediate (M3).The transition state for this reaction, TS2, was located and it has an imaginary frequency of − 492i.Optimized geometries of TS2 and the other species involved in this mechanism are shown in Fig. 9.The E a for this step is 88.8 kJ/mol and ΔH is − 12.6, as depicted in Fig. 10-a.The barrier is relatively low, and the reaction is slightly exothermic which makes hydrogen formation through this step quite favorable.
The intermediate M3 is relatively stable with a spin of 0. So, it can abstract another • CH 3 radical to form M4. The transition state for this reaction, TS3, has a spin of 1 as it links two free-radical systems.TS3 has an imaginary frequency of − 665i corresponding to C -C bond stretching as expected.Interestingly, this step is almost barrierless (− 1.7 kJ/mol) and the reaction is highly exothermic (-281.4kJ/mol, Fig. 10-b), which makes the formation of M4 highly feasible.The latter is an ethoxy radical with a spin of 1, that can undergo a 1,3H-shift to the adjacent O atom to form ethylene and a new intermediate, M5.In fact, M5 and M1 are two different isomers for the same intermediate because they have the same periodic structure with different position of the H atom.The potential energy diagram for this step is shown in Fig. 10-c.The barrier is estimated at 128.6 kJ/mol, and the reaction is endothermic by 117.6 kJ/mol.
Based on the above results, we suggest that methane dehydrogenation with the MoO 3 surface can be initiated by H-abstraction, followed by free-radical attacks to mainly form • CH 3 , H 2 , and ethylene.Depending on the reaction environment, further reactions involving the • CH 3 and • CH 2 CH 3 radicals can proceed to form higher hydrocarbons or aromatics.

Conclusion
The purpose of this study was two-fold; namely: exploring the adsorption of CH 4 on MoO 3 (0 1 0) surface and proposing a mechanism for hydrogen formation through CH 4 dehydrogenation reaction.Through initial benchmark studies, the most fitting level of theory for conducting the adsorption studies was found to be PBE-D/DZP.The final geometries obtained from the adsorption studies revealed interactions in the hydrogen bond range.The hollow position in SC1 supercell was found to be the most favorable site for adsorption as it had the lowest E ads value.In SC2, the bridge and hollow positions yielded favorable energies.The adsorption experiments led to the conclusion that O δ--C δ- repulsions have more significant role in CH 4 adsorption on the catalyst surface than O δ--H δ+ attractions.
The proposed reaction mechanism for the CH 4 dehydrogenation on the MoO 3 (0 1 0) surface has revealed new pathways to form hydrogen gas, ethylene molecule, and methyl radicals.Such species can undergo further secondary reactions to form higher molecular weight molecules.The novelty of our mechanism lies in the higher basis set used which consequently led to more reliable activation barriers.The mechanism is also significant as it employed periodic conditions and the dehydrogenation was conducted on a slab as opposed to separate clusters.This study lays the foundation for future research into the dehydrogenation of CH 4 on catalytic surfaces using higher basis sets as well as detailed work related to adsorption studies and possible doping effects on molybdian/zeolite-based catalysts.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 6 .
Fig. 6. a) Optimized α-MoO 3 surface used to study the reaction mechanism, optimized geometries of TS1 (b), and the intermediate M1 (c) obtained at PBE-D/TZP level of theory.

Fig. 8 .
Fig. 8. Potential energy level diagram for methyl radical formation when methane reacts with α-MoO 3 surface obtained at PBE-D/TZP level of theory (ZPE corrections are not included).

Fig. 10 .
Fig. 10.Potential energy level diagram for methane reactions with α-MoO 3 surface obtained at PBE-D/TZP level of theory (ZPE corrections are not included).
This research is funded by Qatar University's grant number QUCG-CAS-21221.CRediT authorship contribution statement I. Badran: Conceptualization, Data curation, Validation, Writingreview & editing, All authors reviewed the results and approved the final version of the manuscript.N. Riyaz: Data curation, Validation, Writing review & editing, All authors reviewed the results and approved the final version of the manuscript.N. A. Shraim: Validation, writing, review and editing, All authors reviewed the results and approved the final version of the manuscript.N. Nassar: Conceptualization, writing, review and editing, All authors reviewed the results and approved the final version of the manuscript.