Metal dimer sites in ZSM-5 zeolite for methane-to-methanol conversion from first-principles kinetic modelling: is the [Cu–O–Cu]2+ motif relevant for Ni, Co, Fe, Ag, and Au?

Direct methane-to-methanol conversion is a desired process whereby natural gas is transformed into an energy-rich liquid. It has been realised at ambient pressure and temperature in metal ion-exchanged zeolites, where especially copper-exchanged ZSM-5 has shown promising results. The nature of the active sites in these systems is, however, still under debate. The activity has been assigned to a [Cu-O-Cu]2+ motif. One remaining question is whether this motif is general and also active in other metal-exchanged zeolites. Herein, we use first-principles microkinetic modelling to analyse the methane-to-methanol reaction on the [Cu-O-Cu]2+ motif, for Cu and other metals. First, we identify the cluster model size needed to accurately describe the dimer motif. Starting from the [Cu-O-Cu]2+ site, the metal ions are then systematically substituted with Ni, Co, Fe, Ag and Au. The results show that activation of Ag and Au dimer sites with oxygen is endothermic and therefore unlikely, whereas for Cu, Ni, Co and Fe, the activation is possible under realistic conditions. According to the kinetic simulations, however, the dimer motif is a plausible candidate for the active site for Cu only. For Ni, Co and Fe, close-to-infinite reaction times or unreasonably high temperatures are required for sufficient methane conversion. As Ni-, Co- and Fe-exchanged ZSM-5 are known to convert methane to methanol, these results indicate that the Cu-based dimer motif is not an appropriate model system for these metals.


Introduction
The known reserves of natural gas at the end of 2014 were 190 trillion cubic meters, which are enough to sustain global consumption for more than half a century. 1 Unfortunately, much of the reserves are located in isolated areas, like offshore oil wells, and as a consequence a large amount of this valuable resource is today flared into carbon dioxide and water. It is estimated that 3.5%, or about 143 billion cubic metres, of the natural gas globally extracted was flared in 2012. 2 Natural gas consists mainly of methane along with minor amounts of higher alkanes, such as ethane, propane and buthane. 3 Owing to the low carbon emission per energy unit of methane, natural gas is often regarded as a transitional energy source while the energy systems of society are trans-formed to more sustainable solutions. 4 Furthermore, methane is an important raw material for other upgraded products. 5,6 Hence, methods to convert methane into a liquid, e.g., methanol, are highly desirable, as these would enable easier transportation and thereby allow this valuable resource to become accessible to the chemical industry.
At present, industrial conversion of methane to methanol utilises a two-step process, where methane (natural gas) is transformed into synthesis gas using steam reforming, 7 which is later used in methanol synthesis. 8 These energy-intensive processes require large centralised plants, which are poorly compatible with the distribution of methane reserves. Direct conversion of methane to methanol would be a considerably more energy-efficient route. However, this reaction is associated with intriguing challenges and has so far not been industrially realised. 6,9 One of the main challenges is to activate the highly stable methane for oxidation but still prevent complete oxidation of the reactants, which is thermodynamically favorable. 10 Another issue is the activation of molecular oxygen and the formation of a reactive oxygen species that methane can react with. In short, it is a selectivity problem where methane should be oxidised, but not over-oxidised. Ion-exchanged zeolites with different ions (e.g. Cu, Ni, Co and Fe) have been shown to convert methane to methanol under ambient conditions. [11][12][13][14][15][16] One of the most studied systems for this reaction is copper-exchanged ZSM-5. 11,[16][17][18][19][20][21][22][23][24][25][26] The first experimental observation of successful methane-tomethanol conversion under ambient conditions, i.e., atmospheric pressure and a temperature of 448 K, for Cu-ZSM-5 was reported by Groothaert et al. 11 However, the experiment required O 2 activation at a much higher temperature (723 K) and the introduction of a water/acetonitrile mixture for the sequential methanol extraction, which makes the process industrially unattractive. Recent theoretical studies 27 have indicated that water reduces the methanol desorption barrier, which could be the underlying reason for adding an extraction solvent. Although a continuous cycle has recently been successfully demonstrated, 28,29 an efficient catalyst for methane-to-methanol conversion is still lacking. For progress, it is desirable to scrutinise the active site(s), to understand the oxygen activation and to clarify the influence of extraction-solvation molecules like water in these systems.
There have been many attempts to determine the active sites in metal ion-exchanged zeolites. 11,13,[17][18][19][20]22,24,27,28,[30][31][32][33][34][35][36][37][38][39][40][41][42][43][44][45] For Cu-exchanged zeolites, a [Cu-O-Cu] 2+ motif has been put forth as a possible candidate for the active site. This complex was identified by a UV-vis spectroscopic signature at around 22 700 cm −1 , and the appearance of this signature correlated with the methane-to-methanol conversion. 11,23 The [Cu-O-Cu] 2+ motif has been extensively characterised, both experimentally and theoretically. 23,26,46 The need for a [Cu-O-Cu] 2+ site for methane-to-methanol conversion is, however, not conclusive. For instance, recently Narsimhan et al. 28 were able to experimentally demonstrate a continuous catalytic methane-to-methanol cycle in Cu-exchanged ZSM-5 without the appearance of the UV-vis signature at around 22 700 cm −1 . Furthermore, in recent theoretical studies 18,47 trinuclear Cu clusters were suggested as active sites for Cumordenite and Cu-ZSM-5, while Kulkarni et al. 48 suggest a Cu monomer as an active site in addition to the [Cu-O-Cu] 2+ and Cu 3 O 3 motifs. Although there exists a controversy as to the exact nature of the active site(s), the use of zeolites seems to be important to prohibit the complete oxidation of methane and steer the selectivity towards methanol. Zeolites are sometimes viewed as inorganic analogues to the naturally occurring enzyme methane monooxygenase (MMO), 49,50 which can convert methane to methanol under ambient conditions and is reported to have a porous structure and Fe or Cu double metal ion active sites, similar to zeolites. 51 Zeolites are formed by corner-sharing SiO 4 and AlO 4 tetrahedra or T-sites. Such sites can be joined in many different ways, creating a multitude of structures, including e.g. mordenite, chabasite, and MFI. 52,53 ZSM-5 is an example of the latter. By replacing Si 4+ with Al 3+ , different cations, such as H, Cu, Fe, or Co, can be introduced to compensate for the change in charge in the zeolite structure. For any given Si : Al ratio, there exists a distribution of Al T-sites. 54 Depending on the number of Al T-sites per cation, their oxidation states can vary. The common oxidation states for Cu are 1+ and 2+, 55 indicating that one cation can in principle compensate for the charge of one or two Al T-sites. 56 The distribution of Al and the number of Al pairs can be controlled to some degree by tuning the synthesis method. 57,58 However, studies have shown that for a Cu dimer site, and especially the [Cu-O-Cu] 2+ motif, there is a clear energetically favourable configuration within the ZSM-5 structure at the crossing between the straight and sinusoidal channels. 46 An equally clear candidate for the active site in other metal-exchanged zeolites, e.g. Fe-, Ni-and Co-ZSM-5, has yet to be identified. One remaining question to address is whether the [Cu-O-Cu] 2+ motif is also an active site in other metal-exchanged zeolites. This question is well suited for first-principles studies and can be addressed by exchanging Cu with known active elements (Ni, Co, Fe), and further scrutinised by also including Ag and Au, which have the same electronic configuration as Cu. Henceforth, we will denote the motif as [M-O-M] 2+ , which is the structure obtained when replacing Cu with one of the other metals. Furthermore, the level of description of active sites in the ZSM-5 structure needs to be investigated. Small cluster models have been used, 23,46 and the question remains whether such small models adequately capture the effect of a pore structure.
Here, we investigate the [M-O-M] 2+ motif in the ZSM-5 structure and its implications on the reaction kinetics of methane-to-methanol conversion by combining firstprinciples methods and micro-kinetic modelling. To determine the model size needed to accurately describe the reaction on the dimer motif, we study cluster models of different sizes for the Cu-ZSM-5 structure, and compare these to the periodic crystal structure. The focus is on the kinetics of the removal of the oxygen linking the two Cu atoms in the [Cu-O-Cu] 2+ motif, which is the α-oxygen. The mechanism for adding the α-oxygen is not studied here as it is not fully understood. One possibility is the formation of a peroxo intermediate 25 with subsequent removal of the first and second oxygen atoms. The removal of the first oxygen has been proposed 20,22 to be assisted by a Cu + spectator site. If this is the case, our model describes the removal of the second oxygen, which is expected to be rate-limiting. The influence of water is not scrutinised, as experiments are typically carried out without water in the reaction feed. 11,21,59 After finding a suitable cluster model system, we address the question of how generic the [M-O-M] 2+ motif is, via replacing the Cu by Ni, Co, Fe, Ag or Au, and investigate the transient and steadystate kinetics of methane-to-methanol conversion.

Method
The first-principles calculations were performed using density functional theory (DFT), as implemented in the DMol3 software. [60][61][62] The one-electron Kohn-Sham orbitals were expanded in an all-electron double-numeric basis set together with polarization functions (dnp), and the cut-off radius was set to 5 Å. The density was converged until the difference was less than 10 −6 a.u. between two subsequent self-consistent steps. The Coulomb potential was obtained by projection of the charge density onto angular dependent weighting functions centered at each atom. The Kohn-Sham equations were solved self-consistently using an integration technique based on weighted overlapping spheres centered at each atom. In the fully periodic model, only the gamma point was considered in the reciprocal space integration over the Brillouin zone. The exchange-correlation was approximated using the Perdew-Burke-Ernzerhof (PBE) spin-polarised formula. 63 A sensitivity analysis was performed by separately changing the energy of each intermediate and transition state along the reaction pathway by ±0.25 eV and calculating the response in the kinetics. The analysis shows that the conclusions are not affected, and that the PBE level of accuracy is sufficient.
Transition states (TS) were found using the climbing image nudge elastic band (NEB) method, 64 and the ionic positions were optimised until the maximum force was less than 0.05 eV Å −1 for the intermediates and the transition states. For the Cu system, the quadratic synchronous transit (QST) method together with a linear synchronous transit (LST) algorithm as implemented in the DMol3 software was also used to find the transition state for methane dissociation.
Entropic contributions to the rate constants of the reaction steps were calculated from vibrational frequencies obtained by displacing the atoms in the [M-O-M] 2+ site and the C and H atoms belonging to the methane molecule, i.e. in total, eight atoms were displaced. Central differences were used to obtain the Hessian matrix with a displacement of 0.01 Å. Low frequencies were truncated at 100 cm −1 for numerical stability.

Structure modelling
Several structure models, based on the experimental structure, 65,66 were evaluated in order to accurately describe the [Cu-O-Cu] 2+ motif. For all cluster models, the terminal Si atoms were saturated with H atoms in order to remove the dangling bonds and mimic the ZSM-5 zeolite framework. The Si-H distance was set to 1.491 Å, which is the relaxed distance of the SiH 4 molecule, in the direction of the replaced Si-O bond.
The smallest model, consisting of 20 T-sites (small model), is very close to the one that was earlier suggested by Tsai et al. 46 The Al atoms were placed in the crossing between the straight and the sinusoidal pores in the ZSM-5 framework and separated by two T-sites from each other (see the ESI †). This configuration was found to be the most energetically favourable. 46 A second model (expanded model) was created by including two T-sites at the bottom of the pore channel, i.e. beneath the active site (in total, 22 T-sites). Finally, T-sites belonging to the top of the pore were included to form a ring model (in total, 40 T-sites). The periodic unit cell (periodic model) with 96 T-sites was also studied with the Al atoms in the same configuration.
The two Cu atoms were placed on top of the Al T-sites, similarly to Tsai et al. 46 For the ring and the periodic models, the favourable configuration was found when the two Cu atoms without any oxygen were in adjacent hollow sites on the channel wall. For the small and expanded models, the empty Cu-Cu site was obtained from the [Cu-O-Cu] 2+ site by removing the α-O atom located between the Cu atoms and relaxing the system. For further details, see the ESI. † For Ni, Co, Fe, Ag and Au, the same structures as for the Cu-ZSM-5 ring model were used but with the metal atoms exchanged. With the exception of the H atoms in the Si-H bonds, all atoms were allowed to relax during the geometry optimisation.
Relevant zero-point corrected energies were calculated with the [M-O-M] 2+ site plus methane in the gas phase as reference, i.e.
where A-Z represents reaction intermediates on the active site in the zeolite, A = CH 4 + O, CH 3 + OH or CH 3 OH, either in the adsorption or transition state geometry, whereas O-Z and CH 4 represent the separate oxygen-activated ion-exchanged ZSM-5 cluster (the [M-O-M] 2+ site) and methane in the gas phase, respectively.

Micro-kinetic modelling
The kinetics of the methane-to-methanol conversion was investigated using a micro-kinetic model with input parameters solely based on the DFT results. The rate constants of the desorption and activated conversion of adsorbed reactants were expressed using the conventional transition state theory (TST), 67 where Z ‡ and Z are the partition functions for the transition state and the initial state, respectively, and ΔE is the zeropoint corrected energy difference. For the non-activated adsorption of methanol, the rate constant is expressed according to the conventional collision theory (TST yields in this case the same expression with a sticking coefficient of one and a gas-phase-like transition state), where m is the molecular mass of methanol, λ = h/(2πmk B T) 1/2 is the thermal wavelength, and A S is the site area (this area is compatible with the cross sectional area of the pore ≃10 Å 2 ). For the adsorption-desorption equilibrium constant, we have where ΔG is the corresponding free-energy difference, Z a and Z g are the partition functions for the adsorbed and the gasphase states, and ΔE d is the zero-point corrected energy difference between the gas-phase and the adsorbed states. From eqn (3) and (4), or from (2) directly, we can obtain an expression for the rate constant of methanol desorption, where the numerator contains the rotational and vibrational partition functions for the transition state, i.e. the bare zeolite ( for the vacant site) and gas-phase methanol (Z rot and Z vib , which is considered as an ideal gas 68 ), and the denominator contains Z vib a which is the partition function for the adsorbed state (including the metal ions).
Starting from the metal dimer site with an oxygen atom, the [M-O-M] 2+ motif, the most direct mechanism for methane-to-methanol conversion is investigated. Specifically, we consider dissociative adsorption of methane with the formation of adsorbed -OH and -CH 3 groups, subsequent recombination of these groups resulting in the formation of adsorbed methanol (*CH 3 OH), and a methanol desorption step. To calculate the time scale characterising the corresponding transient kinetics of this reaction, we first neglect the back-adsorption of methanol from the gas phase. In this case, the kinetic equations describing the reaction steps are: where p CH 4 is the methane pressure, k 1 ± , k 2 ± and k 3 ± are the rate constants corresponding to the first step (methane dissociation), second step (methyl-hydroxyl recombination), and third step (methanol desorption), in the forward direction towards methanol (+) and the backward direction (−), θ O is the fraction of sites with oxygen, θ OH,CH 3 is the fraction of sites containing dissociated methane species with -OH in between and -CH 3 adsorbed onto one of the metal ions, and θ CH 3 OH is the fraction of sites covered with the adsorbed methanol species. These linear equations can be integrated analytically.
The corresponding expressions are, however, cumbersome and we have used numerical integration. Eqn (6)- (8) correspond to the situation when methanol is removed from the reactor containing the zeolite catalyst and methanol adsorption is negligible. To clarify the role of the latter process, we have also calculated the steady-state kinetics including methanol adsorption by modifying eqn (8) to include a methanol adsorption term as (9) where p CH 3 OH is the methanol pressure and the fraction of empty sites is θ

Results and discussion
The gas-phase C-H distances in methane were calculated to be 1.10 Å which compares well with the experimental value of 1.09 Å. 69 For methanol, the C-H, C-O, and O-H distances were calculated to be 1.10 Å, 1.43 Å, and 0.97 Å, respectively. These values compare well with the experimental distances of 1.10 Å, 1.43 Å and 0.96 Å, respectively. 70 The formation energy of methanol from methane and oxygen was calculated to be −1.08 eV, which compares reasonably well with the experimental value of −1.28 eV. 71

Structure model
The size of the cluster model was converged with respect to the oxygen and methanol adsorption energies and compared to the values for the fully periodic Cu-ZSM-5 zeolite. The oxygen adsorption energies are calculated to be −2.13 eV, −1.73 eV and −1.59 eV with increasing cluster size (see Fig. 1). Compared to the adsorption energy of the periodic model which is −1.62 eV, it is clear that the smaller and expanded models overestimate the adsorption energy of oxygen. The methanol adsorption energy does not have the same monotonic decrease. Instead the calculated values are −1.52 eV, −1.18 eV and −1.50 eV for the small, expanded and ring models, respectively, which should be compared to −1.63 eV for the periodic model. Although absolute adsorption values are important, the relative energy difference between reaction intermediates affects the kinetics even more. Here, we see that the difference between the oxygen and methanol adsorption energies changes from 0.61 eV to 0.55 eV, 0.09 eV, and 0.01 eV for the small, expanded, ring, and periodic models, respectively. The difference between the latter two values is small, and we judge the ring model to be a good compromise between accuracy in energies (left vertical axis) and computational cost (right vertical axis). It is clear that the smallest cluster size, which is closest to the one previously used, 46 is not sufficient. It is also clear (see Table 1) that the ring model agrees well with the periodic structure concerning such properties as Cu-O-Cu angle, Cu-Cu distance and Mulliken charge of the α-oxygen. The results obtained for the ring and periodic models differ by less than 0.3%, while for the small and expanded models the difference to the periodic model is around 2% with respect to the Cu-O-Cu angle and Cu-Cu distance, and 7% with respect to the charge of the α-oxygen. The conclusion from this is clear: the small and expanded models are not sufficient to describe the [Cu-O-Cu] 2+ motif.

Energy landscapes
Using the ring model for the dimer site, the oxygen adsorption energies for Cu, Ni, Co, Fe, Ag and Au are calculated considering one half of an O 2 or N 2 O molecule as a reference (see Table 2). The formation of α-oxygen in the Ag and Au zeolites is calculated to be energetically unfavourable, while the opposite holds for Cu, Ni, Co and Fe. Based on these findings, we do not consider the Au-and Ag-dimer sites in the microkinetic model, as oxygen activation of these sites is unrealistic.
Using the ring model, the energy landscapes were calculated with the considered mechanism and different metal ions (see Fig. 2). It is important to note that these energy landscapes start from an already oxygen-activated site.
We consider the geometrically shortest reaction route for methane dissociation (TS1), through what resembles a surface-stabilised transition state with a CH 3 -group bound to one of the copper ions. There exists another path through a CH 3 radical-like transition state, similar e.g. to that found by Latimer et al. 72 for a copper-dimer site in mordenite. Using the climbing image NEB method, we did not find the radical barrier. However, with the alternative LSTQST method we were able to find the radical-like TS for Cu. This barrier, which was confirmed by vibrational analysis, is 0.39 eV lower than the surface mediated barrier. Although the radical-like barrier is lower than the surface-stabilised barrier for Cu, the sensitivity analysis shows that it is only for the [Cu-O-Cu] 2+ motif that the first TS is a rate-controlling step. For Ni, Co, and Fe, the results are not affected by changing the methane dissociation barrier.
The second barrier (TS2) is the recombination of the methyl group (bound to one of the metal ions) and the hydroxyl group (in between the two metal ions) into adsorbed methanol. This adsorbed methanol species can then desorb in the last step of the reaction mechanism.
We see that Cu has smaller barriers than the other metals: 0.85 eV for TS1 (0.45 eV for the radical-like TS1) and 0.78 eV for TS2, compared to 0.61 and 1.98 eV for Ni, 0.83 and 2.36 eV for Co, and 1.37 and 2.30 eV for Fe. For Cu, the energies of the two different TS1 found for the methane dissociation, in particular the surface-stabilised TS, agree reasonably well with previously reported results, e.g. 0.7 eV (ref. 23) and 0.96 eV (ref. 18). The reaction is also the least endothermic for Cu with

Micro-kinetic modelling
Starting from oxygen-covered [M-O-M] 2+ sites and using realistic experimental conditions such as a methane partial pressure of 1 mbar and a temperature of 523 K, we have calculated the transient reaction kinetics by employing eqn (6)- (8).
In particular, Fig. 3a) shows the fraction of empty sites as a function of time. Owing to the simplicity of our reaction mechanism, this fraction is proportional to the amount of desorbed methanol molecules. The time scale of the reaction can in this case be quantified using the time τ 1/2 , which is the time needed to convert half of the α-oxygen and methane into gas-phase methanol. This is shown for Cu, Ni, Co and Fe as a function of temperature in Fig. 3b). The time scale is reasonable for Cu, but too long for Ni, Co and Fe. There exists, however, experimental evidence for methane-to-methanol conversion on Ni-, Co-and Fe-ZSM-5 samples. [12][13][14] This means that for Ni, Co and Fe, either the proposed reaction mechanism is not correct, or the dimer motif is not the active site.
The steady-state reaction kinetics have been studied using eqn (6), (7) and (9) for a methanol pressure of 10 −6 bar. The corresponding steady-state fractions of empty sites are shown in Fig. 4 as a function of temperature. The [Cu-O-Cu] 2+ site converts methane to methanol at reasonable temperatures (around 500 K), while Ni, Co and Fe require unrealistically high temperatures. Note that the zeolite framework breaks at around 1400 K. 73 We see that the two different barriers for Cu influence the transient kinetics in Fig. 3, but do not affect the steady-state kinetics in Fig. 4. This confirms our conclusion from the transient kinetics that the [M-O-M] 2+ site is not a plausible candidate for the active site(s) for Ni, Co and Fe.

Sensitivity analysis
By changing each energy separately in the energy landscape by ±0.25 eV for Cu, Ni, Co and Fe, we have performed a sensitivity analysis of the micro-kinetic model. The outcome demonstrates that our results are robust, and that a different energy description (e.g. including van der Waals interactions) will not change the conclusion that the dimer site is reasonable for Cu, but not for Ni, Co or Fe. In this context, note that our energy change, ±0.25 eV, is of the same order of magnitude as the shift in adsorption energies attributed by Göltl et al. 74 to van der Waals interactions in zeolites. For further details on the results of the sensitivity analysis, see the ESI. † Furthermore, the sensitivity analysis shows which energies are important for the overall reaction in general and the reaction half-time τ 1/2 in particular. It demonstrates that the kinetics for all the elements (Cu, Ni, Co and Fe) are dependent on the energy of the oxygen-activated [M-O-M] 2+ site. The kinetics for Cu is the most sensitive with respect to the TS1 energy, while for Ni and Co, the TS2 energy is more important. For Fe, the stability of the bare ion site is the most important, as the barriers are small compared to the overall energy cost for producing methanol from methane. We note that for Cu, the radical pathway to TS1 affects τ 1/2 , but that the order (Cu, Ni, Co and Fe) will not be affected by either barrier. This analysis confirms that Cu has the shortest reaction half-time,  followed by Ni and Co, and that the half-time for Fe is unrealistic.

Conclusions
We showed that a large cluster model size is needed in order to accurately describe the [Cu-O-Cu] 2+ motif in the ZSM-5 structure. In order to investigate if the widely studied [Cu-O-Cu] 2+ motif is a relevant candidate for the active site in other ion-exchanged ZSM-5 systems, the energy landscape for methane-to-methanol conversion was calculated for Cu, Ni, Co, Fe, Ag and Au. The oxygen activation of Ag and Au was found to be energetically unfavorable. Thereby the dimer motif is not likely for these elements. The transient and steadystate behaviour of the proposed micro-kinetic model for the Cu, Ni, Co and Fe sites was analysed. The results show that for Cu, the dimer site is a plausible candidate for the active site, with a reasonable reaction time at realistic temperatures. Since no activity under reasonable conditions was found for Ni-, Co-and Fe-dimer sites, which experimentally are known to be active for methane-to-methanol conversion, the dimer site is judged to be not the correct candidate for the active sites in these systems. The strong oxygen binding of, for instance, Fe, indicates that a more highly oxidised Fe dimer motif is a more plausible candidate for the active site.
The sensitivity analysis confirms that the results are robust and relatively insensitive to errors in the energies, for instance with respect to the surface-stabilised or radical-like transition state for methane dissociation. It also shows that the reaction kinetics are the most sensitive to errors in the methane-dissociation barrier for Cu, while the sensitivity is higher with respect to the second barrier in the cases of Ni and Co. For Fe, the overall energy difference between methane and methanol is what determines the reaction time.