Mechanistic Characterization of Zeolite-Catalyzed Aromatic Electrophilic Substitution at Realistic Operating Conditions.

Zeolite-catalyzed benzene ethylation is an important industrial reaction, as it is the first step in the production of styrene for polymer manufacturing. Furthermore, it is a prototypical example of aromatic electrophilic substitution, a key reaction in the synthesis of many bulk and fine chemicals. Despite extensive research, the reaction mechanism and the nature of elusive intermediates at realistic operating conditions is not properly understood. More in detail, the existence of the elusive arenium ion (better known as Wheland complex) formed upon electrophilic attack on the aromatic ring is still a matter of debate. Temperature effects and the presence of protic guest molecules such as water are expected to impact the reaction mechanism and lifetime of the reaction intermediates. Herein, we used enhanced sampling ab initio molecular dynamics simulations to investigate the complete mechanism of benzene ethylation with ethene and ethanol in the H-ZSM-5 zeolite. We show that both the stepwise and concerted mechanisms are active at reaction conditions and that the Wheland intermediate spontaneously appears as a shallow minimum in the free energy surface after the electrophilic attack on the benzene ring. Addition of water enhances the protonation kinetics by about 1 order of magnitude at coverages of one water molecule per Brønsted acidic site. In the fully solvated regime, an overstabilization of the BAS as hydronium ion occurs and the rate enhancement disappears. The obtained results give critical atomistic insights in the role of water to selectively tune the kinetics of protonation reactions in zeolites.

Introduction

Acid-catalyzed (de)alkylation reactions of aromatic substrates constitute a key class of organic chemistry processes. Presumably, one of the most remarkable examples is represented by the zeolite-catalyzed alkylation of benzene with ethene to produce ethylbenzene,[1,2] which is the first step in the synthesis of styrene monomers for polystyrene production.[3] Quite some interest was initially also directed to the use of ethanol as an alkylating agent; however, the low prices of fossil resources in the previous century made this alternative poorly appealing.[4] Recently, on the other hand, the relevance of the use of ethanol is growing together with the search for environmentally friendly processes, as ethanol may be obtained from renewable resources.[5] Moreover, zeolites have also been shown to effectively catalyze the dealkylation of the alkylphenolic monomers derived from lignin valorization processes to produce simpler and more useful molecules such as phenol and olefins.[6,7] Therefore, zeolite-catalyzed (de)alkylation reactions of aromatic substrates are regaining a great deal of attention as an effective tool in the conversion and valorization of biomasses to commodity chemicals.

Mechanistically, the prototypical alkylation of benzene with ethene or ethanol in the H-ZSM-5 zeolite has represented the main case study in the computational investigation of zeolite-catalyzed alkylation reactions.[8−15] There is consensus that two main reaction pathways are possible: (i) a stepwise mechanism (TS0, TS1, and TS3 in Figure ) and (ii) a concerted mechanism (TS2 and TS3 in Figure ). In the former, ethene or ethanol is first chemisorbed on the zeolite walls as a surface ethoxide species (SES, Figure b). The SES then acts as an electrophile in a typical electrophilic aromatic substitution reaction (SEAr), which generates a Wheland complex (also known as σ complex, Figure c) as the intermediate, i.e., a protonated arenium ion. In the concerted mechanism, ethene or ethanol is directly activated by the zeolite Brønsted acid site (BAS) and undergoes the SEAr reaction to attack the aromatic substrate, again forming the Wheland complex. The latter can then transfer a proton back to the zeolite framework, giving the final ethylbenzene product (Figure d).

Figure 1

Schematic depiction of all reactions investigated in this work. (Top) Possible mechanistic pathways for the ethylation of benzene with ethanol (orange arrows) and ethene (blue arrows) are shown. From left to right, the pristine zeolite BAS (a) can directly react with ethene or ethanol (TS0) to give a chemisorbed surface ethoxide species (b) which, in its turn, can attack benzene (TS1) to possibly form the ipso-protonated Wheland complex (c). Alternatively, benzene and ethene or ethanol can directly react to again form the ipso-protonated Wheland complex (TS2). The latter can then deprotonate (TS3), leading to the final ethylbenzene product (d). Ethylbenzene can also be further reprotonated by the framework on the para carbon, forming a second Wheland complex (e); this last step was investigated in the presence of n = 0, 1, 3, 6 water molecules (TS3nw).

Interestingly, the formation of the Wheland complex—not only during benzene ethylation but in all SEAr reactions—is still a matter of debate.[16] While in homogeneous catalysis its existence has been recently questioned,[16] the microporous environment of the zeolite is known to substantially stabilize charged species,[17,18] thus making its presence more likely. From a computational perspective, initial studies on cluster models indicated that the deprotonation step was occurring spontaneously and the Wheland complex was not a relevant reaction intermediate.[8,9] However, later studies showed that its stability is strongly dependent on the size of the cluster adopted in the calculation, becoming an effective minimum on the potential energy surface (PES) once the zeolite environment was better accounted for through the use of bigger cluster models.[11] As could be expected, with the subsequent passage to more realistic periodic models of the zeolite framework, the Wheland complex naturally emerged again as a stable product state of the SEAr step, although its stability has not been deeply investigated so far.[12−14] On the other hand, while studying the effect of the zeolite topology on ethylbenzene transalkylation, Corma and co-workers found that the reaction can start from the ipso protonation of the substrate with an expected free energy increase of 44 kJ·mol–1 in the channel intersection of H-ZSM-5 at 573 K.[19,20] The authors also showed a remarkable dependence of the Wheland complex stability with respect to the chosen framework topology and location of the substrate within it.

Some additional insights can also be derived from the methanol-to-hydrocarbons process, where methylation of the aromatic hydrocarbon pool species also proceeds through a SEAr mechanism.[21] Using static DFT calculations with higher order corrections, Fečík et al.[22] recently found the Wheland intermediate as a minimum in the PES when benzene is methylated with methanol or dimethyl ether. The deprotonation was found to proceed with a small barrier of 47 kJ·mol–1 and a large energetic gain of 127 kJ·mol–1, indicating a very short lifetime for the intermediate.

Experimentally, the existence of the Wheland complex has been recently supported by Chowdhury and co-workers.[23] Through a series of operando UV–vis and NMR spectroscopy measurements, they were able to identify the Wheland complex and the SES as key intermediates present during the alkylation of benzene with ethanol over H-ZSM-5.

To thoroughly investigate the reaction mechanism at operando conditions, i.e., at elevated temperatures and in the presence of water, it is necessary to go beyond the standard static approach where only a few points on the potential energy surface are considered. In addition, the simple harmonic approximation normally adopted to evaluate free energies has proven to be often inadequate to correctly describe the large anharmonic motions of the adsorbates in the zeolite pores.[24−28] This especially reflects in the quality of the reaction entropy estimate and the mobility of the reactive species whose weight rapidly grows with the reaction temperature.[29,30]

In this work, we used enhanced sampling molecular dynamics (MD) simulations to overcome the intrinsic limitations of static methodologies and investigate the alkylation of benzene with ethene and ethanol at realistic operating conditions. First, two-dimensional multiple-walkers[31] well-tempered[32] metadynamics was used to freely explore the free energy landscape of the SEAr reaction with ethene and the SES as ethylating agents, proving that the Wheland intermediate spontaneously appears as a minimum at operando conditions. Second, we fully characterized all of the alkylation transition states with umbrella sampling (US) and used classical transition-state theory to retrieve accurate kinetic constants for all reaction steps (which are graphically shown in Figure and explicitly listed in Table ). In that way, we showed that the use of enhanced sampling can remarkably change the activation energy of mobile transition states and possibly alter the mechanistic conclusions derived from static calculations. Finally, we investigated the para protonation of ethylbenzene (being energetically more favorable than the ipso one in the SEAr) in the presence of various loadings of water in the zeolite unit cell. 0 (TS30w), 1 (TS31w), 3 (TS33w), and 6 (TS36w) water molecules were considered. In this way, we showed that water can act as a proton-shuttling agent and assist the protonation reaction. At the lowest coverage, this effect can increase the rate of proton exchange by about 1 order of magnitude. When higher loadings are considered, on the other hand, the BAS gets solvated by the water cluster as a hydronium ion. This effect tends to remarkably stabilize it, making the proton transfer kinetics slower and similar again to the anhydrous case at the higher considered loading.

Table 1

Nomenclature of the Reaction Steps Investigated in This Work

name	description	reactants	products
TS0OH	SES formation from ethanol	ethanol	SES + water
TS0=	SES formation from ethene	ethene	SES
TS1	SES electrophilic attack on benzene	SES + benzene	ipso-protonated Wheland intermediate
TS2OH	ethanol electrophilic attack on benzene	ethanol + benzene	ipso-protonated Wheland intermediate + water
TS2=	ethene electrophilic attack on benzene	ethene + benzene	ipso-protonated Wheland intermediate
TS3i	ipso-protonated Wheland intermediate deprotonation	ipso-protonated Wheland intermediate	ethylbenzene
TS3pnw	para protonation of ethylbenzene with water	ethylbenzene + n water	para-protonated Wheland intermediate + n water

Computational details

Catalyst Model

To represent the three-dimensional MFI pore structure of the H-ZSM-5 zeolite catalyst, we adopted a periodic model with a unit cell containing 96 tetrahedral Si atoms connected by oxygen bridges (Figure S1). As is commonly done,[33,34] the catalytic active site was modeled by substituting a Si4+ atom at the T12 position (located at the channel intersection) with Al3+ to give an Si/Al ratio of 95 and, therefore, isolated active sites. The negative charge created by the substitution was compensated with the addition of a proton on one of the oxygen atoms adjacent to the defective site (Ozeo,1 in Figure S1). In this way, the proton, i.e., the actual active BAS for the alkylation reaction, finds itself located at the intersection between the straight and sinusoidal channels of the MFI pore structure, ensuring maximal accessibility for the substrate. The initial location of the BAS is of limited importance as all oxygens in the first coordination sphere of the Al defect are equivalent in the collective variable definition (vide infra), allowing the BAS to spontaneously interact with the energetically most favorable site.

To model the reactions, a single molecule of each involved species was manually placed in the proximity of the BAS as initial structure for the MD simulations, using our in-house developed software Zeobuilder.[35] This means, for instance, a benzene and an ethene molecule were used for the reactant state of TS2= or a single ethanol molecule was used for the reactant state of TS0OH (Table ).

Ab Initio Molecular Dynamics

All of the simulations performed in this work are based on enhanced sampling techniques, which rely on biased ab initio MD to derive the free energy profile for a chemical reaction (vide infra). This allows us to properly account for entropic effects by accurately simulating the dynamic behavior of the adsorbate molecules as well as accounting for the catalyst flexibility at operating conditions. Born–Oppenheimer MD simulations were performed using the CP2K software[36] (version 5.1) within the density functional theory (DFT) framework. We adopted the PBE exchange-correlation functional[37] in its parametrization revPBE[38] (due to its improved performance for reaction energies calculations[39]), coupled with Grimme’s D3 dispersion correction[40] to account for long-range dispersive interactions. A triple-ζ quality basis set including valence and polarization functions was adopted for all atoms in the atom-centered Gaussian orbitals + plane waves (GPW)[41,42] basis set approach used by CP2K. GTH pseudopotentials[43] were used to smooth the electron density in the proximity of the nuclei, and the plane waves energy cutoff was set to 350 Ry. The time step for the integration of the equations of motion was set to 0.5 fs. All simulations were conducted in the NPT ensemble using a chain of five Nosé–Hoover thermostats[44,45] and an MTK barostat[46] to control temperature and pressure, respectively, that were set to 573 K and 1 atm to reproduce the experimental conditions.[23] The NPT ensemble was chosen to allow maximum catalyst flexibility, but due to the rigid nature of the H-ZSM-5 structure, the small fluctuations in the unit cell parameters are expected to have little influence on the results.

Collective Variables Definition

Since in regular MD simulations the chances of sampling activated events with a barrier larger than few times kBT (like most chemical reactions) is extremely small, enhanced sampling techniques are used to bias the system along a set of properly chosen collective variables (CVs) and retrieve the activation energy of the process of interest. A CV is selected to monotonically vary in the function of the reaction progression from a value associated with the reactant state to a different value associated with the product state. All CVs considered in this work are function of the atomic coordinates of the system. Being interested in chemical reactivity, an ideal CV should be able to describe the formation/rupture of bonds while, at the same time, allow to account for the chemical equivalence of some atoms (e.g., the six aromatic carbons of benzene). Both goals can be achieved by using a linear combination of coordination numbers (CNs). The CN between two groups of atoms α and β is defined aswhere rij is the distance between atom i, belonging to group α, and atom j, belonging to group β. Unless differently specified, NN = 6 and MM = 2NN. The function argument of the double summation tends to 1 if rij < r0 and quickly drops to 0 otherwise. In our simulations, the parameter r0 is selected to be roughly equal to the length that a bond between two atoms of the α and β groups assumes in the transition state region. In that way, there will be a contribution of ∼1 to CN(α; β) if i and j are bounded (rij < r0TS) and of ∼0 otherwise, thereby making CN(α; β) approximately equal to the number of bonds between group α and group β.

Almost all the reactions investigated in this work involve the simultaneous formation and rupture of bonds between at least three (sometimes four) group of atoms. However, obtaining converged free energy surfaces (FESs) in more than two dimensions is extremely costly and often unfeasible.[47] To reduce the dimensionality of the problem, we selected our final collective variables to be a linear combination of CNs. In practice, two types of linear combinations were used: the first one (CVA) is the simplest and is used to describe a reaction in which one atom (whether or not it belongs to a larger molecular fragment) is transferred from a group of atoms to another. This is the case, for instance, in TS3, where a single proton is transferred between the ipso aromatic carbon and the zeolite oxygens. For this kind of process, a simple difference of CNs can appropriately describe the reaction (eq ).In practice, an atom belonging to group β (simply the proton in the TS3 example) is transferred from group α (the ipso aromatic carbon) to group γ (the four oxygen atoms around the Al site) or vice versa, adding 1 to one of the CNs while removing 1 to the other. Therefore, this type of CV will vary of about two units between reactants and products (see Figure S2 and CV2 in Figure ).

Figure 2

(a, b) Definition of the coordination numbers (CNs) and collective variables (CVs) used in the metadynamics simulations of TS1/TS3 (a) and TS2=/TS3 (b). (c) Schematic depiction of the expected features of the two-dimensional MTD free energy surface as a function of CV1 (describing the electrophilic attack on benzene) and CV2 (describing the deprotonation of the Wheland complex). The location of reactants (R), products (P), and Wheland complex (W) is shown. The orange arrows follow a hypothetical path going through the formation of a stable Wheland complex, while the blue one shows a concerted conversion from R to P.

More complex is the specific case in which ethene is reacting, either to form the SES intermediate (TS0=) or directly attacking benzene (TS2=). In this case, two processes are occurring simultaneously: a proton is transferred from the zeolite to one of the ethene carbon atoms while the other carbon attacks the nucleophilic moiety (either another zeolite oxygen or the benzene molecule). To avoid an excessive increase in the number of dimensions, the following linear combination of CNs was used as collective variable:TS2= can be taken as an example to explain this choice (see Figure b). First, ethene is activated by the zeolite BAS, that protonates it. Similarly to the previous example with TS3, the reaction can be described with a difference in coordination numbers (CN(β; γ) – CN(γ; δ) in eq ). In this case, β includes the two aliphatic carbons while γ the four aliphatic hydrogens and the BAS, as they become indistinguishable once the protonation has occurred. CN(β; γ) is then going to vary from about 4 (the number of C–H bonds in ethene) to 5 (the number of C–H bonds in an ethyl fragment). δ includes the four oxygen atoms surrounding the Al site of the framework, and therefore, CN(γ; δ) goes from 1 (BAS on the framework) to 0 (BAS on the ethyl fragment). Overall, CN(β; γ) – CN(γ; δ) is then equal to about 3 in the reactant region and 5 in the product region. If α is chosen to include the six aromatic carbons, CN(α; β) can describe the C–C bond formation between benzene and the ethyl fragment, with an expected variation from 0 to 1. The dimensionality of the FES can then be reduced by summing the two partial CVs but, since CN(β; γ) - CN(γ; δ) varies of two units while CN(α; β) only of one, the former is first multiplied by 0.5. This should allow for a smoother variation of the CV from reactants to products. A graphical depiction of the aforementioned steps can be seen in Figure S3.

Every collective variable used in the advanced sampling simulations either falls in one of the two previous categories or simply is a single coordination number. More specific information is reported in the following sections.

Well-Tempered Metadynamics

Well-tempered metadynamics (MTD) and umbrella sampling (US) simulations are the two advanced sampling techniques that were chosen to explore the FES of benzene ethylation by biasing the system in a one- or two-dimensional CV space. For both techniques, the bias was applied to the MD simulations using PLUMED[48] (version 2.4.0) as dependency of CP2K. Initially, MTD was used to explore the FES of the SEAr step in benzene ethylation with a SES (TS1 in Figure ) and with ethene (TS2=) and assess whether the Wheland complex spontaneously appears as reaction intermediate. To make the formation of the Wheland complex possible—but not mandatory—a two-dimensional CV space is needed in which the first CV describes the electrophilic attack on the benzene ring while the second the deprotonation of the Wheland complex (Figure c). In the case of TS1, the electrophilic attack can be described with a difference of coordination numbers as shown in eq , by defining α = Car (the six aromatic carbons), β = Cal (the reactive alkyl carbon), and γ = Ozeo (the four oxygen atoms in the first coordination sphere of the Al defect, Figure a). For TS2=, the more complex linear combination of CNs described by eq must be used with, as explained before, α = Car, β = Cal (this time considering both the alkyl carbons), γ = Hal (the alkyl hydrogens including the BAS), and δ = Ozeo (Figure b). For TS3, a difference of CNs (α = Car, β = Har, i.e., the six aromatic hydrogens and γ = Ozeo) can again be used. By using MTD to explore this 2-dimensional CV space, the stability of the Wheland complex can be elucidated. Indeed, if its formation is spontaneous during the reaction, a minimum in the FES will be sampled at the corresponding CVs values (orange arrows in Figure c), while a direct interconversion between reactants and products will be observed otherwise (blue arrow in Figure c).

Concerning the methodology, two multiple-walkers[31] well-tempered[32] metadynamics[49,50] simulations were performed. In this approach, the bias potential is constructed on the fly by spawning a Gaussian-shaped hill at regular time intervals, centered on the average CV position explored by the walker in a certain number of previous simulation steps. Every walker feels, at every simulation step, an overall potential given by the sum of all walkers’ spawned hills up to that moment. Six walkers were run in parallel, three initially located in the reactant basin (benzene + ethene for TS2= or benzene + SES for TS1, see Table ) and three in the product (ethylbenzene) basin. The hills were spawned in the CVs space every 50 fs, with an initial height of 5 kJ·mol–1. Such height was then gradually rescaled according to the well-tempered recipe, with an initial BIASFACTOR—as defined by PLUMED—of 10, subsequently increased to 15 to improve convergence. The mass of all hydrogen atoms was set to 2 to improve the stability of the simulations. The simulations were stopped once the hills height dropped below 1 kJ·mol–1, and a few barrier recrossings were observed among the various walkers. The final estimate of the FES was then obtained by inverting the sum of the hills spawned by all walkers. A set of walls was required to improve barrier recrossing or avoid the occurrence of side reactions, the details of which are given in the Supporting Information, Section S2.4.

Umbrella Sampling

Reaching a satisfactory convergence of the FES estimate in two-dimensional MTD is extremely difficult and requires a prohibitively long computational time. For this reason, we subsequently used one-dimensional umbrella sampling[51,52] (US) simulations to derive the FES of all the benzene ethylation reaction steps, both with ethene and ethanol (see Figure ). Moreover, US was also used to investigate the role of water in the para protonation of ethylbenzene when various loadings are considered (TS3nw, n = 0, 1, 3, 6). The adopted CVs are fundamentally analogous or precisely the same to the one presented for MTD in Figure , belonging to one of the two categories described by eqs and 3. An exception is represented by TS3nw (n ≠ 0), in which a single CN between the ethylbenzene para carbon and all hydrogen atoms involved in the protonation process was sufficient to sample the reaction. A full list of the US CVs is reported in Table S1 and graphically depicted in Figure S4.

In a US simulation, various quadratic potentials (the “umbrellas”) are set along the selected CV from the reactants to the products basin and a MD simulation is then run in each of them. The bias has the formwhere V(CV) is the bias of the ith umbrella, κ is its spring constant, and CV0, is the collective variable value at which it is centered. The κ and CV0, choice was guided by previous literature reports[30,53] and tuned to balance a uniform sampling between reactant and products while keeping low the number of required umbrellas. Being an equilibrium technique and, therefore, intrinsically more stable than MTD, the hydrogen mass in US was set back to 1. The FES estimate was obtained by combining the simulations through the weighted histogram analysis method[54,55] (WHAM), as implemented in our in-house developed ThermoLIB library.[56] A full list of the umbrella parameters used in the simulations for all transition states is reported in the Supporting Information, Tables S2–S9.

While using a linear combination of coordination numbers allows us to perform 1-dimensional US simulations, which is much faster to converge than multidimensional US, the risk of poorly exploring some important regions of the phase space increases the more degrees of freedom are “squeezed” in a single linear combination. To ensure that an appropriate sampling of the relevant phase space regions in every reaction was achieved, the 1-dimensional free energy profiles have been expanded in terms of their constituting coordination numbers using the statistical analysis tools of ThermoLIB,[56] to obtain 2-dimensional FESs. More details on the procedure can be found in the Supporting Information, Section S2.3. While for most reactions a good sampling of the reaction path was observed, it was noticed that some possibly important regions around the transition state of TS0= were poorly sampled. To solve the problem, some extra 2-dimensional umbrellas were added to reach a satisfactory coverage of the transition state region and the complete 2D FES was then reprojected onto the original 1D CV (more details are reported in Section S2.3 of the Supporting Information).

As for the MTD case, most of the simulations required some walls to prevent undesired side reactions or facilitate convergence by improving barrier recrossing. Full details are reported in the Supporting Information, Section S2.4.

Calculation of the Phenomenological Reaction Barriers

As recently shown by some of us,[53,57] the information derived from an US simulation can be combined to classical transition state theory (TST) to remove the dependency on the chosen CV and retrieve reliable kinetic constant for the reaction of interest. In this framework, the reaction kinetic constant can be written asin which β = 1/kBT and CV⧧ is the value of the collective variable corresponding to the reaction transition state. ⟨|∇⃗ CV|⟩ CV is the ensemble average of the CV gradient with respect to the mass-weighted coordinates of the system, computed when the CV is restricted at the transition state value. From the kinetic constant, it is then possible to use Eyring’s equation in order to retrieve a so-called phenomenological reaction barrier no longer dependent on the CV choicewith h being Planck’s constant. Statistical analysis was used to compute the 95% confidence interval for the kinetic constant and the phenomenological barrier, more details can be found in Section S2.5 of the Supporting Information.

Results and Discussion

Probing the Wheland Complex Formation within the Zeolite Channels

To assess whether the Wheland complex is part of the reaction intermediates during the ethylation of benzene at operando conditions, we investigated the SEAr step with 2-dimensional multiple-walkers MTD considering both the SES (TS1/TS3) and ethene (TS2=/TS3) as alkylating agents to probe the effect of a mobile vs chemisorbed electrophilic agent.

By adopting two specific CVs, one describing the electrophilic attack while the other the deprotonation step (vide supra, Figure c), the system can evolve from reactants to products with the formation of the Wheland complex being possible but not mandatory. With each of the 6 parallel walkers in the two simulations running for about 70 ps, a total of more than 420 ps of total simulation time per mechanism was necessary. The final FESs are shown in Figure .

Figure 3

Free energy surfaces related to the alkylation of benzene with an SES (a) and with ethene (b), as obtained from MTD simulations. Isoenergetic lines are placed every 10 kJ·mol–1. Black crosses show the minimum free energy path (MFEP) connecting reactants and products, obtained according to the procedure of ref (58). Free energy values for the stationary points along the MFEP are reported in kJ·mol–1. For the CVs definition, see Figure .

Both profiles look qualitatively very similar. Interestingly, it is clearly visible how the SEAr proceeds in two distinct steps: first, the electrophile attacks the benzene ring (TS1 and TS2= in Figure a,b, respectively) with the formation of the Wheland complex. The raw reaction barrier, obtained from the minimum free energy path,[58] is higher for TS2= than for TS1 (117 vs 72 kJ·mol–1, respectively). This can be expected due to the preactivated nature of the SES with respect to ethene. The differences in activation energy will be further discussed in the following paragraph while analyzing the US results.

The Wheland intermediate is rather elusive and quickly deprotonates to ethylbenzene, restoring its aromaticity. Indeed, its corresponding minimum (top right of Figure a,b) is very shallow and a small barrier (10–25 kJ·mol–1, depending on the considered profile) is associated with its deprotonation, which is also largely exergonic (80–71 kJ·mol–1).

The differences observed in the TS3 barriers for the two cases are to be expected, as converging 2-dimensional free energy surfaces with MTD is extremely expensive. On the other hand, the purpose of these simulations was to freely explore the free energy surface of the SEAr reaction, from which it clearly emerged that the Wheland complex has to be considered as reaction intermediate. For this reason, 1D umbrella sampling simulations were subsequently performed on all the reaction steps, to retrieve accurate barriers and unravel the full reaction mechanism when ethene and ethanol are used as alkylating agents.

Reaction Mechanism of Benzene Ethylation

An overview of all the reactions investigated with US is shown in Figure and in Table . Six separate US simulations were performed, two for the formation of the SES (TS0= and TS0OH) and one for its electrophilic attack on the benzene ring (TS1), two for the direct ethylation of benzene (TS2= and TS2OH), and finally, one for the Wheland complex deprotonation (TS3). As previously stated, a full list of the adopted umbrellas and the walls adopted in each simulation can be found in the Supporting Information, Sections S2.2 and S2.4.

As previously explained, all reaction profiles were deprojected as a function of the single CNs constituting the final CV to ensure that the path connecting reactants and products was adequately sampled. All of the raw free energy profiles and the respective 2-D expansion are shown in Section S3.1 of the Supporting Information. Only in the case of TS0= was it found that the transition state region had possibly relevant portions of the phase space not properly sampled by the 1D umbrellas. For this reason, extra 2-dimensional umbrellas were added to improve the sampling and the final 2D FES reprojected on one dimension to obtain the final phenomenological barrier (Section S2.3).

A full overview of the reaction free energy profile can be seen in Figure a. Note that the alignment between reactant and product state of different reactions is done for graphical purposes only, as even for reactions with the same reactant/product state (for instance the ipso-protonated Wheland intermediate without coadsorbed molecules is formed both in TS1 and TS2=) there is no absolute guarantee that the same phase space was explored in both simulations.[59]

Figure 4

(a) Complete free energy profile for the ethylation of benzene with ethene and ethanol, as derived from ab initio umbrella sampling simulations. The numbers close to the curves report the respective reaction phenomenological barrier (in kJ·mol–1) together with the respective 95% confidence interval (also graphically shown by the shaded regions around the curves). The alignment of reactant and product states between different reactions is for graphical reasons only (see main text). (b) Cross-section of the zeolite model for the various reactions as seen along the sinusoidal channel, showing a representative snapshot arbitrarily extracted from the umbrella on top of the transition state. All zeolite atoms, except for the Al tetrahedra, are in white for the sake of clarity.

Starting from the stepwise path (TS0–TS1), the formation of the SES is the rate-determining step for both ethene and ethanol, which present a barrier of 101 ± 7 and 115 ± 3 kJ·mol–1, respectively. The formation of the SES from ethanol is therefore more difficult than from ethene as, even if the opposite extremes of the confidence interval are considered, a minimum difference of 4 kJ·mol–1 is present between the two barriers. At the reaction conditions, such difference would still imply that the formation of the SES from ethene is more than twice as fast as from ethanol. This is in line with previous reports in the literature,[13] in which it was shown that ethanol can adsorb more strongly on the BAS thanks to the formation of an H bond with its hydroxyl group, thereby stabilizing the reactant state. A strong interaction between the BAS and the hydroxyl group of the ethanol is also observed in the umbrella corresponding to the reactant state of TS0OH, where a hydrogen bond between ethanol and the BAS persists for the whole simulation time (Figure S14). The hydroxyl group hydrogen can also form a second hydrogen bond with the zeolite oxygens that is, however, significantly weaker (Figure S15). While potential walls are present in the simulations to explicitly prevent the molecules from leaving the active site, ethene appears to be quite more mobile and prone to diffuse away (Figure S16). Once formed, the SES can readily react with benzene, with a much lower barrier of 80 ± 5 kJ·mol–1, to produce the ipso-protonated Wheland intermediate as shown in the MTD simulations.

In the direct ethylation (TS2), ethene or ethanol are directly activated by the BAS and attack the benzene ring. Also in this case ethanol presents a slightly higher barrier than ethene (108 ± 8 vs 102 ± 3 kJ·mol–1, respectively), however, the difference is rather small certainly given the large error bars. These results are potentially counterintuitive, as both TS0 and TS2 consist in an activation of the electrophile (ethene or ethanol) with concerted attack on an electron-rich moiety (the zeolite oxygens or the benzene carbons). Therefore, one could expect a proportional difference in activation energies between ethanol and ethene, while for the former attacking benzene seems to be easier than attacking the framework. An explanation for this can be found in the transition-state geometries. Indeed, to form the SES (TS0OH in Figure b), ethanol is completely protonated by the BAS forming an EtOH2+ cationic species, which must then rotate to expose the electrophilic carbon toward the framework and undergo the SN2 reaction. On the other hand, when reacting with benzene, the proton transfer to the hydroxyl group can occur gradually while the reaction proceeds, as no large reorientation is needed to form the C–C bond (TS2OH in Figure b). This can also be seen by expanding the free energy profile in term of a new collective variable, encoding the proton transfer from the zeolite to the hydroxyl group (Figure S17). In the case of TS0OH, the BAS has been fully transferred from the zeolite to the hydroxyl group when the transition state is reached, while for TS2OH the transfer is still progressing and the partially positive H2O moiety remains in interaction with the zeolite framework. Similarly to what has already been observed in the methanol-to-hydrocarbon process, the presence of extra protic molecules could assist the proton transfer from the framework to the reacting ethanol molecule, facilitating the formation of the SES.[60,61] Therefore, the differences between ethanol and ethene for the two mechanisms could significantly change at high ethanol or water loadings.

The reaction barriers are in line with previous literature reports (for a full comparison see Table S13). While it is known that activation energies obtained with improved level of theory are in general higher,[12] an interesting comparison can be made with the static results reported by Wang et al.,[13] where PBE-D3 on a periodic H-ZSM-5 model is also used. The barriers for TS0=, TS0OH, and TS1 are in very good agreement, while a significant deviation can be seen for TS2= (∼−17 kJ·mol–1) and TS2OH(∼+18 kJ·mol–1). This can be explained considering that TS0 and TS1 are quite rigid transition states, involving either the formation of a bond with the zeolite framework (TS0) or the transfer of a framework-bounded species to the hindered benzene molecule (TS1). On the other hand, in TS2, the ethylating agent interacts contemporarily with the benzene and the BAS and quite a large mobility of the species is observed (Figure S18). Since the geometries explored by the molecules during the dynamic simulations do not dramatically deviate from the optimized transition state geometries of Wang et al.,[13] the relatively large difference in the computed barrier could be likely attributed to entropic effects arising from the large anharmonic motions of the adsorbates in the zeolite pores.

The change in activation energy between static and dynamic simulations can have an important effect on the preferred mechanism. While it is known that microkinetic modeling is needed to properly include all the reaction variables,[13,62] this is outside the scope of this work. On the other hand, the ratio between the forward kinetic constants of the stepwise and concerted pathways can already provide some interesting mechanistic insights. By considering the activation energies computed by Wang et al.,[13] TS2= is about an order of magnitude faster than TS0= while the opposite is true for ethanol (in line with the full results of the microkinetic model). On the other hand, by considering our barriers, we learn that both the stepwise and concerted mechanisms proceed with a similar speed when the same ethylating agent is considered, with ethanol likely favoring the concerted path (despite the partial overlap between the confidence intervals). It must be pointed out that Chowdhury et al.[23] observed the SES formation using an equimolar amount of benzene and ethanol in the reacting mixture. This is a quite high ethylating agent concentration, which for normal industrial applications is kept about 5 times lower than benzene to suppress the formation of polyalkylated products.[63] According to our results, the formation of the SES in the reaction environment does not derive from a strong energetic difference between the two mechanisms—as found by Wang et al.—but more likely from the statistical unlikeliness of having benzene and ethanol simultaneously coadsorbed on the BAS, which becomes more significant at low benzene concentrations. This highlights the importance of the use of advanced sampling methods in zeolite catalysis, where certain reactive events can exhibit significant entropic effects that are not possible to fully capture with the traditional static approach.[64]

Once the alkylation has occurred, the newly formed ipso protonated Wheland complex can easily deprotonate to give the final ethylbenzene product, with a barrier of only 18 ± 1 kJ·mol–1 and a large energetic gain (78 kJ·mol–1). These results point toward an extremely short lifetime for the Wheland complex, making it present in minimal concentrations in the reaction environment. It must be pointed out, on the other hand, that the ethylbenzene product carries an electron-donating substituent which should, in principle, facilitate the protonation of the ortho and para positions on the aromatic ring. Moreover, when ethanol is used as ethylating agent, water is formed during the reaction. Its effect on the protonation kinetics is rather unexplored and, therefore, we also investigated the para protonation reaction of ethylbenzene with various water loadings in the zeolite (0, 1, 3, and 6 molecules per unit cell). In this way we assessed whether–under specific conditions–the formation of Wheland complexes could become more favorable.

Water Influence on the Wheland Complex Formation

It is well-known that when water adsorbs in the H-ZSM-5 pores it can strongly interact with the BAS. When the number of water molecules per BAS exceeds 3, the proton can even become fully solvated as hydronium ion.[65−67] It goes without saying that this phenomena can have an impact on proton transfer reactions and, therefore, on the formation of the Wheland complex in the zeolite pores.

To assess the effect of water on the formation kinetics and stability of the Wheland complex we performed four separate US simulations with varying amount of water in the zeolite unit cell, namely 0, 1, 3, and 6 water molecules. Considering that the maximal amount of water that can adsorb on the BAS before condensation at 298 K is around 7–8 molecules in H-ZSM-5,[68] the chosen coverages should be representative for a wide range of water partial pressures. As previously stated, we focused on the para protonation of the ethylbenzene product, which should be significantly more favorable than the ipso protonation constituting the last step of the ethylation reaction.

The raw reaction profiles for the ethylbenzene para protonation in the presence of various amounts of water are shown in Figure a–c. The forward and backward phenomenological reaction barriers are listed in Table . Starting from the anhydrous case, the protonation on the para carbon exhibits—as expected—a significantly lower activation energy than the ipso one, going from 96 to 76 kJ·mol–1. The backward barrier, on the other hand, is computed to be exactly the same (18 kJ·mol–1), thereby making the para-protonated Wheland complex 20 kJ·mol–1 more stable than the ipso-protonated one. By looking at the effect of water, the forward barrier experiences a strong drop, from 76 ± 2 to 62 ± 3 kJ·mol–1, when one water molecule is introduced in the unit cell. By increasing the amount of water to 3 molecules per unit cell, the barrier increases to 68 ± 1 kJ·mol–1, reapproaching then the anhydrous one with the highest considered loading of 6 water molecules per unit cell (72 ± 2 kJ·mol–1). The backward barrier, on the other hand, experiences a much less pronounced variation, from a minimal value of 11 kJ·mol–1 for 3 water molecules per unit cell to a maximum of 18 kJ·mol–1 for 0 and 1 water molecules per unit cell.

Figure 5

(a–c) One-dimensional free energy profiles obtained from the umbrella sampling simulations of the p-ethylbenzene protonation in the presence of 1 (a), 3 (b), and 6 (c) water molecules per unit cell. (d–f) Percentage of occurrence of a certain water cluster size in proximity of the active site as a function of the reaction collective variable. (g–i) Two-dimensional free energy surfaces obtained by expansion of the original one-dimensional profile showing the free energy as a function of the distance between the para carbon and the proton location (see main text) and the minimum distance between the first-coordination sphere O atoms around the Al defect and the proton location. (j–l) Cross-section of the zeolite model, as seen along the sinusoidal channel, showing a representative snapshot arbitrarily extracted from the umbrella in proximity of the transition state. All zeolite atoms, except for the AlO4 tetrahedra, are in white for the sake of clarity.

Table 2

Forward and Backward Phenomenological Barriers (kJ·mol–1) for the Para Protonation of Ethylbenzene in the Presence of Various Amounts of Water in the Unit Cell

	no. of water/unit cell
	0 H2O	1 H2O	3 H2O	6 H2O
ΔF⧧f	76 ± 2	62 ± 3	68 ± 2	72 ± 2
ΔF⧧b	18 ± 2	18 ± 3	11 ± 2	14 ± 2

To better understand the reasons underlying these variations in the forward barriers, we analyzed the size of the water cluster being actively involved in the reaction as a function of the reaction CV. In other words, we determine for each simulation step in all umbrellas the amount of water molecules that are in close proximity of the active site, thereby giving an indication if the water—which is free to move around in the simulation cell—is actively taking part to the reaction. More information about the definition of the cluster size can be found in Section S5.1 of the Supporting Information. The cluster sizes were collected as a function of the reaction CV, and their fraction of occurrence is shown in Figure d–f. Some interesting observations can be made: first, it can be noticed that the cluster size tends to be significantly smaller in the product region than in the reactant one. For instance, in the simple case of 1 water molecule per unit cell, the fraction of samples in which the water is far from the active site (cluster size of 0) tends to increase going toward the products and the same can be said for both the 3 and 6 water molecules cases (dark blue bars in Figure d–f). This is caused by the strong affinity between the water molecules and the BAS which, in the reactant state, are free to interact. When the proton is transferred to the ethylbenzene molecule, on the other hand, its ability to form hydrogen bonds with water is strongly inhibited and the latter tends then to diffuse away from the Al defect. With water poorly interacting or even leaving the active site region in the product state, it can be expected that not much influence will be seen in the reactivity of the Wheland complex toward deprotonation, in line with the low variation in the backward phenomenological barriers (Table ).

Interestingly, it is also possible to notice how the transition state region for 3 water molecules (Figure e) is associated with a preference for larger cluster sizes than the reactant state. In the latter, cluster sizes of 1, 2, and 3 water molecules are similarly sampled in the simulation, but in the central region a large prevalence of a 3 molecules cluster is observed. In the reactant state, the BAS interacts strongly with one water molecule but, at the relatively high simulation temperature, interwater interactions are not strong enough to keep the cluster together. By pushing the system to react, on the other hand, the proton is forced toward the ethylbenzene and a transient hydronium ion is formed when water act as proton transfer bridge (Figure k). This hydronium ion is a strong H-bond donor and can be stabilized by the partial solvation offered by the remaining two water molecules, which remain tightly bounded over most of the simulation time.

For the 6 water molecules per unit cell the situation is quite different. In the reactant state there are enough water molecules to fully solvate the BAS as hydronium ion, thereby creating a reactive cluster mainly consisting of 3–5 water molecules (Figure f). By forcing the system toward the transition state, the hydronium ion is forced at the edge of the water cluster in order to interact with the ethylbenzene (Figure l). This reduces the possibility of being surrounded by the remaining water molecules and a smaller predominant cluster size of 3 is observed in the central region of the CV range.

While analyzing the size of the reactive water cluster gives valuable information about the chemistry of the system, we also adopted a second—more quantitative—approach to investigate the role of water on the reaction energetics. The 1-dimensional free energy profile was expanded in terms of two new collective variables, namely the distance between the ethylbenzene para carbon and the atom carrying the extra proton in the system and the minimum distance between the oxygen atoms in the first coordination sphere of the Al defect and, again, the atom carrying the extra proton. The latter is defined using a combination of the definitions proposed by Pérez de Alba Ortíz et al.[47] and Grifoni et al.[69] (more detailed information can be found in the Supporting Information, Section S5.2). The two-dimensional FESs for the various water loadings can be seen in Figure g–i. Starting from the lower loading of 1 water molecule per unit cell (Figure g), two separate paths can be seen going from the bottom (min Ozeo–proton distance = 0, proton on the zeolite) to the left (Cpara–proton distance = 0, proton on the ethylbenzene para carbon). The first one (p1 in the figure) corresponds to the anhydrous protonation, in which the BAS directly jumps on the ethylbenzene without water mediation (corresponding to the states with a cluster size of 0 around the transition state region in Figure d). The second path (p2), on the other hand, indicates that the water molecule can also act as proton transfer medium between the zeolite and the ethylbenzene (Figure j). The opening of this second possible protonation path which, as visible, is quite wide and similar in energy to the anhydrous protonation, is the reason for the large decrease in forward barrier going from 0 to 1 water molecule per unit cell.

Going toward higher water loadings, two things are visible. First, the anhydrous protonation path becomes less prominent as, with more water molecules in the catalyst unit cell, it is less and less likely for all of them to diffuse away from the active site, especially in the reactant region (compare with Figure e,f). Second, the number of states where the proton is solvated by the water as hydronium ion increases drastically, and in the case of 6 water molecules per unit cell, the proton can diffuse quite far from the active site while jumping from water to water. The creation of more and more stable states in the reactant region–due to the large mobility of the proton at high water loadings, is responsible for its increased stabilization and, ultimately, for the progressive increase of the forward reaction barrier. The transition state region, on the other hand, is qualitatively independent of the water loading. There, the hydronium ion must move toward the edge of the water cluster to interact with ethylbenzene (Figure k,l). This limits the possible interactions between the other water molecules and the hydronium ion and only 2–3 of them remain close by, as highlighted by the lack of cluster sizes greater than 4 in the transition state region of Figure f. Apart from the possibility of forming a few extra hydrogen bonds with the hydronium ion, the presence of extra water molecules in the catalyst does not significantly change the main features of the protonation transition state. Therefore, the changes in the reaction barrier can be mostly attributed to the stabilization of the reactant state deriving from the solvation of the hydronium ion when the number of adsorbed water molecules goes beyond one.

According to our simulations, water is expected to play a key role in modulating protonation reaction kinetics in zeolites. When looking at the forward reaction rate constant as a function of the water content (Figure ), an increase of more than 1 order of magnitude is expected when one water molecule is introduced in the reaction environment. The rate then progressively decreases with increasing amounts of water as the solvation of the BAS stabilizes the reactant state. An interesting comparison can be made with the rate of H/D exchange for benzene, which should proceed through a Wheland-like transition state. Chen et al. found through NMR spectroscopy measurements that the presence of trace amounts of water can increase the speed of benzene H/D exchange at room temperature, but only for the relatively low Si/Al ratio of 15 and not a higher ones of 40.[70] At higher—but still substoichiometric—water concentrations the exchange rate was found to be slower in all cases. Our computational model is certainly more in line with the Si/Al = 40 case consisting of isolated active sites, where a slowing down of the rate was experimentally observed for any amount of adsorbed water. It must be pointed out, however, that water is known to not uniformly distribute on the zeolite active sites and rather form heterogeneous clusters where some of the BAS are solvated while other basically anhydrous.[71] This means that the boost in activity that calculations suggest for 1 water molecule per BAS is unlikely to be measurable experimentally and rates similar or lower to the anhydrous case are going to be observed instead. At lower Si/Al ratios, cooperative effects between proximal BAS could significantly change the chemistry of the system and such effect is certainly worth of future investigation.

Figure 6

Forward kinetic constant (kf) and equilibrium constant (Keq) for the para protonation of ethylbenzene as a function of the number of water molecules in the zeolite unit cell. Dotted lines are present to guide the reader’s eye only.

Since, as shown before, the reverse phenomenological barrier remains quite similar for all water loadings, the equilibrium constant of the protonation reaction follows a similar trend as the forward kinetic constant (Figure ). In practice, water can speed up the protonation while interacting with the BAS, but then diffuses away once the proton has been transferred, thereby not changing the deprotonation kinetics. This reflects in an increase of the equilibrium constant from 6 × 10–6 in the anhydrous case to 8 × 10–5 in the presence of one water molecule. With 3 and 6 water molecules in the unit cell, the equilibrium constant decreases back to 7 × 10–6 and 4 × 10–6, respectively. Even with the most favorable conditions, the equilibrium constant remains always in strong favor of neutral ethylbenzene, suggesting that the Wheland complex will exist in traces inside the catalyst. The Wheland intermediate is thus suggested to exist as an elusive intermediate; however, it remains an open question in how it is sensitive to UV–vis in such low amounts. The predicted lifetimes found in our study are moreover an upper bound for its actual concentration, as pure GGA functionals, like the revPBE-D3 employed here, tend to overstabilize charged species.[72,73] As final remark, it is also possible that other types of active sites, for instance Extra Framework ALuminum (EFAL) species, could significantly increase the BAS acidity through cooperative effects[74] and, with it, the concentration of Wheland intermediates.

Conclusions

In this work, we report for the first time a complete mechanistic investigation of the benzene ethylation reaction with both ethene and ethanol using enhanced sampling molecular dynamics techniques to capture realistic operating conditions. Moreover, we thoroughly investigated the effect of various water loadings on the formation of Wheland complexes in the catalyst. At the chosen reaction conditions (623 K, 1 atm) the ethylation reaction is shown to proceed with analogous rates both through the stepwise mechanism (with the intermediate formation of a surface ethoxide species) and the concerted mechanism. The latter is, however, made unlikely by the restrain that both benzene and the ethylating agents must find themselves in proximity of the active site at the same time. In both cases, a transient Wheland complex is formed as reaction intermediate, which is however rather short living.

To understand the role of water on the protonation kinetics, we considered the more favorable para protonation of the ethylbenzene product and we analyzed it in the presence of 0, 1, 3, and 6 water molecules per zeolite unit cell. We showed that water can actively act as proton transferring agent, lowering the activation energy for the protonation reaction and increasing the rate of about 1 order of magnitude when 1 molecule per BAS is considered. At higher coverages, the BAS is partially or fully solvated by the water cluster. Such solvation strongly stabilizes the reactant state and is associated with a decrease in the protonation rate, which becomes basically as low as the anhydrous case with 6 water molecules per unit cell. The insights presented here provide further mechanistic details on the role and effect of water as a proton-transfer agent in zeolite-catalyzed reactions.