Selective conversion of methane to methanol has remained a major scientific and practical challenge for several decades and has attracted significant attention recently due to the increased availability of natural gas. The major challenge for this reaction is to only activate the first C-H bond in methane and prevent its over-oxidation to thermodynamically more stable CO/CO2. Therefore, an indirect industrial process has been developed, where methane is converted to syngas 1 , which is then transformed to longer chain hydrocarbons or methanol 2,3 . The methane steam reforming conversion takes place at high temperatures and pressures and is economically not viable at remote extraction sites, where such technology could not become available. To achieve on-site conversion of methane to more valuable and easier to transport methanol, it would be desirable to find a direct route for this reaction, operated at mild conditions. In nature, this is achieved by the enzyme methane monooxygenase, where Cu-and Fe-oxo clusters serve as the active site [4][5][6] . Inspired by this enzymatic catalysis, extensive research has been focused on reproducing similar active sites in heterogeneous catalysts [7][8][9] . In particular, Cu-exchanged zeolites have drawn significant attention in this context [10][11][12][13][14][15][16][17][18][19][20][21][22][23] .
Due to strict conversion/selectivity limits in the catalytic conversion of methane to methanol 24 , a stepwise conversion process shows greater promise in comparison with the continuous-flow process 25 . In the stepwise process, the catalyst is activated in an oxidizing atmosphere, followed by exposure to methane. Contact of the activated catalyst with methane triggers the formation of methoxy species, which are subsequently extracted as methanol by using water vapor.
After the initial discovery of the activity of transition-metal (TM) exchanged zeolites, several oxidants 11,[26][27][28] and zeolite structures have been tried for this reaction 17,22,26,29 .
However, despite significant efforts, the nature of the active site is still under debate. This is best illustrated for Cu-exchanged Mordenite, where different studies report the presence of either Cu-oxo dimers 10,16,26,[30][31][32] or trimers 18 as active sites for this chemistry. Further, several studies for different zeolite structures indicate a change in the active sites as a function of experimental conditions 20,30,33,34 . At the same time, full experimental characterization of a complex system like a zeolite is challenging. These difficulties are rooted in the presence of a distribution of active sites, which show overlapping signals in Raman spectroscopy and UV-vis spectroscopy 29,35 . Even combining information with EXAFS data on coordination environments and bond distances requires a series of assumptions about the active site 33,34 .
First-principles based modeling can support and extend these experimental assignments by systematically considering an extensive set of possible active site structures and by calculating their relative thermodynamic stability 36,37 . This information can then be summarized in phase diagrams, which reveal the most stable phase under various experimental conditions [38][39][40] . Combining the information from phase diagrams with knowledge of the anchoring point distribution, enables derivation of the site distribution present in the zeolite material under specific experimental conditions 37,41 . In the conversion of methane to methanol using Cu-exchanged zeolites, attempts so far have aimed at describing the relative stability of different Cu sites in Mordenite 42 , ZSM-5 19 or SSZ-13 29,43 . These studies, however, considered only few possible active sites, which does not describe the complex distribution of sites expected to be present in such a zeolite system.
The ideal candidate zeolite-framework to study the site-speciation of Cu is SSZ-13 (Cu- SSZ-13). This is a zeolite in the Chabazite framework with a highly symmetric primitive unit cell, which limits the number of possible Cu configurations 37 . This material has been shown to efficiently convert methane to methanol 11,17,21,29,33,34 and a significant body of work on characterization of the active sites in this material already exists 29,33,34 . Lastly, Cu-SSZ-13 has been extensively studied in the context of deNOx-SCR 44,45 . In fact, it is possible to extend the methodology developed to describe active sites during deNOx-SCR 39 in order to gain a comprehensive understanding of Cu-SSZ-13 in the methane to methanol conversion.
In this work, we use first-principles modeling to study the nature of the active sites in Cu-SSZ-13 during the stepwise conversion of methane to methanol. We develop a detailed theoretical model based on the chemical potential of Cu in various monomer, dimer and trimer structures in the zeolite matrix and use phase diagrams to determine the thermodynamically preferred state of Cu. We find that, depending on the local Al distribution, either Cu2O2H2 or Cu2OH, is stabilized during catalyst activation. When methane is introduced, methanol adsorbed to Cu monomers is formed, which desorbs when the water pressure in the system is increased. Our findings are capable of rationalizing experimental EXAFS measurements reported in the literature 33 .
One of the major challenges in converting methane to methanol is to prevent methane overoxidation to thermodynamically preferred CO2. In catalytic conversion, selectivity/conversion limits have been derived for this chemical reaction and several strategies have been suggested to circumvent these problems 24 . For transition metal exchanged zeolites, in particular, the stepwise conversion of methane to methanol using O2 as the oxidant has drawn significant attention 25 . The catalyst is first activated in an O2 containing atmosphere at high temperature to generate the active sites. In a subsequent step, the system is cooled down and exposed to methane. Interaction with methane at this process step triggers the formation of surface methoxy species (i.e., oxygenates of CH4 such as methanol), which cannot desorb. Further oxidation of methanol would require its desorption from the Cu site, followed by regeneration of the active sites. The regeneration of the active sites would require oxidative conditions, but since oxygen is absent during methane exposure, no further active sites are formed, and stably adsorbed surface methoxy species are protected from overoxidation. In a final step, methanol is extracted at lower temperature by introducing water vapor to the system.
For the studies reported in the literature, the fundamental steps remain similar, i.e., activation is followed by methane exposure and methanol extraction. However, the exact experimental conditions can vary between different research groups, which can influence methanol yields 22 . For Cu-SSZ-13, a systematic analysis of the impact of various experimental parameters on the performance of each step in this process has been reported in the literature 33 . Here, we will focus on conditions described as optimal for Cu-SSZ-13 Typically, at points 1, 3 and 5, conditions change abruptly (see Fig. 1), while the system is kept at specific conditions for several hours for points 2, 4, and 6 (see Fig. 1). At some points of the process (e.g., after 3, 4, and 5), the system is exposed to inert gases (such as He) to remove excess gas pressures from previous exposure. Inert gases do not influence site speciation and are, therefore, omitted from further discussion below.
While the externally applied gas pressures are known, residual pressures of the other gases, caused by impurities in the supplied gases or small amounts of gases adsorbed in the zeolite nanopores from previous exposure, will also be present. The exact values will depend on the experimental setup and the detailed process protocol. However, since residual gas pressures are impossible to measure experimentally, we use approximate values for them.
If not explicitly stated otherwise, we assume a logarithmic O2 gas pressure ln(P O 2 /P 0 )=-7, a logarithmic H2O gas phase pressure ln(P H 2 𝑂𝑂 /P 0 )=-7 and after exposure to CH4, a residual logarithmic CH4 pressure ln(P CH 4 /P 0 )=-9. All pressures are calculated with respect to a reference pressure P0 of 1 bar. The list of all parameters at each point along the stepwise conversion process for methane to methanol is given in Supporting Information, Table S1.
Careful testing revealed that the choice of residual pressures has only a subtle effect on results. Wherever variation of the residual pressures affects the results, we explicitly discuss them.
The Thermodynamic Model:
Structures included in the model: We focus on the Cu exchanged zeolite SSZ-13. SSZ-13 is a material with the chabazite structure, the zeolite framework structure with the smallest primitive unit cell. This feature reduces the complexity of the system, while still retaining all basic features associated with zeolites, which makes this material a nearly ideal testsystem for zeolite catalysis. 48 , all of which are shown in Fig. 2b. To study Cu speciation within the SSZ-13 zeolite framework, we consider an extensive set of plausible Cu sites. We start our discussion with mononuclear Cu (Cu1) sites included in our model, which have been extensively discussed in the literature 36,38,39,[49][50][51] , followed by the Cu-dimer and Cu-trimer sites.
In Cu1 species can be stabilized in the 6R or the 8R (exchange sites are shown in Fig. 2 (c)).
During the stepwise conversion of methane to methanol, the zeolite is exposed to a series of different gas-phase environments. It has been reported that, in particular, the presence of H2O changes Cu coordination and leads to Cu-H2O complexes at lower temperatures 38,39 . We, therefore, explore all the different possibilities to form Cu1(I) and Cu1(II) located in a 6R and 8R, for all Al configurations shown in Fig. 2 (b), and allow for the coordination of Cu1 with up to six H2O molecules, which will allow for the formation of Cu1-hexaaqua species. A detailed discussion of all Cu1(II) has been given in the literature and we rely on the related published structures 39 . For Cu1(I), we follow a similar strategy, and all the optimized structure files are provided in the Supporting Information. Only the most stable Cu1(I) and Cu1(II) structures for a given number of adsorbed H2O molecules are included in our thermodynamic model. In subsequent phase diagrams, these structures are labeled as: Cu(I)+nH2O and Cu(II)+nH2O for Cu1(I) and Cu1(II) sites, respectively, with n adsorbed H2O molecules each.
In the literature, several different Cu dimers and trimers have been suggested to be active in the conversion of methane to methanol 18,26,29,42 . Similar to the Cu1 case, Cu in dimers and trimers will bind to O atoms adjacent to framework Al, but these structures are generally too large to fit into the double 6R structure of a unit cell. We, therefore, constructed four different local Al configurations reaching over adjacent unit cells, two of which are located in the same 8R (see D-A and D-B in Fig. 3 (a) and 8R in Fig. 3 (b)) 29,35 .
The other two Al configurations allow for Cu dimers bridging a 6R and an adjacent 8R (D-C and D-D in Fig. 3 (a) and 6R/8R Fig. 3 (b)). We construct Cu2OyHz (y=1, 2; z≤y) dimers, for all four Al configurations (D-A though D-D). For the Cu2O2H2 stoichiometry we furthermore considered the formation of associated monomers, which have been suggested
as active sites for methane to methanol conversion in Mordenite 52 and Zeolite Omega 53 .
Additionally, the presence of Cu trimers bound in 8R structures has been suggested for Mordenite 18 . We, therefore, include Cu3O3Hz (z≤3) clusters for the D-A and D-B Al configurations. Accordingly, and herein, the different clusters are denoted by X-CuxOyHz, where X stands for A through D, which denotes Al configurations D-A through D-D, while
x, y and z describe the cluster stoichiometry. All structures were optimized using the PBE-TS 54,55 density functional and periodic boundary conditions as implemented in VASP 56,57 . Several spin states were probed for each configuration. The spin ground state is used for further analysis and relative energies for the different spin states are given in Supporting Information, section S2. All optimized structures are shown in Fig. 4; structural files are provided in Supporting Information.
In our thermodynamic model, we assume that -due to the long activation and reaction times -the system reaches its thermodynamic equilibrium during each step. We, therefore, rely on the thermodynamically most stable of all considered configurations and neglect potential intermediate states of the system. Whenever necessary, we will explicitly discuss cases where intermediate states might be important to understand the performance of Cu-SSZ-13 in the stepwise conversion of methane to methanol.
To develop a thermodynamic model, we rely on structures obtained from static structural optimization using PBE-TS 54,55 . While functionals using the generalized gradient approximation are known to lead to reasonable structural guesses (i.e. systematic errors of ~1-2% in terms of bond length 58 ), comparison with post Hartree Fock/post DFT methods such as MP2 59 or the Adiabatic-Connection Fluctuation-Dissipation Theorem in its Random Phase Approximation (RPA) 60 , shows significant errors for reaction energies in the deprotonation of isobutene in ZSM-5 and in the conversion of methane to methanol over Fe-oxo sites in SSZ-13 61,62 . Furthermore, RPA has shown excellent performance in describing the adsorption of alkanes in protonated zeolites 63,64 . We, therefore, use RPA to calculate total energies for each PBE-TS optimized configuration. To derive the finite temperature Gibbs Free Energies reported here, we incorporate vibrational zero point energies, and entropic corrections from static vibrational corrections calculated using the PBE-TS functional 39 .
Since a Cu-SSZ-13 particle is not connected to an external Cu reservoir, we assume that after ion exchange, the number of Cu atoms in the zeolite matrix is constant. At any point in time, a distribution of Cu sites is present in the material and the total energy of the system is the sum over the energies of all Cu atoms and every time a Cu atom changes its coordination environment and moves from one type of active site to another, the total energy of the system changes. Additionally, only a modest increase in methanol yield with activation time 33 indicates that Cu cations move quickly in the zeolite matrix compared to the time-scales typically associated with the various steps in the experimental procedure.
Therefore, Cu will always be present as the configuration in its thermodynamically preferred position with the lowest chemical potential μ Cu (i.e., the lowest energy per Cu atom). Since we are interested in a reaction environment containing O2 and H2O, we calculate μ Cu for sites anchored in unit cells containing 2Al as
where, G indicates the Gibbs free energy of the structure in the superscript, where CuxOyHz-zeo indicates CuxOyHz bound to the zeolite framework and 2H-zeo refers to a zeolite framework passivated by Brønsted protons. The gas phase chemical potentials of O2 (𝜇𝜇 𝑂𝑂 2 ) and H2O (𝜇𝜇 𝐻𝐻 2 𝑂𝑂 ) are obtained by correcting the energies obtained from static electronic structure PBE-TS calculations using vibrational zero-point corrections, translational, vibrational and rotational entropies, and pressures. When only one Al atom is present in the unit cell (for the 1Al Cu1 case), the formula for μ Cu is
Throughout this work, we will use μ Cu to calculate phase diagrams. In a realistic zeolite system, Cu will occupy the most stable Al configuration and as soon as it is filled, the next most stable Al configuration will be occupied 36,39,65 , which we account for in our study.
The exact nature of the observed Cu species will therefore depend on the distribution of local Al configurations, which is expected to vary based on the synthesis and material parameters across the zeolite samples 66 .
Here, we systematically study the impact of the Al distribution on the observed Cu sites.
More specifically, to identify the thermodynamically preferred states for Cu in the zeolite framework, we focus on the influence of two parameters on the phase diagram: (i) the
All calculations were performed using the Vienna Ab-Initio Simulation Package (VASP) 56,57 , a plane wave code using PAW pseudo potentials 67 , adapted by Joubert and Kresse 68 . All calculations were performed with an energy cut-off for plane waves of 420 eV and were restricted to the Γ-point. The basic unit cell parameters for periodic calculations are given in the literature 69 and the zeolite volume was set to 830 Å 3 . As described in the main text, all structures were optimized using Density Functional Theory in the parameterization of Perdew, Burke and Ernzerhof 54 . van der Waals interactions were introduced using the Tkatchenko-Scheffler force field 55 . Prior to optimization a 5 ps molecular dynamics simulation at 500 K using the Andersen thermostat 70 was performed and three different structures obtained after equal simulation times from the last 1 ps of this simulation were optimized. In all calculations the spin-states of the different clusters remained fixed and the spin-state leading to the lowest energy was used for further analysis.
All minimum spin states are given in Supporting Information, section S2. For these spin states RPA calculations were performed 60,62 . Here the energy cut-off was increased to 600 eV and the energy cut-off for the response function was set to 250 eV. RPA calculations were restricted to the spin ground state determined in DFT calculations. In a subsequent step, harmonic vibrational frequencies were calculated at the PBE-TS level by numerical differentiation of the forces using a second-order finite difference approach with a step size of 0.01 Å. Vibrational and translational entropies as well as zero point vibrational corrections for gas phase molecules and zeolite unit cells were calculated using the code thermo.pl 71 , a code provided by the National Institute of Standards and Technology. The impact of including translational entropies for the zeolite unit cells on 𝜇𝜇 𝐶𝐶𝐶𝐶 , the chemical potential of Cu, is discussed in the supporting Information, section S3. To remove unphysical translational modes and low energy vibrational modes, only vibrational frequencies above 50 cm -1 were considered in the analysis. For computational efficiency we reproduced dimer configurations D-A through D-D, which reach over two adjacent double six O-ring structures, in a single unit cell using periodic boundary conditions. To confirm the validity of this approach, we compared energies for dimers in the Al configuration D-C for a single and double primitive unit cell and found that energies per
Cu atom calculated at the PBE-TS level lie within 6 kJ/mol (see Supporting Information, section S3). Such small energy differences indicate that constructing dimer and trimer structures in a single primitive unit cell is a good approximation. Gas phase molecules were modeled in a 10 Å x 10.1 Å x 12 Å box.
Phase diagrams:
The first step in the stepwise conversion of methane to methanol is the activation of the catalyst through its exposure to the oxidant for generating the active sites. Here two variables play a crucial role, namely the O2 pressure (P O 2 ) and temperature (T). The T/P O 2 phase diagrams are presented in Fig. 6 Points marked with numerals: 1, 2, 3, and 4, correspond to specific points in the stepwise methane to methanol conversion, as defined in Fig. 1. A detailed legend is given in Fig. 5.
So far, we have studied phase diagrams where only one possible dimer/trimer Al configuration and one monomer configuration was included. However, in a realistic system, multiple Al configurations will be available. Here, initially Al configurations leading to most stable CuxOyHz species will be occupied. While extensive discussion about relative stability of Cu1 in SSZ-13 exists in the literature 36,38,39,48,65 , the relative stability of dimers has been discussed on a limited basis. We, therefore, study phase diagrams including all four Al configurations for dimer formation (D-A through D-D) and Cu1 bonded to the 1Al configuration (see Fig. S3) and find that A-Cu2OH and A-Cu2O2H2 are the most stable Cu dimers (Fig. S3) when all dimer exchange sites are available. In a realistic system, only a finite number of Al configurations D-A will exist, therefore, after they are filled with dimers, the next most stable sites will be formed 37 . We remove Al configuration D-A from the phase diagrams and find that B-Cu2O2H2 is the second most stable site at low T/high P O 2 , while D-Cu2OH is most stable at high T/low P O 2 . When removing either D-B or D-D we find that the other site becomes dominant, which indicates that C-Cu2OH is the least stable dimer structure and will therefore be formed last.
Using these phase diagrams, we can follow the sites during the first four steps of the stepwise conversion (see Fig. 1, Fig. 6, and Fig. S2). Initially, the system is at 323 K under O2 pressure of 1 bar (point 1 in Fig. 1). For all 2Al configurations (Fig. (a)-(c)) is preferred over dimer formation. However, we cannot exclude that dimer decomposition is slow, compared to the reaction times, and therefore dimers might initially still be present at these conditions.
It is interesting to see that at an activation temperature of 773 K several local Al configurations lead to the formation of Cu1(II). The activation temperature for zeolites has been a topic of discussion 11,22,33 . To accelerate dimer formation the highest possible temperature that still allows for ideal performance is preferred. However, at some temperature, monomer formation becomes thermodynamically favorable. Typically, an upper limit of 723 K is assumed, but Pappas et al. report an increase in methanol production at 773 K 33 . In our analysis at this temperature, the formation of monomers was preferred for two of the 2Al configurations (2Al-A (Fig. 6 (b) and 2Al-D (Fig. S2 (c)). However, we relied on a residual O2 pressure (P H 2 O ) of e -9 . When studying the P H 2 O dependence during activation for Al configuration D-A and Cu1 in 2Al-A (see Fig. S4), we find that an increase in P H 2 O shifts the transition to Cu1(II) to higher temperatures, i.e., a slightly higher P H 2 O will allow for higher activation temperatures during dimer formation.
The nature of the active sites in Cu-SSZ-13 during catalyst activation has been studied in detail using in situ EXAFS measurements 33 , and the following qualitative observations were made: (i) a minority of Cu sites do not form dimers during the activation process and
(ii) the coordination number of a fraction of Cu atoms increases from three to four when Cu-SSZ-13 is activated in O2 at 773 K and 473 K. As we will show in the following, both observations are well reproduced from the phase diagrams constructed here. Fig. 6 (b) and Fig. S2 (c) show that monoatomic Cu1 in the 2Al-A and 2Al-D configurations is more stable than dimers or trimers at high T and high P O 2 . Therefore, Cu1 in the 2Al-A and 2Al-D configurations will not form dimers during the high temperature activation process, which in agreement with the experimental observations. Furthermore, for Al configuration D-A, threefold coordinated Cu in the dimer A-Cu2OH will be formed at 773 K, but fourfold coordinated Cu in A-Cu2O2H2 will be most stable for activation at 473 K. Cu in the next most stable dimers B-Cu2O2H2 and D-Cu2OH will not change its coordination number.
Therefore, for activation at 773 K, Cu hosted in dimers in Al configuration D-A will be in a three-fold coordinated structure, while after activation at 473 K, Cu in Al configuration D-A will be in a four-fold coordinated structure, which agrees perfectly with experimental observations documenting an increase in Cu coordination number at lower activation temperatures.
EXAFS measurements containing information about the coordination shells of Cu and distinct coordination peaks at distances of 1.86±0.05 Å, 1.97±0.04 Å, 2.72±0.02 Å and 3.41 Å have been observed 33 . To reproduce these values, we studied the structural information of all the dimers and compare them to experimentally measured values (Fig. 7). We split the information into distances between Cu-O in the dimer, Cu-O in the framework, Cu-Cu and Cu-Al atoms. We find that distances between Cu and the O atom for hydroxyl groups in the dimers, correspond (within error) to 1.86 Å, while the While structural agreement between the stable sites and experiment is encouraging, the presence of other sites such as CuOH, Cu2O2, Cu2O, or Cu trimers, which have been previously suggested in the literature to be the active sites in the methane to methanol conversion 10,16,18,29,72 , is still possible. We therefore discuss bond-lengths for these sites in 𝜇𝜇 𝐻𝐻 2 𝑂𝑂 (𝑇𝑇,𝑃𝑃 𝐻𝐻 2 𝑂𝑂 )
, where, definitions provided earlier are extended here by w, the number of CH3 groups in the system, P CH 4 , the CH4 pressure and µ CH 4 , the chemical potential of CH4 in gas phase.
Using this definition, it is possible to perform an analysis similar to the one performed for the zeolite activation phase during exposure to oxidant.
The last step in the stepwise conversion of methane to methanol is the extraction of methanol by introducing water vapor to the system (see Fig. 1). Accordingly, we focus on a T/P H 2 O phase diagram in the presence of methanol (see Fig. 10, Fig. S7 in Supporting Information; a detailed legend for these figures is given in Fig. 5). We find that for Cu1 in the 2Al configurations (see Fig. phase diagrams (see Fig. 9 and Fig. S6), where adsorbed methanol was not stabilized for the corresponding Al configurations.
As discussed earlier, the methodology presented here almost perfectly reproduces experimentally measured bond lengths and coordination numbers reported in the experimental literature 73 The presence of hydroxylated dimers raises significant questions relevant to methane activation. Methane activation over oxo-sites has been extensively discussed and is often proposed to proceed via an abstraction and rebound step 74,75 or a concerted mechanism 76 .
In zeolites only Cu-OH has been discussed as a potential hydroxylated active site 72 .
However, the site geometries presented in this work imply different mechanisms (see Scheme 2) and based on our results we hypothesize the following: For Cu2(OH) stoichiometric arguments point towards H2 formation during the activation of CH4 (Scheme 2 (a)). Indeed, experimental measurements indicate the formation of small quantities of H2 during methane exposure 28,33 . On the other hand, for the Cu2O2H2 sites a similar mechanism activating two CH4 molecules and forming two methanol molecules is possible One of the most surprising observations made in this work is that at no point Cu-OH, a site that is believed to be active in the selective catalytic reduction of NOx 38,45 and has been suggested to be active in the conversion of methane to methanol 72 , is stabilized. It remains to be seen whether the presence of NH3, as encountered in deNOx-SCR, can stabilize Cu-OH over hydroxylated Cu-dimers.
yellow Si, blue Cu, white H and blue-grey Al, respectively. *At 1 and 6, all H2O molecules not directly coordinated to Cu are omitted for clarity.
In particular, the prediction that hydroxylated dimers are thermodynamically preferred and at no point Cu-oxo sites, such as Cu2O, Cu2O2, or Cu3O3, are stable, contradicts common assumptions about the nature of active Cu centers in zeolites and also enzymes active in the conversion of methane to methanol. In the future, it will be interesting to see to what extent our findings can be utilized to arrive at an improved understanding of this reaction and how we may be able to optimize the conditions for the stepwise process of methane to methanol conversion.
Supporting Information: Supporting Information contains information about the studied reaction conditions along the stepwise conversion of methane to methanol, the preferred spin states for Cu dimers and trimers, the effect of translational entropy, the effect of unit cell size, a graphical legend to the phase diagrams and phase diagrams for Al configurations 2Al-B through 2Al-E. Furthermore, phase diagrams for the relative stability of dimers in the different exchange sites is given. Additionally, all dimer/trimer/Cu1(I) structure files are provided.