How Chain Length and Branching Influence the Alkene Cracking Reactivity on H-ZSM-5.

Catalytic alkene cracking on H-ZSM-5 involves a complex reaction network with many possible reaction routes and often elusive intermediates. Herein, advanced molecular dynamics simulations at 773 K, a typical cracking temperature, are performed to clarify the nature of the intermediates and to elucidate dominant cracking pathways at operating conditions. A series of C4-C8 alkene intermediates are investigated to evaluate the influence of chain length and degree of branching on their stability. Our simulations reveal that linear, secondary carbenium ions are relatively unstable, although their lifetime increases with carbon number. Tertiary carbenium ions, on the other hand, are shown to be very stable, irrespective of the chain length. Highly branched carbenium ions, though, tend to rapidly rearrange into more stable cationic species, either via cracking or isomerization reactions. Dominant cracking pathways were determined by combining these insights on carbenium ion stability with intrinsic free energy barriers for various octene β-scission reactions, determined via umbrella sampling simulations at operating temperature (773 K). Cracking modes A (3° → 3°) and B2 (3° → 2°) are expected to be dominant at operating conditions, whereas modes B1 (2° → 3°), C (2° → 2°), D2 (2° → 1°), and E2 (3° → 1°) are expected to be less important. All β-scission modes in which a transition state with primary carbocation character is involved have high intrinsic free energy barriers. Reactions starting from secondary carbenium ions will contribute less as these intermediates are short living at the high cracking temperature. Our results show the importance of simulations at operating conditions to properly evaluate the carbenium ion stability for β-scission reactions and to assess the mobility of all species in the pores of the zeolite.

Introduction

Zeolite-catalyzed alkene cracking is omnipresent in many petrochemical processes for the production of fuels and light olefins.[1−4] Recently, the increase in shale gas recovery and interest in renewable feedstocks have reinforced the search for new technologies to fulfill the rising propene demand. As a result, propene on-purpose processes have become an economically viable alternative for the conventional steam cracking technology.[5−9] In this context, catalytic alkene cracking has gained importance (e.g., to increase light olefin yields in product streams with a high content of less valuable C4–C8 alkenes).[9,10] Despite the wide application of alkene cracking, the nature of the reactive intermediates remains unclear.[3,11,12] Experimentally tracking individual reactions at the high operating temperatures (773–833 K) is a difficult task due to the occurrence of elusive intermediates and many side reactions.[13−18] Nowadays, molecular modeling techniques have matured to a level where zeolite catalyzed reactions can be modeled at true operating conditions, using advanced sampling techniques which allow to reconstruct the free energy surface.[18]

Herein, we unravel the stability and reactivity of typical alkene cracking intermediates at operating conditions with first-principle-based molecular dynamics methods. We show that static methods relying on a limited number of points on the potential energy surface inadequately account for the conformational mobility and temperature effects, hence giving a limited view on the zeolite catalyzed transformations under study. For alkane cracking, some molecular dynamics-based modeling studies clearly showed the importance of accounting for true operating conditions to unravel the nature of reactive intermediates and predict accurate product selectivities.[16] Buc̆ko et al. studied propane cracking in H–CHA and found that alkanes move rather freely in the zeolite pores at elevated temperature.[19] Later, this was also confirmed by Jiang et al., who studied the effect of temperature on the adsorption of short alkanes in SSZ-13.[20] Bell and co-workers showed for n-pentane and n-hexane cracking in H-ZSM-5 that more realistic product selectivities were obtained by means of innovative quasiclassical trajectory simulations.[21,22]

Upon alkene adsorption at the Brønsted acid site (BAS) of the catalyst, four different intermediates can be identified, as depicted in Scheme . A physisorbed van der Waals complex (1) is characterized by only dispersion interactions between the alkene and the zeolite walls. If the alkene double bond interacts with the acid proton, a physisorbed π-complex (2) is formed. Two chemisorbed intermediates can occur upon protonation of the double bond, either a carbenium ion (3) or a framework bound alkoxide (4). Especially the nature of the chemisorbed intermediates is strongly debated.[16,23−26] Despite the availability of many studies on the stability of carbenium ions and their impact on β-scission reactions, only a limited number of modeling studies are available on carbenium ion chemistry at operating conditions.[27,28]

Scheme 1

Four Possible Adsorption States of a 2-Alkene

1: van der Waals complex, 2: π-complex, 3: carbenium ion, 4: alkoxide.

The high reactivity of alkenes hampers the experimental characterization of the often short-lived intermediates.[16,18] Aside from a series of NMR[29−35] and FT-IR[26,36−41] studies, the nature of chemisorbed alkenes has mainly been studied theoretically. Early computational studies were carried out on small zeolite cluster models, consisting of only a few T atoms.[42−49] These spectroscopy experiments and cluster calculations both identified alkoxides as stable intermediates. Recently, Sarazen and Iglesia addressed the stability of bound species during alkene reactions on solid acids by a combined in situ infrared and density functional theory (DFT) study.[50,51] Carbenium ions, on the other hand, were suggested to act as transition-state species rather than cracking intermediates.[25,52,53] However, spectroscopy studies are typically carried out at relatively low temperatures. Since entropy effects disfavor the formation of framework-bound alkoxides,[27] different conclusions may arise at cracking temperatures. Furthermore, small cluster calculations lack a proper description of the zeolite confinement which may complicate alkoxide and facilitate carbenium ion formation.[16,17,54−57]

The existence of stable carbenium ions was proposed in a number of computational studies on extended cluster models or periodically repeated unit cells.[24,27,28,58−61] Fang et al. concluded that the accommodation of the hydrocarbons inside the zeolite framework is a crucial factor in determining carbenium ion stability.[61] The accommodation of carbenium ions is in addition influenced by the length of the adsorbed olefin.[62,63] Upon comparison with small cluster models, it was indicated that steric constraints destabilize bulky alkoxides relative to carbenium ions. Boronat et al. also stated that alkoxide stability is sensitive to the framework type as well as the position and local environment of the acid site.[24,55] Next to the zeolite topology, the stability of the intermediates and the cracking activity are also influenced by the acid strength of the zeolite.[64−67] Chu et al. demonstrated that the reactivity of dimerization and cracking reactions can be significantly enhanced by increasing the acid strength. The sensitivity to acid strength is especially important for intermediates with ionic character.[64]

The adsorption states of isobutene have received considerable interest. Upon isobutene protonation, either a tertiary butyl carbenium or a tertiary butoxide is formed.[24,68−72] Tuma and Sauer were the first to include entropy effects to periodic PBE calculations on H-FER.[68] They assessed the importance of the entropy effect in comparing stability differences at finite temperatures. At temperatures of 120 K or higher, the tert-butyl carbenium ion was shown to be more stable than the tert-butoxide.[68,69] This observation was later confirmed for a set of different zeolites by Nguyen et al.[70] Dai et al. obtained evidence for the existence of a tert-butyl carbenium ion in a combined computational and spectroscopy study.[72]

Recently, we investigated the nature of linear and branched C4 and C5 intermediates in H-ZSM-5 at reaction conditions using a combination of static DFT calculations and molecular dynamics simulations.[27,28] For the linear species, alkoxides and π-complexes were shown to be stable intermediates at low temperature (323 K). At a higher temperature of 773 K, though, alkoxides are very unstable and become short-lived inside the zeolite pores. Most modeling studies performed so far presumed the existence of alkoxide species even at high cracking temperatures,[50,58,73−75] whereas our modeling efforts point toward the importance of carbenium ions at operating conditions for catalytic cracking.[27,28] For the branched species, the tert-butyl carbenium ion was identified as the most stable intermediate, at both low and high temperature. Tertiary alkoxides were shown to be unstable, even at low temperature.[27,28]

Alkene cracking occurs through a complex network of oligomerization, isomerization, H-transfer, and cracking reactions. β-Scission of carbenium ions is generally accepted as the olefin cracking mechanism on solid catalysts.[3,11,56,76−78] Upon β-scission, a free alkene and a smaller chemisorbed carbenium ion are formed. These primary products can also be subject to further alkylations and subsequent cracking.[12] Fast isomerizations prior to the cracking reactions allow the formation of several carbenium ion types (primary, secondary, and tertiary) with different relative stabilities.[12,56,76] β-Scission reactions are therefore often classified according to the carbocation type involved in the reactant and product state (vide infra).

Several experimental studies focused on characterizing and predicting the olefin cracking product distribution.[75,79−88] Standl et al. developed a single-event kinetic model capable of correctly estimating olefin cracking selectivities.[85,86,89] Rate coefficients for each β-scission mode were considered to be independent of the chain length.[89] Nevertheless, Buchanan et al.[12] showed that C5–C8 olefin cracking rates over H-ZSM-5 at 783 K increase with carbon number. For longer olefins, the number of available isomers becomes larger; hence, more energetically favorable β-scission modes, involving branched, tertiary carbocations, are accessible.[12] However, pore size restrictions in H-ZSM-5 may prevent the formation of highly branched isomers.[12,75−77,90−95]

In previous modeling studies,[45,58,73,74,96] alkene cracking pathways have mainly been investigated with static DFT calculations. For C6–C8 olefin cracking in H-ZSM-5, Bhan and co-workers found that activation energies are determined primarily by the substitution order of the transition state and to a lesser extent by the substitution order of the reactant.[74] The highest activation energies were observed for β-scission modes involving a transition to a primary carbocation, while the lowest are found for β-scissions with a transition to a tertiary carbocation. However, the static modeling approach, considers only a single adsorbate configuration on the potential energy surface, while the energy surface for alkene adsorption is quite flat and contains many nearly isoenergetic local minima.[28] Also, at high temperature, entropy contributions may be underestimated considerably.[16,19,97,98] Many studies a priori postulated stable alkoxide intermediates as cracking reactants[58,96,73,74] and deduced the governing cracking modes based on the stability of alkoxides and intrinsic β-scission barriers determined from static DFT calculations.

In this work, we aim at gaining insight into the alkene cracking mechanism in H-ZSM-5 at operating conditions by using advanced molecular dynamics simulations. We investigate how the stability of carbenium ions in the reactant, transition and product states affect crucial β-scission rates at 773 K, a typical cracking temperature. To obtain a proper understanding of the stability of the intermediates formed upon C4–C8 alkene adsorption at the Brønsted acid sites, a selection of linear, mono-, di-, and tribranched alkenes with varying chain length are investigated as depicted in Scheme . Primary carbenium ions are not retained due to their highly unstable nature.[12,74,85,86] Ab initio molecular dynamics simulations, which inherently account for entropy effects and adsorbate mobility, are carried out to determine the lifetime of the various intermediates. In a next step, this knowledge is used to study the reactivity of linear and branched C8 intermediates toward monomolecular alkene cracking. Starting from the previously identified stable intermediates we estimate the reaction kinetics at operating conditions for five specific β-scission modes (Scheme ). Herein, we adopt the nomenclature, introduced by Weitkamp et al.[99] to distinguish between the elementary β-scission modes. Our approach relies on enhanced sampling molecular dynamics techniques to reconstruct the free energy profile at cracking conditions. A critical comparative assessment is made with activation energies at 0 K, predicted by static DFT calculations. Overall, the study yields new insights into the nature and reactivity of plausible intermediates and the impact of carbenium ion stability on the cracking pathways at operating temperatures.

Scheme 2

Adsorbed π-Complex and Carbenium Ion Intermediates of the Linear, Monobranched, and Multibranched C4–C8 Alkenes

Scheme 3

Five Possible β-Scission Modes for Monomolecular C8 Alkene Cracking Pathways

Computational Methods

Zeolite Model

All calculations are performed on a fully periodic model of H-ZSM-5, containing 289 atoms and a single acid site per unit cell. This way, the confinement of the zeolite pores can be properly accounted for. H-ZSM-5 exhibits the MFI topology, characterized by perpendicularly intersecting straight and sinusoidal 10-ring channels, forming a three-dimensional medium-sized pore system. A Brønsted acid site is created by substituting a Si atom at the T12 position by an Al atom and a charge-compensating proton. The active site is located at the intersection of the straight and sinusoidal channel, allowing maximal space for the adsorbed guest molecules. The time-averaged dimensions of the orthorhombic unit cell are obtained from a preliminary 5 ps molecular dynamics run in the NpT ensemble at 773 K and 1 bar (cf. Table S1).

Static Calculations

To evaluate activation barriers and reaction energies for a set of β-scission reactions, static DFT calculations are employed using the Vienna Ab Initio Simulation Package (VASP 5.4.1).[100−103] First, transition states for the β-scission reactions were searched by the improved dimer method.[104] To localize the reactant and product states corresponding to the identified transition state, the atoms were slightly displaced along the imaginary normal mode corresponding to the transition state motion. Then, the reactant and product ions were relaxed into the ground state with the conjugate gradient method. For all calculations, the Γ-point approximation for the sampling of the Brillouin zone was applied. Each stationary state has been optimized using the gradient corrected PBE functional[105] with additional Grimme D3 dispersion corrections[106] and the projected augmented wave method (PAW).[107,108] The plane wave kinetic energy cutoff was set to 600 eV while the self-consistent field (SCF) convergence criterion was fixed at 10–5 eV. The unit cell volume was kept fixed during the geometry optimization to eliminate the effect of lattice strain and to ensure the same truncated basis set is used. Energies are reported at the PBE-D3 level to allow proper comparison with the molecular dynamics simulations. The nature of the optimized local minima and transition states was verified by a normal-mode analysis. To obtain the normal modes, the partial Hessian vibrational analysis (PHVA)[109−111] was applied on the adsorbate and an 8T cluster of the framework, centered around the acid. Thermal corrections were estimated on the basis of the harmonic oscillator (HO) approximation. Ultimately, we determined the reaction enthalpies, entropies, and Gibbs free energies at the cracking temperature of 773 K using the in-house developed TAMkin software.[112]

Ab Initio Molecular Dynamics Simulations

Molecular dynamics simulations were performed to assess the nature and lifetime of alkene cracking intermediates. The dynamic behavior of some intermediates may depend on the reaction conditions such as temperature and pressure. In contrast to static calculations, framework flexibility and entropy effects can be inherently accounted for. All ab initio MD simulations have been carried out with the CP2K software package (CP2K 3.0),[113] using the revPBE[114] level of theory with additional Grimme D3 dispersion corrections,[106] a DZVP basis set and GTH pseudopotentials.[115] A combination of Gaussian and plane waves (GPW)[116,117] basis functions with a cutoff energy of 320 Ry are used. The self-consistent field (SCF) convergence criterion was set at 10–6 eV. The time step for integrating the equations of motion was set at 0.5 fs. The simulations were performed in the NVT ensemble at a temperature of 773 K, controlled by a chain of five Nosé–Hoover thermostats.[118,119]

After an equilibration run of 5 ps, a production run of 100 ps was performed to obtain a sufficient sampling of the phase space. When isomerizations or other rearrangements occurred in the equilibration run or at the beginning of the production run, the simulation was repeated thrice with reinitialized starting conditions. This will ensure that the observed transition is not merely a rare event, caused by the random choice of initial velocities.

The identification of the various possible intermediates that are sampled in the MD simulation was computed on the basis of a distance criterion. The adsorbate is classified as a physisorbed π-complex if the distances between the acid proton and the double bond carbon atoms are both smaller than 2.85 Å; otherwise, it is considered as a van der Waals complex. The sampled state is considered a chemisorbed carbenium ion if all distances between a hydrogen and an acid site oxygen (O) are larger than 1.25 Å. To track the position of the positive charge in the carbenium ion state, the trivalent carbon atom is identified by assigning each hydrogen to the closest carbon. An alkoxide is sampled if a covalent C–O bond with the framework smaller than 1.9 Å is formed.

Metadynamics Simulations

The probability of sampling rare events, like isomerization or cracking reactions in regular MD is typically very low. To enhance the sampling of these activated transitions, accelerated sampling techniques are required. Metadynamics simulations (MTD) are applied to assess the free energy difference between the elusive intermediates.[120,121] We applied the same procedure as described in ref (27). The sampling of specific reaction paths on the free energy surface requires the definition of a collective variable (CV), uniquely describing the reaction coordinate.[120−122] MTD simulations are run with the same settings as for the MD simulations, i.e., in the NVT ensemble at actual cracking conditions (773 K). At regular instants during the simulation, Gaussian hills are spawned along the collective variable which is defined based on coordination numbers (CN):in which the sum runs over two sets of atoms i and j with r the interatomic distance between atoms i and j and r0 a reference distance. The parameters nn and nd are set at 6 and 12, respectively. All simulations in this work use a single CV (1D MTD): CN(Ch-Hh) – CN(Hh-O) to describe the proton transfer from the zeolite (z) to the hydrocarbon (h), shown in Scheme . A reference distance of 1.25 Å is selected which is close to typical transition state bond distances for (de)protonation reactions.

Scheme 4

Collective Variable (Difference of Coordination Numbers) Used in the Metadynamics Simulations of the Protonation of a 2-Alkene into a 2-Alkyl or 3-Alkyl Carbenium Ion

To ensure a clear distinction between reactant, product, and transition state, isomerizations were prevented; hence, the same carbon atom would get protonated each time. To this end, the simulation is constrained to the particular area of interest by applying quadratic walls. Hills with a height of 2 kJ/mol and a width of 0.035 are added every 100 time steps. Upon each recrossing of the barrier, the hill height was reduced by half to improve the convergence of the free energy profile. The latter is reconstructed on the basis of the sum of all spawned Gaussian hills along the reaction coordinate. The simulations were continued until the addition of extra hills no longer changed the resulting free energy profile. Finally, activation free energy barriers ΔG‡ were computed as the difference between the maximal free energy in the transition state (TS) region and the minimal free energy in the reactant (R) or product (P) valley.

Additional technical details on the MTD simulations can be found in Section S2 of the SI.

Umbrella Sampling

In metadynamics simulations, the reaction profile is obtained through multiple crossings of the activated transition state, separating the reactant and product valley. Accordingly, this requires the existence of clearly defined, stable minima to allow the forward and backward reactions to be properly sampled. However, in the case of β-scission reactions the formed reaction products are often unstable and prone to rapid rearrangements. Since these reorganizations are not uniquely described by the selected collective variables, the sampling of the backward reaction may be hindered. Therefore, the more robust umbrella sampling (US) method[123,124] is used to determine the free energy profiles for the selected alkene cracking reactions. With this technique, the conformational space is subdivided into a number of windows along the reaction coordinate that are sampled individually during the free energy simulations. To the best of our knowledge, this is the first application of umbrella sampling to study zeolite-catalyzed reactions.

All umbrella sampling simulations have been performed with CP2K as MD engine, interfaced with the PLUMED module.[125] The protonation of a linear alkene into a carbenium ion is sampled using a coordination number between the acid site oxygen atoms and the hydrocarbon hydrogen atoms, CN(O-Hh) as collective variable (Scheme a). To properly describe the cracking reaction, a coordination number between all atoms involved in the β-scission; e.g., the carbon atoms in α and β position of the positive charge, CN(Cα-Cβ), is defined as collective variable (Scheme b). The total range of the collective variable, q, is divided into a number of windows. For each window, biased MD simulations are run in parallel. The initial configuration for each window is randomly selected from a moving bias potential simulation, describing the entire collective variable range of the reaction. To ensure a sufficient sampling along the entire reaction coordinate, harmonic bias potentials, centered around the equilibrium value q0, are employed.

Scheme 5

Collective Variables (Coordination Numbers) Used during the Umbrella Sampling Simulations of (a) the Protonation of 2-Alkenes and (b) the β-Scission of a Carbenium Ion

The free energy difference in each window can be obtained from the sampling distribution along the specified collective variable. The weighted histogram analysis method (WHAM),[126,127] implemented by the Grossfield Laboratory,[128] is applied to reconstruct the free energy profiles by combining the results of each window. An overview of the used collective variables, the umbrella positions and bias strengths in the US simulations can be found in Table S2 of the SI.

Results and Discussion

Influence of Chain Length on Carbenium Ion Stability

The stability of adsorbed alkene intermediates is—next to temperature and branching—influenced by the chain length of the adsorbate. Longer alkyl chains are expected to adsorb more strongly due to the increased contribution of dispersion interactions. Entropy is also expected to become more important for longer chains.[50] To assess the impact of the chain length effect, we investigated the dynamic behavior of a series of linear and branched C4–C8 alkene intermediates (see Scheme ) with MD simulations at 773 K. In our previous work, we showed that both linear and branched pentoxides are unstable inside the zeolite pores at elevated temperature. Instead, they immediately transform into carbenium ions.[27,28] Based on these observations, we can safely assume that longer alkoxides will not exist either. Therefore, we only performed MD simulations starting from the physisorbed alkene π-complex and the chemisorbed carbenium ion.

Linear Alkenes

For the linear alkenes, we considered the homologous series of 2-alkenes, ranging from 2-butene to 2-octene. In the course of the simulations, spontaneous alkene protonation and deprotonation reactions are occasionally observed. Since the transitions between neutral alkene and carbenium ion species are sampled in regular MD simulations, we can evaluate their relative stability by comparing the lifetime of these intermediates. Figure summarizes the resulting time fractions during which the alkene π-complex, the alkene van der Waals complex, and the carbenium ion are sampled in the various MD simulations. Starting from the carbenium ion configuration, two independent simulations (A and B) are conducted.

Figure 1

Fractions of the sampled intermediates during the 100 ps MD simulations on the linear C4–C8 species at 773 K, with either the π-complex or carbenium ion as starting configuration. Two independent carbenium ion simulations are denoted as A and B (*data taken from ref (27)).

Figure reveals some trends which are specific for the linear alkenes. First, the combined fractions of the physisorbed states (π-complex and van der Waals complex) range from 86 to 100% in every simulation, showing that the neutral alkene is significantly more stable than the carbenium ion for the entire C4–C8 range. The van der Waals complex is slightly favored over the π-complex, but frequent transitions (500–1000 per simulation) take place between the two physisorbed states. No clear trend with increasing chain length can be identified for the ratio of π-complex to van der Waals complex. Despite the additional stabilization from the π-H interaction in the π-complex, the more loosely bound van der Waals complex is entropically favored at high temperature. To study the mobility and the entropic contributions of the various species in the zeolite pores in depth, we postprocessed our MD data to determine for each adsorbate the translational part to the entropy in the zeolite pore system and to construct 2D scatter plots visualizing the mobility of the various species and 1D probability distributions along the x, y, and z directions. The y direction corresponds to the direction of the straight channel, whereas the direction of the sinusoidal channel lies in the xz plane. More details are given in Section S5 of the SI. Furthermore, for 2-butene and 2-octene, we disentangled the contributions from the π-complex and van der Waals complex. The resulting mobility plots for 2-butene are shown in Figure .

Figure 2

Scatter plot of the center of mass of (a) the 2-butene π-complex and (b) the 2-butene van der Waals complex inside the channel system of H-ZSM-5, projected on the xy plane. The purple sphere represents the location of the acid site.

It is immediately clear that the van der Waals complex has a much larger mobility compared with the π-complex and travels along relatively long paths in the straight channel of the zeolite during the ca. 70 ps in which the van der Waals complex is sampled. The entropic stabilization of the van der Waals complex due to translational movements in the straight channel is about 7 kJ/mol higher than the entropic stabilization of the π-complex. The translational entropy contributions for all species are summarized in Table S5. Similar mobility plots for 2-butene, 2-octene, isobutene, and 2,3,4-trimethyl-2-pentene are shown in Figure . The distinction between the various adsorbed states has not been made anymore, as the latter exist for the majority of time in their carbenium ion state (vide infra). The mobility of 2-butene and 2-octene are substantially different. 2-butene enters more deeply into the sinusoidal channel, whereas the longer 2-octene chain prevents the adsorbate to enter the sinusoidal channel efficiently in a simulation time span of 90 ps. Instead, 2-octene travels further along the straight channel of the zeolite, which yields an additional entropic stabilization of about 14 kJ/mol for the 2-octene van der Waals complex compared to the π-complex.

Figure 3

Scatter plot of the center of mass of (a) 2-butene, (b) 2-octene, (c) isobutene, and (d) 2,3,4-trimethyl-2-pentene inside the channel system of H-ZSM-5, projected on the xy plane. The purple sphere represents the location of the acid site.

Second, (de)protonation reactions occasionally occur during the course of the simulations. The probability of sampling a protonated alkene is small but not negligible at higher chain lengths. Approximate lifetimes for the various carbocations are shown in Figure S5. For longer alkenes, the maximum carbenium ion lifetime increases due to inductive and hyperconjugative stabilization which allow for a better accommodation of the positive charge. During the sampling of the cationic state, low-barrier hydride shifts are frequently observed. This effect also contributes to the enhanced stabilization of longer alkyl cations. While only two carbon atoms can bear the positive charge in a n-butyl carbenium ion, the charge can effectively be delocalized over six carbon atoms in an octyl carbenium ion. Isomerization occurs either via a 1,2-hydride shift, a 1,3-hydride shift, or a 1,5-hydride shift for the longer alkyl chains. The positive charge distribution among the carbon atoms of each alkyl carbenium ion is visualized in Figure . Owing to the short sampling times of the cationic state, the reported distributions do not represent the equilibrium. Nevertheless, they clearly show that the positive charge preferably resides on the central carbon atoms of the chain.

Figure 4

Delocalization of the positive charge (%) along the alkyl chain during the carbocation sampling of the 2-alkyl carbenium ion simulations (simulation A and B combined).

The transformation of a physisorbed complex into a carbenium ion is an activated process, even at higher temperatures. Cracking modes starting from linear carbenium ions (D2 mode, see Scheme ) would first have to overcome the barrier for carbenium ion formation. To quantify the free energy difference between the carbenium ion and the physisorbed alkene more accurately and to determine free energies of activation for their interconversion, metadynamics simulations are carried out for the protonation of 2-alkenes.

The uncertainty on free energy barriers from a single metadynamics simulation can be quite large.[27] Therefore, to obtain reliable free energy estimates, each simulation was repeated 5 times and the average and spread on the resulting activation barriers and reaction free energies are reported in Table . The 2-alkene protonation barrier (ΔG‡fwd) has a constant value of 50 kJ/mol approximately, regardless of the chain length. The free energy barrier for deprotonation (ΔG‡bwd) is systematically lower than the protonation barrier. The relatively low deprotonation barriers, combined with the overall higher stability of the alkenes are in line with the relatively small fraction of the time in which carbenium ions are sampled in the MD simulations. The free energy differences (ΔGr) corroborate the observation that the positive charge is better accommodated at central chain positions. Both the C4 and C5 carbenium ions are around 25 kJ/mol higher in free energy than the corresponding alkenes. For the C7 and C8 species, a clear distinction can be made between the 2-alkyl and 3-alkyl carbenium ion. While the 2-alkyl carbenium ions are about 21 kJ/mol higher in free energy than the corresponding 2-alkenes, the 3-alkyl carbenium ions are only 10 kJ/mol less stable than the 2-alkenes.

Table 1

Average Free Energy Activation Barriers and Reaction Free Energies at 773 K for the Protonation of the Linear C4–C8 Alkenes into the Corresponding Carbenium Ions from Metadynamics Simulationsa

		ΔG‡fwd (kJ/mol)	ΔG‡bwd (kJ/mol)	ΔGr (kJ/mol)
C4	2-butene → 2-butyl cation*	52 ± 4	23 ± 2	26 ± 5
C5	2-pentene → 2-pentyl cation*	49 ± 4	26 ± 2	23 ± 4
	2-pentene → 3-pentyl cation*	49 ± 4	22 ± 6	27 ± 5
C6	2-hexene → 2-hexyl cation	48 ± 3	30 ± 4	18 ± 1
	2-hexene → 3-hexyl cation	46 ± 3	30 ± 4	17 ± 7
C7	2-heptene → 2-heptyl cation	47 ± 3	23 ± 4	24 ± 2
	2-heptene → 3-heptyl cation	47 ± 3	38 ± 3	9 ± 4
C8	2-octene → 2-octyl cation	50 ± 3	26 ± 5	23 ± 7
	2-octene → 3-octyl cation	48 ± 3	38 ± 3	10 ± 3

The symbol * indicates data taken from ref (27).

In principle, free energy differences between the various adsorbed states can also be roughly deduced from regular MD simulations provided that sufficient transitions are sampled. The free energy differences between the linear alkenes and carbenium ions from the regular MD simulations are presented in Table S8 in SI. Although the carbenium ion sampling time is rather short, the ΔG estimates from regular MD and MTD simulations are in good agreement. For all C5 to C8 MD simulations with a sufficiently large carbenium ion sampling time, the alkene/carbenium ion free energy differences lie in the same range as those predicted by the MTD simulations.

Branched Alkenes

Branched alkenes contain tertiary carbon atoms which may significantly impact the stability order of the intermediates. In earlier work, we have shown that the tertiary carbenium formed upon isobutene and 2-methyl-2-butene adsorption is a stable intermediate at 773 K.[27] In this section, we evaluate the influence of chain length on C4–C8 alkenes with one or more methyl branches.

Figure shows the fractions of the states sampled in MD simulations at 773 K for the monobranched 2-methyl-2-alkenes (cf. Scheme ). Irrespective of chain length, carbenium ions from branched alkenes are clearly much more stable than from linear alkenes. In contrast to the linear alkenes, no clear trend with higher carbon number can be discerned. Throughout the largest part of the simulations, tertiary carbenium ions are sampled. Hydride shifts to secondary carbenium ions are rarely observed. The positive charge resides almost exclusively on the tertiary carbon atom, reflecting the higher stability of tertiary carbenium ions. This is also evidenced by the lifetime of these carbenium ions across the different simulations, depicted in Figure S6. The observed maximum lifetime of the tertiary carbenium ions is much larger than for any of the linear, secondary carbenium ions.

Figure 5

Fractions of the intermediates sampled during the 100 ps MD simulations on the monobranched C4–C8 species at 773 K, either with the π-complex or carbenium ion as starting configuration (*data taken from ref (27)).

All branched 2-methyl-2-alkenes are also stable, although the physisorbed alkene fraction is much smaller compared to linear alkenes. The π-complex comprises only a very small fraction, indicating that the van der Waals complex is by far the most stable physisorbed state. Clearly, the bulkier branched alkenes are sterically more hindered making it more difficult to approach the acid site to form the π-complex.

Similar observations can be made for C6–C8 alkene isomers with multiple side chains. Figure displays the fractions of the sampled intermediates in MD simulations on a set of di- and tribranched alkenes at 773 K. Carbenium ions are again the dominating intermediates due to the presence of tertiary carbon atoms. The smallest tribranched 2,3,4-trimethyl-2-pentyl species exists solely as carbenium ion, regardless of the starting configuration. Immediate protonation takes place at the beginning of the alkene simulation. Despite the presence of neighboring tertiary carbon atoms, almost no H-shift isomerizations are observed for the multiple branched chains. The barrier for isomerization between different tertiary carbenium ions appears to be much higher than for secondary carbenium ions. In contrast to the linear chains, for which isomerization partly provided stabilization, tertiary carbenium ions are stable, regardless of the positive charge position along the chain. All carbenium ions for the branched species remain positioned at the intersection of the straight and sinusoidal channel as is clear from the mobility plots for isobutene and 2,3,4-trimethyl-2-pentene, shown in Figure c,d respectively. Due to the cationic nature, these species remain relatively close to the positive acid site to avoid excessive charge separation.

Figure 6

Fractions of the sampled intermediates during the 100 ps MD simulations on the multiple branched C6–C8 species at 773 K, either with the π-complex or carbenium ion as starting configuration.

For C6+ alkenes, geminal methyl substituted intermediates, i.e., with a quaternary carbon atom, can also be formed. However, our MD simulations indicate that these species are extremely unstable. In all simulations of the homologous series of 2,2-dimethyl-3-alkyl secondary carbenium ions (Scheme ), immediate rearrangements are observed. The 2,2-dimethyl-3-butyl carbenium ion instantly isomerizes via a methyl shift into the 2,3-dimethyl-2-butyl carbenium ion. The larger 2,2-dimethyl-3-pentyl and 2,2-dimethyl-3-hexyl carbenium ions both stabilize through a rapid H-shift, followed by a β-scission reaction, forming the tert-butyl carbenium ion and propene/butene. Opposed to these extremely reactive secondary carbenium ions, the tertiary and highly branched 2,4,4-trimethyl-2-pentyl carbenium ion remains stable for the entire simulation time. However, previous studies have reported that the formation of species with geminal methyl substitutions or more than two methyl substitutions is prohibited by the pore dimensions.[76,85,93−95]

In conclusion, at typical cracking temperatures, linear alkenes are predominantly found as π-complex or van der Waals complex in the zeolite pores. Linear carbenium ions are short-living, although their lifetime increases with chain length, since the positive charge can be delocalized over more carbon atoms. Alkenes with one or multiple branches occur almost exclusively as carbenium ions with a relatively low mobility compared to linear alkenes.

Reactivity of C8 Cracking Intermediates

β-scission of carbenium ions is widely accepted as the predominant mechanism for alkene cracking in zeolites. Next to the stability of the cationic reactants, the stability of the product cations, formed in the transition state, is a key factor determining the reactivity.[74,75] From the MD analysis, physisorbed alkenes and tertiary carbenium ions were identified as the most stable species for linear and branched alkenes respectively. In this section, the reactivity of the governing C8 intermediates toward catalytic cracking of different β-scission modes is assessed. Three possible cracking modes for the highly stable tertiary carbenium ion reactants are considered, as shown in Scheme : a transition from a tertiary carbenium to a primary carbenium ion (mode E2), to a secondary carbenium ion (mode B2) or to a tertiary carbenium ion (mode A). For linear alkenes, only a secondary to primary carbenium ion transition (mode D2) can occur. Transitions from a secondary carbenium ion to a secondary carbenium ion (mode C) or to a tertiary carbenium ion (mode B1) are not explicitly considered as our MD simulations indicate that these cracking precursors tend to immediately isomerize to more stable tertiary carbenium ions. We therefore assume that isomerization occurs much faster than cracking and that the contribution of these β-scission modes will be minimal. For the B2-1 reaction, both static and molecular dynamics simulations are performed to evaluate the impact of true operating conditions on the cracking activity. For this particular reaction, this was possible since tertiary carbenium ions are still relatively stable at low temperature,[27] and hence, they can be localized with static calculations.

Cracking Barriers from Static Calculations

First, static DFT calculations were performed at the PBE-D3 level of theory to quantify the intrinsic activation barrier for the cracking of a tertiary 2,4-dimethyl-2-hexyl carbenium ion (reaction B2-1). This tertiary, dibranched cation was found to be a very stable intermediate, and hence, it can be considered as a probable cracking reactant. Upon β-scission, isobutene and a 2-butyl cation are formed. This reaction is an example of the B2 cracking mode in which a transition of a tertiary into a secondary carbenium ion takes place. To illustrate that a multitude of transition states exist and that the zeolite environment plays a crucial role in the stabilization of the transition state, four differently oriented transition state geometries around the BAS were localized, as shown in Figure . More information on the particular transition states and their geometrical characteristics is given in Section S10 of the SI.

Figure 7

Sinusoidal channel view of the static geometries for the four B2-1 transition state structures.

In the optimized transition states, the 2-butyl cation interacts with both the isobutene C=C bond and the nearest basic zeolite oxygen. In the product state, the 2-butyl cation immediately binds to the nearest framework oxygen, thus forming a 2-butoxide. Indeed, static calculations previously showed that 2-butoxide is substantially more stable than a 2-butyl carbenium ion.[27] However, at cracking temperatures, entropy contributions disfavor the formation of stable alkoxides, which is confirmed by the MD simulations on C4 species (Figure ), where the 2-butyl carbenium ion stabilizes by transforming into a π-complex.

The four different cracking reactions also show large energetic variations. Table lists the enthalpy and free energy barriers at 773 K. We notice large variations of almost 80 kJ/mol on the forward activation barriers. This trend originates mainly from an enthalpy effect, but is enhanced by the entropy contributions. Upon studying the same B2-1 reaction, Mazar et al. also found a large variation (82 kJ/mol) in activation energy between two distinct reaction configurations.[74]

Table 2

Enthaply, Entropy, and Free Energy Barriers at 773 K for the Four Conformations of 2,4-Dimethyl-2-hexyl Carbenium Ion Cracking (Reaction B2-1)

structure	ΔH‡fwd (kJ/mol)	ΔS‡fwd (J/mol/K)	ΔG‡fwd (kJ/mol)
1	151	–20	167
2	120	–15	132
3	104	–25	124
4	94	5	90

These results show that many different pathways connect the reactant carbenium ion and cracked products, characterized by largely different activation barriers. The same feature has been reported for propane cracking in chabazite by Buc̆ko et al.[19] Alkene cracking reactions taking place at high temperature have a complex nature, where a restricted number of configurations on the PES are insufficient to capture the complexity of the transformation. In static calculations, only a single configuration of the adsorbed transition state is considered. Furthermore, the 0 K energy surface may differ from the free energy surface at cracking temperatures.[21,28,98] The latter is evidenced by the 2-butoxide formation in the product state. Additionally, the harmonic oscillator (HO) approximation underestimates the mobility of the adsorbates and entropy of the stationary states.[129,130] Consequently, relying solely on static calculations seems insufficient to correctly describe the high-temperature cracking chemistry. Hereafter a more realistic study of the reactivity is performed based on molecular dynamics simulations at operating conditions.

Free-Energy Profiles for Cracking Reactions at Operating Conditions

To properly account for all possible reaction paths as well as finite temperature and entropy effects, the umbrella sampling technique is applied for determining cracking barriers. We first simulated both the β-scission of the 2,4-dimethyl-2-hexyl carbenium ion, reaction B2-1, and the 2,4-dimethyl-4-hexyl carbenium ion, reaction B2-2 (Scheme ). A broad range of transition states exists for these cracking reactions as indicated by the mobility plots shown in Figure a. The red and blue dots represent the position of the center of mass of the two C4 fragments. During the umbrella sampling simulation, a large variety in the relative orientation of the two C4 fragments in the transition state region is sampled. The most visited configurations are localized in the straight channel or in the intersection of the straight and sinusoidal channel and resemble on average mostly the most stable static transition states (3 and 4). Figure c displays a superposition of transition state structures (corresponding to a collective variable value of 0.25) to show the conformational freedom in the activated complex region. Also the orientation of the breaking C–C bond, indicated by the green bars, significantly fluctuates throughout the simulation.

Figure 8

(a) Scatter plot of the center of mass of the transition state window configurations for the B2-1 cracking reaction in the channel system of H-ZSM-5, projected on the xy plane. The dots connected by green bars indicate the static transition state configurations. (b) Scheme of the B2-1 transition state with the isobutene fragment in red, the n-butyl fragment in blue and the breaking C–C bond in green. (c) Superimposed snapshots of four different transition state window configurations. The purple sphere represents the location of the acid site.

The resulting free energy profiles are shown in Figure . For the two B2 type β-scissions, cracking barriers of 69 and 73 kJ/mol are obtained, respectively. These barriers are significantly lower than the lowest estimate predicted by static calculations (90 kJ/mol). At 773 K, the adsorbed species have a high mobility, allowing them to access many different configurations. On average, the sampled transition state configurations are mainly situated in the straight channel or at the intersection, but at longer distances from the framework compared to static calculations (see Figure S23). By considering just a single transition state and reactant configuration, static calculations clearly give an incomplete picture of the cracking barriers.

Figure 9

One-dimensional umbrella sampling free energy profiles at 773 K for the β-scission reactions B2-1 (a), B2-2 (b), E2-1 (c), A-1 (d).

Furthermore, the height of the cracking barrier seems almost independent of the cracking products (isobutene/2-butyl cation or 2-methyl-1-butene/2-propyl cation). These findings indicate that the activation barrier is mainly determined by the type of β-scission mode—in this case a tertiary to secondary carbenium ion transition—rather than the type of alkene products formed in the reaction. The observation that reactions of the same cracking mode have similar free energy barriers is supported by von Aretin et al.[131] They established a microkinetic model capable of correctly estimating alkene cracking product distributions with the assumption that cracking reactions of the same mode have the same rate coefficient.

Upon cracking, unstable secondary carbenium ion products are formed. As evidenced in earlier work,[27,28] the n-propyl and n-butyl cations are unstable and will immediately rearrange. For these cracking products, two different transitions can be identified. The cation can either deprotonate to the framework or transfer a proton to the branched alkene product, thus forming two neutral alkenes or an alkene and tertiary carbenium ion as stable products. The latter rearrangement seems to be preferred as it was spontaneously observed in the US simulations. Static calculations on the two product states also confirmed that the alkene/tert-butyl cation products are about 12 kJ/mol more stable than the alkene/alkene products. Once the cracking barrier is crossed, stabilizing isomerization reactions occur and the produced product fragments quickly diffuse away from each other, resulting in a sharp free energy drop (Figure a). Because of the large configurational space of the separate cracked fragments, the product state is insufficiently sampled to yield physically meaningful backward reaction barriers, hence these are not reported.

Second, we simulated the cracking of a 2-methyl-2-heptyl carbenium ion into isobutene and a 1-butyl carbenium ion (reaction E2-1). This cracking mode has a much higher activation barrier of 112 kJ/mol (Figure c) which can be explained by the formation of an unstable primary cation in the transition state. After cracking, the unstable 1-butyl carbenium ion product instantly rearranges via a rapid hydride shift to a more stable secondary 2-butyl carbenium ion. The latter again undergoes a hydride shift with the formation of n-butene and a tert-butyl carbenium ion. This cracking mode is expected to be less important due to the primary carbenium ion character of the transition state. However, monobranched cracking precursors will exhibit less steric repulsion with the framework compared to highly branched alkenes and will therefore be more probable to occur.

Next, we simulated the cracking of the tribranched 2,4,4-trimethyl-2-pentyl carbenium ion into isobutene and the tert-butyl carbenium ion (reaction A-1), the only possible A-mode cracking for C8 species. In Section 3.1.2, we showed that geminal substituted carbenium ions were more reactive than monobranched carbenium ions. Nevertheless, the 2,4,4-trimethyl-2-pentyl cation was still observed as a stable intermediate. Since two stable products are formed upon cracking, also no product rearrangements are observed. As confirmed from the MD simulations, both states are long-living and coexist inside the zeolite pores. For this cracking mode, a relatively low barrier of 53 kJ/mol (Figure d) is obtained which can be attributed to both the high reactivity of the cracking reactant and the formation of a stable tert-butyl cation in the transition state. However, the formation of highly branched cracking precursors may be prohibited due to steric constraints on ZSM-5.[85,93,95] Indeed, the umbrella sampling simulations show a limited mobility for the 2,4,4-trimethyl-2-pentyl cation, which is mainly situated at the channel intersection.

Finally, we consider the cracking of a linear C8 carbenium ion (D2 mode). Since secondary carbenium ions were shown to be relatively unstable cracking precursors with short lifetimes, the neutral alkene should rather be considered as cracking reactant. Therefore, two-dimensional umbrella sampling is performed with a first CV describing the protonation of 2-octene and a second CV describing the β-scission of the resulting 3-octyl carbenium ion (reaction D2-1). The free energy profile for this cracking pathway is shown in Figure . For the protonation, an activation barrier of 55 kJ/mol and a reaction free energy of 10 kJ/mol are obtained, which is in good agreement with our MTD simulations (Table ). For the cracking, a relatively high intrinsic barrier of 103 kJ/mol is obtained, which can again be explained by the combination of a relatively unstable reactant cation and an unstable product cation. The formation of a primary carbenium ion in the transition state is clearly unfavorable, yet unavoidable for linear chains. In the product state, the primary carbenium ion immediately isomerizes to a secondary carbenium ion. The latter is also expected to quickly stabilize by deprotonating to the framework, although this part of the reaction falls beyond the scope of the current study.

Figure 10

Two-dimensional umbrella sampling free energy profile at 773 K for the protonation of 2-octene and reaction D2-1, β-scission of the 3-octyl carbenium ion.

Our results allow an assessment of the importance of various alkene cracking modes from a theoretical perspective. Cracking mode A is by far the fastest, as it concerns a transition between tertiary carbenium ions. Such a cracking mode is only possible for alkenes with at least eight carbon atoms and three branches. Once formed, these precursors were found to have a reasonable lifetime at 773 K within ZSM-5, but their formation rate might be limited due to steric constraints. Type B2 cracking modes are expected to be very important cracking modes at operating conditions, given the fact that they start from stable carbenium ions and have relatively low intrinsic activation barriers. Cracking modes B1 which start from a secondary carbenium ion and yield a tertiary carbenium ion are expected to be less important because we found that the geminal methyl substituted cracking precursors (with a quaternary carbon atom) are extremely unstable. Cracking mode C, a transition between secondary carbenium ions, is also predicted to be less important as the reactants will quickly rearrange into more stable tertiary carbenium ions. Cracking modes D2 and E2 both require the formation of a transition state with primary carbenium ion character, thus resulting in substantially higher intrinsic barriers. For cracking mode D2, the linear alkene reactant resides freely in the zeolite pores at the high temperature. In order to be reactive for β-scission, it first has to transform into a carbenium ion, which yields an additional contribution to the activation barrier for the D2 mode. The conclusions drawn here are in line with the experimental findings in the seminal papers of Weitkamp.[95,99]

Conclusions

Using advanced molecular dynamics simulations, catalytic alkene cracking pathways in H-ZSM-5 were studied at actual operating temperatures of 773 K. Alkene cracking is known to proceed via a complex reaction network with many possible reaction routes and often elusive intermediates. At first instance the nature and stability of C4–C8 alkene cracking intermediates were determined. Over the course of the simulations, regular transitions were observed between a physisorbed van der Waals complex, a physisorbed π-complex, and a chemisorbed carbenium ion. Alkoxides, on the other hand, were not sampled and appear to be unstable at cracking conditions due to high entropic penalty for forming a covalent bond with the framework. Linear alkenes preferentially occur as a physisorbed π-complex or van der Waals complex. The latter is stabilized due to a higher mobility in the zeolite pore system. The stability of linear, secondary carbenium ions increases with carbon number, as the positive charge can be better delocalized over the carbon atoms along the chain. Linear alkenes first need to overcome the barrier for protonation before cracking reactions can occur. Free energies of activation for interconversion between a physisorbed alkene and a carbenium ion were determined using the metadynamics approach and amount to ca. 50 kJ/mol. For branched alkenes, carbenium ion intermediates were found to be more stable than physisorbed alkenes at cracking temperatures. No clear influence of the chain length could be observed. In contrast to linear carbenium ions, tertiary carbenium ions are undoubtedly very stable, long-living intermediates, which can act as cracking reactants. Secondary carbenium ions with geminal methyl branches though were shown to undergo rapid isomerization or β-scission into more stable carbenium ions.

In a next step, intrinsic barriers for β-scission reactions of mode A, B2, D2, and E2 were determined using the umbrella sampling molecular dynamics technique. Overall, the importance of the various cracking modes was found to be determined by the interplay between the reactant carbocation stability and the transition state (or product) carbocation stability. Regardless of the formed products, similar cracking barriers are expected for reactions of the same β-scission mode (i.e., reactions with the same carbocation transition). Cracking modes A (3° → 3°) and B2 (3° → 2°) were found to be dominant at 773 K, given the higher stability of branched cracking precursors and the lower free energy barriers for β-scission. Cracking modes B1 (2° → 3°), C (2° → 2°), D2 (2° → 1°), and E2 (3° → 1°) are expected to be less important. Mode D2 and E2 involve a transition state in which a unstable primary carbenium ion would be formed. Modes B1, C, and D2 start from metastable, secondary carbenium ion reactants, which will typically undergo rapid isomerizations to more stable tertiary carbenium ions. While cracking mode A has a very low intrinsic barrier, the highly branched cracking precursors are very reactive and their formation may be sterically hindered in the zeolite pores.

The B2 cracking mode, starting from a very stable tertiary carbenium ion, was also studied with static DFT calculations–the standard procedure to determine reactivity differences in current literature. Static cracking barriers were highly dependent on the particular orientation of the reactants and transition states in the zeolite channels. Furthermore, static calculations pointed to the formation of alkoxide products, which were proven to be nonexistent at high temperature. Standard DFT calculations are not capable of correctly predicting the stability of various cracking intermediates and do not account for the conformational freedom along the alkene cracking reaction profile. In contrast, consistent free energy barriers were obtained from umbrella sampling simulations which account for the full conformational mobility of all species at cracking temperatures. Our results clearly demonstrate the importance of simulations at operating conditions to assess the nature and reactivity of alkene cracking intermediates.