Para-xylene (PX) is an important raw aromatic feedstock for the industrial production of terephthalic acid and dimethyl terephthalate, which are both monomers for the commercial production of polymers. Currently, PX is mainly produced by catalytic reforming and toluene disproportionation processes ^[1]. Given that worldwide petroleum reserves represent a dwindling resource, the raw materials used for catalytic reforming are gradually decreasing. Furthermore, toluene disproportionation and xylene isomerization reactions generate large quantities of benzene as an unwanted byproduct. The purification of the PX generated using these routes is rather costly, because of the similarities in the boiling points of PX (411.5 K), meta-xylene (MX, 412.3 K), ortho-xylene (OX, 417.6 K) and ethylbenzene (409.4 K). The catalytic alkylation of toluene with methanol has been explored extensively during the past four decades as a promising alternative to catalytic reforming and toluene disproportionation for the production of highly pure PX.

Several zeolites have been reported as catalysts for the methylation of toluene ^[2-4]. HZSM-5, in particular, has attracted considerable interest, because its pores allow for the diffusion of PX out of the channels around 100 and 1000 times faster than OX and MX ^[5], respectively. The selectivity of zeolite-mediated methylation processes for PX can be enhanced up to 90% by modifying the catalyst with inorganic agents (e.g., Mg, B, Si, P, and Pt) ^[6-11] or optimizing the process variables ^{[1, 2]}. These modifications are mainly used to increase the relative diffusion rate of PX or reduce the number of acid sites located on the external surfaces that can cause undesirable secondary reactions.

In addition to forming xylenes as the major products during the alkylation of toluene with methanol, several liquid products (e.g., benzene, trimethylbenzene and ethylbenzene) and gaseous products (e.g., methane, ethane and propene) are also formed. To increase the selectivity of this process for PX and minimize the formation of byproducts, we need to develop a deeper understanding of the reaction mechanism. Quantum chemical calculations are indispensable when it comes to elucidating reaction mechanisms at the molecular level. For modeling reactions in zeolites, cluster modeling and our own N-layered integrated molecular orbital and molecular mechanics (ONIOM) methods are generally used. However, these methods are inadequate for describing the steric constraints imposed by micropore walls ^[12]. Increasing computational capabilities have led to an increase in the popularity of the periodic method because it provides more quantitative information on the effects of the micropore walls on the guest molecules, enabling confinement effects to be described with greater accuracy. Fundamentally, toluene methylation also represents an important step in the aromatic cycle of the methanol-to-hydrocarbon (MTH) process. A detailed mechanistic and kinetic understanding of this step could therefore provide basic information for the development of new catalytic processes and the optimization of existing operating conditions.

As for the alkylation of toluene with methanol, two distinct routes have been proposed, namely the concerted and stepwise pathways, based on theoretical and experimental works ^[13]. The primary distinction between these two processes is that the stepwise pathway involves the formation of a covalently bound surface-methoxy group as a reaction intermediate, whereas no such intermediate is involved in the concerted pathway. Several computational reports using cluster modelling and the ONIOM method have shown that the formation of surface-methoxy groups by the dehydration of methanol involves intrinsic energy barriers of more than 200 kJ/mol ^{[14, 15]}. Based on these results, the stepwise pathway was long considered implausible. However, the results of recently reported periodic density functional theory (DFT) calculations suggested that much lower energy barriers (100-150 kJ/mol) were possible for the formation of surface-methoxy groups from methanol at different locations on the mordenite zeolite ^[16]. Furthermore, spectroscopic results from in situ FTIR and NMR experiments have convincingly demonstrated the existence of surface-bound methoxy groups, thereby corroborating the fundamental viability of a stepwise mechanism. However, experimental kinetic measurements are much more consistent with the concerted pathway ^{[12, 17]}. Brogaard et al. ^[18] reported that the stepwise pathway prevailed during a typical MTH reaction because of the higher entropy loss associated with the concerted pathway. Li et al. ^[19] computed the reaction mechanism for the methylation of toluene with methanol over HZSM-5 based on electronic energies and the ONIOM method. Vos et al. ^[20] calculated the concerted mechanism for the mordenite-catalyzed methylation of toluene using a periodic method, and found that the electronic energy barrier for the formation of PX was around 20 kJ/mol lower than those of its isomers. However, to our best knowledge, a systematic investigation of the methylation of toluene with methanol over HZSM-5 has not been reported by using periodic DFT.

In this study, we have evaluated the concerted and stepwise pathways for the methylation of toluene using periodic DFT. We have also compared the calculated results with the experimental values to develop our understanding of the mechanism and shape-selectivity for this reaction.

Self-consistent density functional calculations were performed using the Vienna Ab Initio Simulation Package (VASP) ^{[21, 22]}. The projected augmented wave (PAW) method ^{[23, 24]} was used in combination with the PBE exchange-correlation functional ^[25] and a plane wave basis set kinetic energy cutoff value of 400 eV. Periodic boundary conditions were applied in all of the calculations using a real-space grid spacing of 0.10 Å. In the adsorption calculations, the super cell was identical to the unit cell and Brillouin zone sampling was restricted to the Γ-point.

The lattice constants of silicalite-1, the purely siliceous version of ZSM-5, were determined using the approach reported by Svelle et al. ^[26]. The high-temperature structure of ZSM-5, which has an orthorhombic unit cell with lattice constants a = 20.08 Å, b = 19.89 Å and c = 13.37 Å, was used in this study ^[27]. The parameters of the unit cell were optimized using a force threshold of 0.01 eV/Å and an energy cutoff of 800 eV ^[26]. The optimized parameters a = 20.32 Å, b = 20.16 Å and c = 13.46 Å were used in all of the subsequent calculations.

Transition state structures were located using the improved-dimer method and were considered converged when the forces on all of the atoms were less than 0.03 eV/Å ^[28]. Enthalpy and entropy calculations treated the gas-phase species in the rigid rotor harmonic oscillator approximation and employed the harmonic approximation for adsorbates and transition states. A partial vibrational Hessian approach was used for the frequency calculations based on the atoms of the organic molecule, as well as the H atom of HZSM-5, whereas the atoms of the zeolite framework were fixed. This subset of atoms was shown to be sufficient for calculating the thermodynamic quantities of alkanes in HZSM-22 ^[29]. The numerical differentiation of the forces was achieved using a displacement step of 0.015 Å and an energy convergence threshold of 10^-9 eV/electron. Each of the aforementioned atoms was displaced +δ and -δ from their corresponding equilibrium positions along each of the three Cartesian coordinates. Zero-point-energy (ZPE) corrections were calculated using statistical mechanics based on the Boltzmann distribution, and subsequently used to correct the total energies.

DFT methods are unable to properly describe the van der Waals dispersive contribution, which is critical for understanding reactions in zeolites ^[26]. The total energy of a system (E_DFT-D) can therefore obtained by adding a semi empirical dispersion contribution (E_disp) to the self-consistent Kohn-Sham energy ^[30] according to the formula E_DFT-D = E_DFT + E_disp + ZPE. E_disp can be obtained using the dftd3 program, which is proposed by Grimme et al. ^{[30, 31]}.

The Gibbs free energies of all species at temperature T are given by G(T) = E_DFT-D+ dH(0→T) - TS^T, where dH(0→T) is the change in the enthalpy from 0 to T K, and TS^T is the entropy at temperature T. For surface structures, the translational and rotational modes can be replaced with vibrational modes corresponding to frustrated translation and rotation on the surface, and all of the terms in G(T)can be calculated in the harmonic approximation based on the calculated vibrational frequencies. For gas-phase molecules, the entropy and enthalpy terms in G(T) can be calculated using the Shomate equation, and the related parameters are taken from Ref. ^[32].

The rate constant k obtained using transition state theory is defined as follows: k = (k_BTh^-1)exp(-△G₀^≠/RT), where k_B, h and R are the Boltzmann, Planck and universal gas constants, respectively, and △G₀^≠ represents changes to the standard molar Gibbs free energy.

DFT energy calculations for gas-phase molecules were carried out using cubic boxes with edge lengths of 20 Å to minimize interactions with the periodic images ^[33]. A force threshold of 0.01 eV/Å was used in these calculations.

The first step of this investigation was to build an appropriate model of the HZSM-5 zeolite. To create an acidic site in HZSM-5, one of the 96 Si atoms in the unit cell was replaced with an Al atom at the T12 site ^[34] and the resulting negative charge was compensated by a proton bonded to one of the neighboring framework oxygen atoms. Specifically, we selected the Al12-O(H)-Si(3) site because of its location at the intersection of the straight and sinusoidal channels, which provided the maximum available space and created the most accessible active site to accommodate various guest molecules ^{[33, 34]}. The resulting model for HZSM-5 was similar to other reported models ^{[33, 34]}.

The optimized co-adsorption geometries for toluene and methanol are displayed in Fig. 1(a)-(c). The results revealed that methanol was adsorbed on the acidic site by an end-on hydrogen-bond structure, forming two strong hydrogen bonds between the hydrogen on the protonated methanol and the oxygens of the deprotonated acid site, which was in agreement with other theoretical studies ^{[35, 36]}. Toluene was adsorbed at the intersection of the straight and sinusoidal channels, and various co-adsorption geometries with various orientations for the methyl group in the toluene molecule were obtained, which were denoted as Coads_P, Coads_M and Coads_O. In terms of the geometries of Coads_M and Coads_O, the methyl group projected towards the straight channel, and there was an equivalent distance of 3.74 Å between the carbon atom of the methyl group in methanol and the carbon atoms of the benzene ring. In contrast, the distances between the carbon atom of the methyl group in methanol and the carbon atoms of the benzene ring in Coads_P were different. Some of the other geometric parameters in Coads_M and Coads_O were also equivalent, suggesting that their adsorption energies would be similar. As shown in Fig. 2, the calculated co-adsorption free energies of Coads_P, Coads_M and Coads_O were -67, -44 and -42 kJ/mol, respectively.

Different co-adsorption geometries may evolve into various transition states, leading to the formation of PX, MX and OX, as described in Fig. 1(d)-(f), respectively. The methyl group was almost planar with the benzene ring in the transition state structure, and the oxygen atom of the water molecule, the carbon atom of the methyl group and the neighboring carbon in the benzene ring were in a line. The ring carbons closest to the carbon atom of the methyl group were in the para, meta and ortho carbons, leading to the formation of PX, MX and OX, respectively. The lengths of the C-C bonds being formed in these cases were 2.12, 2.10 and 2.09 Å, whereas the lengths of the C-O bonds being broken were 2.13, 2.18 and 2.14 Å, respectively. As shown in Fig. 2, the free energy barriers for the formation of PX, MX and OX were 167, 138 and 139 kJ/mol, which were similar to the electronic energy barriers previously reported by Vos et al. ^[20] (127, 151 and 142 kJ/mol) for the zeolite mordenite and Li et al. ^[19] (123, 144 and 134 kJ/mol) for the zeolite HZSM-5. Our results therefore represent the highest energetic barriers for the formation of PX in HZSM-5 compared with those reported by Li et al. ^[19] and Vos et al. ^[20]. This difference can be attributed to the varied natures of the model and functional. Moreover, different co-adsorption geometries were used in the current study for toluene and methanol, which were optimized for the varied orientations of the methyl group in the toluene molecule. In contrast, the co-adsorption geometries used for the formation of PX, MX and OX were the same in the study of Li et al. ^[19]. The resulting co-adsorption geometry for the formation of PX was the most stable of all of the geometries determined in the current study.

The carbon atom of methanol could attack the benzene ring in toluene to give a series of C₈H₁₁⁺ intermediates (depending on the regioselectivity of attack), which would be rapidly deprotonated, yielding the corresponding xylenes with the regeneration of the Brönsted acidic site. MX was found to be the most thermodynamically favored of the adsorbed products, because it was 36 and 15 kJ/mol more stable than the adsorbed PX and OX products, respectively (Fig. 2). Notably, this result was in agreement with the predictions reported by Li et al. ^[19]. This result could be attributed to the geometry of MX being the best fit for the channel of HZSM-5, in the sense that MX could position its methyl groups in the straight and sinusoidal channels, respectively, thereby lowering its adsorption energy.

The stepwise pathway for the methylation of toluene involved the initial formation of a methoxy group and its subsequent use to achieve the methylation of toluene. A free energy diagram for the formation of the methoxy group is shown in Fig. 3. The first step would involve the adsorption of methanol in the acidic site. According to this model, methanol would be adsorbed in an end-on adsorption mode and stabilized by two H-bonds to HZSM-5 with lengths of 1.50 and 1.37 Å. The free energy for the adsorption of methanol over HZSM-5 was determined to be -22 kJ/mol. According to the transition state structure, the bond lengths of the C-O bonds being broken and formed were 1.93 and 2.09 Å, respectively. The free energy barrier for the dehydration of methanol in this way was determined to be 145 kJ/mol, which is consistent with the electronic energy barriers obtained in HZSM-22 (139 kJ/mol) ^[18] and mordenite (150 kJ/mol) ^[16]. Following the formation of the methoxide intermediate, a single water molecule would be adsorbed in close proximity to this species. However, we considered that the co-adsorbed water molecule would not interact with the reacting species during the stepwise pathway, and it was therefore concluded that the water molecule would be desorbed prior to the subsequent methylation step.

After the formation of the surface methoxy group, a toluene molecule would enter the pore of the HZSM-5 catalyst and be adsorbed on to its surface in close proximity to the methoxy group. As mentioned above, this would result in various geometries with different orientations of the methyl group in toluene, which were denoted as Ads_p, Ads_m and Ads_o (Fig. 4). Although the geometric parameters of zeolite framework were similar in these structures, the distances between the carbon atom of the incoming methyl group and the ring carbons of toluene varied considerably. The closest ring carbon to the carbon atom of the incoming methyl group was in the para, meta and ortho carbon position of toluene, which proceeded to form different transition states, resulting in the desorption of PX, MX and OX as products, respectively. The corresponding adsorption free energies of Ads_p, Ads_m and Ads_o were determined to be -21, -18 and -10 kJ/mol, respectively (Fig. 5).

According to the transition state structures shown in Fig. 4(d)-(f), the methyl groups would be almost co-planar and parallel to the adsorbed toluene molecules. The lengths of the C-O bonds being broken in TS_p, TS_m and TS_o were 2.08, 2.13 and 2.10 Å, whereas the lengths of the C-C bonds being formed were 2.22, 2.17 and 2.19 Å, respectively. As shown in Fig. 5, the free energy barriers for the formation of PX, MX and OX were determined to be 127, 105 and 106 kJ/mol, respectively, suggesting that the formation of PX has the highest energy barrier.

Having overcome the energy barriers, the methyl groups would proceed to attack the toluene molecules, leading to the formation of C₈H₁₁⁺ species. As shown in Fig. 6(a)-(c), the charged intermediates, Int_p, Int_m and Int_o, would be adsorbed at the intersection of the HZSM-5 channels with their methyl groups positioned in close proximity to the acidic site. To obtain the corresponding xylenes, the charged intermediates would have to back-donate a proton to the zeolite framework. As shown in Fig. 6, this process would require the rotation of Int_p, Int_m and Int_o to give Int rot_p, Int rot_m and Int rot_o, which would provide the necessary geometries for proton back-donation. Jobic et al. ^[37] reported that the energy barrier for the rotation around the molecular C₂ axis of benzene inside an 8 Å silicalite was between 25 and 105 kJ/mol, with an average value of 60 kJ/mol. Vos et al. ^[20] estimated that the reorientation energy barrier of p-C₈H₁₁⁺ over mordenite was around 15 kJ/mol, which would be easy to overcome, making the reorientation of p-C₈H₁₁⁺ very fast. With this in mind, it was envisaged that the reorientation energy barriers of Int_p, Int_m and Int_o over HZSM-5 would be easy to overcome. However, the confinement effect of the zeolite framework could potentially prevent the reorientation of these species ^[20].

Hindering the reorientation of these intermediates would result in the movement of proton H2 towards the methyl group because of the effect of the negative charge in the acidic site, which would result in the transition state geometries shown in Fig. 7. According to these geometries, the carbon atom of the methyl group would be simultaneously connected to four protons, which could lead to the cleavage of the C1-C2 bond and the formation of methane. The calculated free energy barriers for the formation of methane from Int_p, Int_m and Int_o were 136, 132 and 134 kJ/mol, respectively, which are much higher than those required for the reorientation of these species (Table 1). The reorientations of the Int_p, Int_m and Int_o would therefore be energetically favored over the formation of methane. Given that the reorientation energy barriers were low, the protons of the intermediates could readily back-donate to the zeolite framework, releasing the corresponding xylenes and recovering the acidic sites. This would complete the catalytic cycle, allowing it to run without deactivation. Lercher et al. ^[38-40] conducted a series of experiments, which showed that methane was mainly formed via the demethylation of polymethylaromatics, with methanol acting as a hydride donor leading to the formation of methane, which is consistent with our results.

As well as being methylated at its para, meta and ortho positions, toluene can also be methylated at its geminal position. The main geometries for the geminal methylation of toluene are shown in Fig. 8. The lengths of C-O and C-C bonds being broken and formed in the concerted pathway were 2.24 and 2.06 Å, respectively, while they were 2.16 and 2.20 Å in the stepwise pathway. The geminal methylation of toluene would lead to the formation of the geminal C₈H₁₁⁺ species shown in Fig. 8(d). The para, meta, ortho and geminal C₈H₁₁⁺ species could all be formed via an alkyl transfer process, which would be strongly hindered based on the confinement effects of the zeolite framework. For the geminal C₈H₁₁⁺ species, it would be almost impossible to position a proton in such a way as to achieve back-donation to the negative bridging oxygen. Strong electrostatic effects between the zeolite framework and the geminal C₈H₁₁⁺ species could therefore lead to the mutual approach of both methyl groups, leading to the formation of ethene or ethane, in a similar manner to the hydrocarbon pool (HCP) mechanism ^[41-44].

The free energy barrier for the geminal methylation of toluene via a concerted pathway was determined to be 183 kJ/mol, which is higher than those of the para, meta and ortho methylation reactions. However, the free energy barrier for the geminal methylation of toluene via a stepwise pathway was 114 kJ/mol, which is similar to those of the corresponding para, meta and ortho methylation reactions. These results therefore indicate that the geminal methylation of toluene would proceeded much more readily via a stepwise pathway. These results also indicate that the para, meta, ortho and geminal C₈H₁₁⁺ species would be generated in equal amounts via a stepwise pathway based on the similarities in their free energy barriers. However, geminal C₈H₁₁⁺ would be unable to back-donate a proton to the zeolite framework or diffuse out of the pores of the catalyst, leading to the formation of ethene or ethane.

The enthalpies and free energies of the adsorption geometries involved in this study are summarized in Table 2. The adsorption enthalpy of methanol was determined to be -128 kJ/mol, which is slightly smaller than that of the experimental values (-115 ± 5 kJ/mol) ^[45]. However, this result was consistent with the value reported by Mynsbrugge et al. ^[46] (-125 kJ/mol), who used the same method as the one used in this study. The calculated free energy for the adsorption of methanol was -22 kJ/mol, which is smaller than the value (8 kJ/mol) reported by Brogaard et al. ^[18] over HZSM-22. This deviation could be attributed to differences in the functional used for the calculation, which, in contrast to the energy barriers, can be quite sensitive to the methanol adsorption energy ^[18].

The rate-determining step for the stepwise pathway was determined to be the dehydration of methanol, which had a free energy barrier of 145 kJ/mol. This value was similar to those obtained for the concerted pathway, making it consistent with the results reported by Li et al. ^[19]. The differences in the energy barriers for the two pathways were therefore too small to allow us to distinguish between them, because diffusion plays a significant role in this reaction. The calculated order of kinetic selectivity for the formation of the different xylenes was para < ortho ≈ meta for both pathways, which is different from those obtained by Li et al. ^[19] and Vos et al. ^[20] based on electronic energies. However, experimentally, this process exhibits a greater selectivity for PX over MX and OX, which could be attributed to the high diffusion coefficient of PX, which is about 100 and 1000 times higher than those of OX and MX above 523 K, respectively ^[5]. Furthermore, MX and OX experience stronger diffusional resistance compared with PX, and the increased diffusional resistance of these species would lead to increases in their residence time relative to PX. Consequently, the selectivity of this process for PX could be improved by enabling the isomerization of the MX and OX to PX ^{[1, 47]}. The activation energies of these systems also suggest that the xylenes would isomerize at a much faster rate than toluene alkylated by methanol ^[48].

The HCP mechanism, which involves the formation of hydrocarbons via the repeated methylation and cracking of olefins and aromatic compounds encapsulated in the pores of a zeolite, is widely accepted as an intrinsic part of the MTH process. Most notably, the HCP mechanism provides a reasonable interpretation of the induction period during the early stages of the MTH process and can be used to explain the formation of the major products of this reaction ^{[41, 43, 44, 49-54]}. The products formed during the methylation of toluene with methanol can also be explained by the HCP mechanism. However, the methylation of toluene with methanol does not require an induction period to form HCP species, because the initial HCP species already exist in the reactants. The methylation of toluene would lead to the formation of methane, ethylene, propene and several other gaseous products via the cracking of olefins from the HCP species, which is consistent with the experimental results reported by Ahn et al. ^[55]. Furthermore, benzene is an important byproduct of the methylation of toluene, which is generally produced by the disproportionation of toluene ^{[47, 56]}. However, the dealkylation of the g-C₈H₁₁⁺ intermediate to give ethene would lead to the formation of a C₆H₇⁺ species, which would be capable of rapidly back-donating a proton to the zeolite framework, forming benzene. Benzene could also be produced by the dealkylation of several other carbocations.

The reaction coordinates for the methylation of toluene with methanol over HZM-5 were investigated using periodic density functional theory for the concerted and stepwise mechanisms. Various adsorption geometries were obtained with a diverse range of orientations for the methyl group of the toluene molecule, which eventually led to the formation of PX, MX and OX. The calculated free energy barriers for the formation of PX, MX and OX via the concerted pathway were 167, 138 and 139 kJ/mol. The rate-determining step for the stepwise pathway was determined to be the dehydration of methanol, which had a free energy barrier of 145 kJ/mol. The free energy barriers for the formation of PX, MX and OX via the stepwise pathway were 127, 105 and 106 kJ/mol, respectively. Moreover, the results showed that the formation of PX was kinetically disfavored by both pathways. C₈H₁₁⁺ species were formed as important intermediates in the concerted and stepwise pathways, which could back-donate a proton to the zeolite framework via a reorientation step or form gaseous products through demethylation. Methane was formed via an intramolecular hydrogen transfer from a ring carbon to the carbon of the methyl group, with calculated free energy barriers of 136, 132 and 134 kJ/mol for the para, meta and ortho C₈H₁₁⁺ species, respectively.

The methylation of toluene at the geminal carbon position would lead to the formation of a g-C₈H₁₁⁺ species, which would be unable to back-donate a proton to the zeolite framework, leading to the formation of ethene or ethane. The energy barrier for the geminal methylation of toluene via the concerted pathway was 183 kJ/mol, which was higher than that of the stepwise pathway, making the stepwise pathway the preferred pathway. The utilization rate of methanol is less than 100% for the methylation of toluene because of the formation of gaseous products from the C₈H₁₁⁺ species.