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.

. Schematic picture showing the passages from a 2-dimensional to a 1-dimensional free energy surface when one of the partial collective variables (CN1-CN2 in this case) varies twice as the other (CN3). Firstly, the initial 2D profile (a) is squeezed along the horizontal axis by multiplying by 0.5. In this way, both the partial CVs vary of one unit (b). Then, the profile is projected on the sum of the two partial CVs (red line in scheme b) giving the final 1D profile (c). The values of the (partial) CVs in the figure have been chosen to mimic the ones actually encountered in our simulations, but in principle any 2D CV space in which one CV varies twice as the other can be projected in 1D with a similar procedure. Table S1. Summary of the collective variables adopted for the simulation of all the transition states of benzene ethylation in the US simulations (see also Figure S4). The CNs used to define the CV are also reported, with their relative r 0 value (Å).

S2.2. Umbrella sampling simulations details
The approximate values that the CV and the CNs are supposed to assume in the reactant (R) and product (P) state of each reaction are also shown (according to the reaction arrows of Figure 1 in the main manuscript). Each CN can include the oxygen atoms in the first-coordination sphere of the Al defect in the zeolite framework (Ozeo), the aliphatic α-carbon of ethanol and SES or both the aliphatic carbon atoms of ethene (Cal), the oxygen of the hydroxyl group of ethanol (Oal) the six aromatic carbons of benzene (Car) or their respective hydrogens (Har) and the other non-aromatic hydrogen atoms (Hal).  Figure S4. CVs used in the umbrella sampling simulations for the formation of a surface ethoxide species from ethene (TS0 = ) or ethanol (TS0 OH ), for the direct alkylation of benzene with ethanol (TS2 OH ), for the protonation of ethylbenzene on a specific position (TS3, ipso and para carbon are considered in this work; in the para case, also the proton already linked to the carbon is considered in the CN as it becomes indistinguishable from the BAS once the protonation has occurred) and for the para protonation of ethylbenzene in the presence of water (TS3p nw with n≠0, 3 water molecules are shown in the figure as an example). TS0 = required some extra 2-dimensional umbrellas to improve the sampling of the transition state region, therefore two separate CVs are reported (see Section S2.3 for more information). For the nomenclature of the group of atoms see also the caption of Table S1.

TS
Each umbrella sampling simulation consists of multiple umbrella windows, distributed between reactants and products along the collective variable. In each window, a quadratic bias potential in the form V i (CV) = κ i 2 (CV − CV 0,i ) 2 is applied, being CV 0 the center of the umbrella and κ its force constant.
These two parameters are reported for each of the investigated transition states in the following tables. In general, low κ values were chosen to sample flat regions of the FES (the minima, in particular) while higher values were chosen for region with a large or unknown gradient of the FES.

S2.3. 2-D expansion of the US FESs.
While the linear combinations of coordination numbers used as collective variables allows to drastically reduce the computational time, making the reaction profile one-dimensional, the risk of poorly exploring important regions of the higher-dimensional FESfunction of the separate CNsincreases with the number of degrees of freedom that are projected on the linear combination. Luckily, by applying some simple statistical analysis, it is possible to expand a one-dimensional profile as a function of 2 or more new collective variables. 1 In formula: Where k B is Boltzmann's constant, T the temperature and F the free energy. This formula allows to expand a 1-dimensional FES function of an initial collective variable (F(CV)) in a new 2-dimensional FES function of two new collective variables (F(CV 1 , CV 2 )). To do so, the conditional probability p(CV 1 , CV 2 |CV) is used, which encodes the likelihood of observing the values CV 1 and CV 2 when, during the simulation, the value CV is visited. In practice, the 2-dimensional free energy F(CV 1 , CV 2 ) will be low when CV 1 and CV 2 are very likely to occur (high p(CV 1 , CV 2 |CV)) in correspondence of a CV value of low free energy (high exp(− F(CV) k B T ⁄ )) and vice versa. This formula was used to expand all the US FESs in two dimensions and check that the path connecting reactant and products was well-sampled (see Section S3.1) but also to convert the FES of TS3p nw as a function of two new collective variables to track the proton position (see Section S5.2). While in most of the simulations the whole reaction path resulted to be well-sampled, TS0 = presented obvious regions around the transition states that were insufficiently visited by the 1D umbrellas, when the free energy was expanded in function of x=0. 5 Figure S3b. Projecting the free energy on the sum of the two partial collective variables (x+y, as done in Figure S3c) can also be seen as an initial 45° rotation of the FES together with some stretching, creating a new 2D space defined by the horizontal axis x+y and its orthogonal direction x-y ( Figure S5a), followed by a projection on the horizontal axis. In this intermediate 2-dimensional space, the 1D umbrellas acting on x+y can be seen as 2D, with a force constant of zero in the orthogonal x-y direction ( Figure S5b). This allows to introduce proper 2-dimensional umbrellas, with force constants different from zero in both the x+y and x-y directions, to sample difficult regions of the FES and subsequently combining them with the original 1-dimensional umbrellas with 2-dimensional WHAM analysis. The obtained 2-dimensional FES can then be projected on the horizontal axis to retrieve a final 1-dimensional free energy profile, using the simple formula:

[CN(Cal;Hal)-CN(Ozeo;Hal)] and y=CN(Cal;Ozeo). Such two dimensional FES would correspond to the schematic one shown in
. When a 2-dimensional free energy surface is projected on the sum of the two collective variables of which it is function (see also Figure S3), the process can also be visualized as a 45° degrees rotation of the space (together with some stretching, as this is not a rigid rotation) followed by a projection on the new horizontal axis, x+y (a). In this new space, the vertical direction will consist of the difference between the two CVs instead of their sum (x-y). 1-dimensional umbrellas acting on the x+y direction can also be seen as 2-dimensional in this new space, with a force constant of 0 in the x-y direction (b) and can then be combined with 'regular' 2D umbrellas (force constant ≠0 in both directions) using a 2-dimensional WHAM analysis.

S2.4. Walls applied during the simulations
By biasing along the two CVs chosen to sample the FES of TS 1 and TS 2 = with MTD, the system could also evolve towards an undesired (meta)stable state while approaching the high energies of the transition state. Moreover, the diffusion of reactant/products also has to be limited in order to make barrier recrossings more likely to occur and improve the FES convergence. For this reason, a series of quadratic walls was applied to specific collective variables in the MTD simulations. A quadratic wall has the form V(CV) = 0.5 • κ(CV − CV 0 ) 2 if CV > CV 0 and 0 elsewhere. A list of the walls applied in the MTD simulations, together with their parameters and purposes, is reported in Table S10. Table S10. Summary of the quadratic walls adopted during the multiple-walkers well-tempered metadynamics simulations, together with their parameters and their purpose. The labelling of the atom species is the same as in Table S1.
The same applies for the umbrella sampling simulations. A complete list of all upper walls (i.e. acting when CV > CV 0 ) adopted in the US simulations, together with their parameters and purposes, is reported in Table S11. For TS1, also a lower wall (i.e. acting when CV < CV 0 ) was needed to prevent unwanted processes to occur. For the sake of clarity, it is reported separately in Table S12.  Few comments can be made on the walls selection; firstly, they can be divided into walls used to improve barrier recrossing and walls used to prevent unwanted side reactions. Starting from the formers, their presence will have an impact on the reaction FES. Indeed, they will drastically reduce the available configurational space for the adsorbed molecules, possibly causing a significant entropic penalty to the free energy (in their absence, the diffusing molecule could in principle explore the whole available volume in the zeolite unit cell). Being more interested in an estimate of the intrinsic reaction kinetics than in a full analysis of the transport properties in the material, we decided to set similar cutoff distances for similar reactions which should allow for a sufficient mobility of the species around the active site while allowing a fast convergence of the FES estimate. This approach is commonly used when studying zeolite-catalyzed reactions with enhanced sampling techniques (see refs. 2-5 as an example). The second type of walls, used to prevent side reactions, should have a minimal impact on the FES. The type of wall and the related parameters were mostly chosen based on chemical intuition and on a trialand-error procedure. For instance, TS3i required to limit the interactions between the reacting proton and all the other aromatic carbons, while for TS3p only walling the meta carbon atoms proved to be enough to avoid protonation of the spectating aromatic carbons (being the para protonation inherently easier than the ipso). The bias due to the walls presence was carefully monitored in all simulations, to ensure that they were hit only seldomly whenever the system was trying to go through an undesired reaction path (because of the random statistical motion of the atoms) and not in a continuous fashion, symptom of an improperly placed wall that would significantly impact the final free energy estimate.

S2.5. Error estimation
As stated in the main manuscript, the kinetic constant of a reaction can be estimated from an umbrella sampling simulation using classical transition state theory with the formula: 3 To obtain a final estimate for the confidence interval associated with the kinetic constant, the error associated with the FES estimate (F(CV)) and with the rate prefactor (A = √1 2πβ ⁄ 〈|∇ ⃗ ⃗ x CV|〉 CV ≠ ) were evaluated separately. An accurate estimate of the reaction free energy surface and the relative confidence interval can be obtained only by fully decorrelating the CV time series obtained from the US simulations. To achieve this, the Pymbar library 6,7 was used, with which we also automatically identified and discarded the simulation equilibration time. For the 2-dimensional umbrellas of TS0 = , the highest statistical inefficiency among the two computed separately for each CV was used to subsample the trajectories. In some umbrellas, very large statistical inefficiencies were computed which led to discard more than 95% of the data. While this approach necessarily increases the uncertainty of the free energy surface, it is also the only way to obtain a reliable and unbiased estimate of it. In order to adequately estimate error bars on the free energy surface obtained from WHAM analysis, we start by interpreting the WHAM equations as a maximum likelihood estimate (MLE). 8 Using the asymptotic normality of the MLE, one can construct the Fisher information matrix and use it to derive the covariance matrix associated with the MLE. 9 However, since the WHAM equations represent an MLE associated with various simulations, we start by computing the Fisher matrix ̿ for each simulation (consisting of simulation steps) separately and combine them afterwards into the global covariance matrix ̿: Note that this formula can be interpreted as simply adding up all available information. The square root of the diagonal elements of the covariance matrix then allows to compute the standard deviation on the probability (and hence free energy) of each bin in the free energy profile. By taking an error bar as two times the standard deviation, we compute a 2-sigma error bar corresponding to a 95% confidence interval (assuming normal distribution of the estimated free energy values, which is indeed the case for large samples sizes according to the asymptotic normality property of the MLE). Since the rate prefactor is evaluated only on a small interval around the transition state value (0.02 in CV units), the system will cross it only seldomly during the simulation, even in the umbrella on top of the TS. For this reason, the data obtained from the trajectories are much less correlated. The confidence interval on the rate prefactor was then estimate using the block averaging method. 10 Both the MLE and block averaging methods are implemented in the ThermoLIB library, which was used for the data analysis of all reaction steps. 11 To combine the two error estimates together, F(CV) and A were assumed to be normally distributed around their average value, with standard deviation computed from the MLE or block averaging methods, respectively. k TST was then recomputed 100,000 times using a Monte Carlo (MC) procedure to pick random values of both F(CV) and A from the respective normal distribution. It must be pointed out that, by selecting random values for F(CV), the position of the transition state (intended as the point with highest free energy between reagents and products) can also change. However, computing the rate prefactor for every possible TS value that can occur in the MC simulation would be extremely time consuming and practically unfeasible. We then decided to fix the average value and confidence interval of A to the ones obtained from the original FES. To validate this approach, TS0 OH was chosen as a case study (the relative FES can be seen in Figure S8). Firstly, we analyzed in how far the transition state position (intended as the point on the FES with higher free energy between reagents and products) can change during the MC run, since each point in the FES is randomly shifted up or down based on the confidence interval (bar plot in Figure S6, left). All values counted for more than 10 3 times lie within 0.1 CV units from the original TS and the maximal and minimal observed variation in the TS position are contained within 0.4 CV units. We then computed the rate prefactor on various equally spaced points within this interval (purple dots in Figure S6, left). While at the original TS a value of ~3·10 12 s -1 is computed, at the extremes of the intervals a decrease to ~2.5·10 12 s -1 can be seen. Such difference would correspond to a decrease in the kinetic constant of about k TST (1) k TST Once the rate constant and the related phenomenological barrier (both forward and backward) were computed for each MC step, we considered the average value of the barrier as final estimate with an associated 95% confidence interval given by twice the standard deviation. While this is possible since the values of the phenomenological barrier are normally distributed, the rate constant presents a lognormal distribution (Figure S6, right). For this reason, we computed the confidence interval as exp(〈log(k)〉) /exp(2 • σ(log(k))), exp(〈log(k)〉) • exp(2 • σ(log(k))).

S3. Benzene alkylation mechanism
S3.1. 2D US FESs As previously mentioned (Section S2.3), the US 1-dimensional FESs were expanded in 2 dimensions representing the CNs part of the linear combination defining the 1D CV. For CVs with format CV = CN(α; β) − CN(β; γ), the 2D FES is reported as function of CN(α; β) and CN(β; γ) while for CVs with format CV = CN(α; β) + 1 2 [CN(β; γ) − CN(γ; δ)] as function of CN(α; β) and CN(β; γ) − CN(γ; δ), but also as function of CN(β; γ) and CN(γ; δ) to indeed show that the sampling is adequate along every separate coordination number. In the following figures the final 1-dimensional free energy profile with associated error bars and the expanded two-dimensional free energy surface(s) for the various steps of the benzene ethylation reaction are shown. In the case of TS0 = , because of the additional 2-dimensional umbrellas, the final 2D FES is shown together with its projection along the horizontal axis.

S5.1. Cluster size involved in the reaction
When simulating the para protonation of ethylbenzene in the presence of water, the n (1,3 or 6) water molecules were manually placed in proximity of the active site at the beginning of the simulations. No restraints were imposed on their positions and therefore they were able to freely diffuse in the unit cell.
To understand whether water remains close to the active site, possibly playing a role in the protonation reaction, the size of the water cluster in proximity of the Al defect was monitored in the simulations. For each umbrella, 1 frame every 5fs was considered and, for each of such frames, a cutoff radius of 5Å around the 4 O atoms in the first coordination sphere of the Al defect (Ozeo) was set. Then, the presence of at least one water oxygen (Ow) within this cutoff radius was assessed. If no Ows were found, the cluster size was set to 0 and the following frame was considered. Otherwise, the cluster size was set to 1. From the Ow that is closest to Ozeo, a second cutoff radius of 4Å was defined and all other Ow in such cutoff counted and their number added to the cluster size. Then, the newly counted Ows became the center of a new cutoff radius to look for not yet included Ows and so on, until no further Ow were found ( Figure S19). For every water loading, the results of all umbrellas were joined and the percentage of occurrence of the various cluster sizes as function of the reaction CV reported (Figure 5d-f in the main manuscript). Figure S19. Schematic depiction of the H-ZSM-5 channel intersection, as seen along the sinusoidal channel. The approximate location of ethylbenzene during the simulations of TS3p nw is also shown. Three water molecules are randomly drawn in the intersection to illustrate the definition of water cluster size in proximity of the active site. Firsly, water 1 is found to be the closest to Ozeo and within the 5Å cutoff radius (blue region in the figure). Then, it is found that water oxygen 2 lies within 4Å from water oxygen 1. Finally, no further Ow are found within 4Å from water 1 or 2 (green circles in the figure). A water cluster size of 2 would then be associated with a simulation snapshot resembling this schematic picture (color code: H white, C grey, O red, Al pink and Si yellow).

S5.2. Expansion of the reaction profile in two dimensions
To quantitatively describe the energetic effect of the water-mediated proton hopping in the protonation reaction at various loadings, the free energy profile was expanded in term of two new collective variables, defining the proton-Ozeo and proton-Cpara distances. To define such fictitious proton atom, a combination of the collective variables recently proposed by Pérez de Alba Ortíz et al. 18 and Grifoni et al. 19 was used. Firstly, each proton in the system that could be involved in the protonation reaction (the BAS, all the water protons and the proton bounded to the ethylbenzene para carbon) is attributed to its closer proton acceptor (the water oxygens, Ow, the ethylbenzene para carbon, Cpara, and the oxygens in the first coordination spere of the Al, Ozeo) using the formula: In which i and m index the proton acceptors (A) and j runs over the protons. i is the position of atom i and λ is a parameter chosen to be 6 in this case. Thanks to this formulation, W i roughly equals the number of protons that are closer to atom i than to any other proton acceptor in the system. Indeed, the argument of the main summation tends to 1 if proton j is closer to atom i than to any other proton acceptor and quickly drops to zero otherwise. If a proton is exactly in between two proton acceptors, it will contribute to W i by 0.5 for both. For each i th proton acceptor an offset value W 0,i can be defined, from which the fictitious proton position can be computed: By using γ = 14, this weighted average will ensure that the position of the fictitious proton atom P corresponds with the position of the proton acceptor whose number of closer protons (W i ) is greater than its offset value W 0,i , i.e. the atom that carries the extra proton. Here, we selected an offset value of 0 for the four Ozeo, of 2 for all Ow of 1 for Cpara. Then, as final collective variables, we considered the minimum distance between the proton and Ozeo (CV 1,P = min(| P − O zeo |)) and the distance between the proton and Cpara (CV 2,P = | P − C para |). The former will be equal to zero when the proton is attached to the zeolite framework and greater than zero otherwise, while the latter will be equal to zero when the ethylbenzene has been protonated and greater than zero otherwise ( Figure S20). Figure S20. Schematic depiction of the H-ZSM-5 channel intersection, as seen along the sinusoidal channel. The adsorbed ethylbenzene and some water molecules (3, in this example) are also included. The variation of the two collective variables used to monitor the proton position when different atoms are protonated is shown. a. When the extra proton is located on the zeolite framework, CV1,P=0 (the location of the fictitious proton atom corresponds with the Ozeo on which the BAS is located) and CV2,P>0. b. When the proton is solvated by the water forming an hydronium ion, both CV1,P and CV2,P are greater than 0. c. When the ethylbenzene has been protonated, CV1,P>0 and CV2,P=0 (the location of the fictitious proton atom corresponds with Cpara). Color code: H white, C grey, O red, Al pink and Si yellow.