Unravelling Channel Structure–Diffusivity Relationships in Zeolite ZSM‐5 at the Single‐Molecule Level

Abstract The development of improved zeolite materials for applications in separation and catalysis requires understanding of mass transport. Herein, diffusion of single molecules is tracked in the straight and sinusoidal channels of the industrially relevant ZSM‐5 zeolites using a combination of single‐molecule localization microscopy and uniformly oriented zeolite thin films. Distinct motion behaviors are observed in zeolite channels with the same geometry, suggesting heterogeneous guest–host interactions. Quantification of the diffusion heterogeneities in the sinusoidal and straight channels suggests that the geometry of zeolite channels dictates the mobility and motion behavior of the guest molecules, resulting in diffusion anisotropy. The study of hierarchical zeolites shows that the addition of secondary pore networks primarily enhances the diffusivity of sinusoidal zeolite channels, and thus alleviating the diffusion limitations of microporous zeolites.

by repeating the latter part, i.e., the cycle of ethanol extraction and recovery by diethyl ether and ethyl acetate. The purity was checked by 13 C NMR. TPAOH was used as received for the synthesis of boriented zeolite crystals. TPAOH was used as received for the synthesis of b-oriented crystals.
Synthesis of a-oriented and b-oriented silicalite-1 seed crystals. For the synthesis of the a-oriented silicalite-1 crystals with a leaf-like shape, a solution with a ratio of 6 TEOS: 0.75 timer-TPA 3+ : 3.75 KOH: 1425 H2O was prepared after stirring for 8 h. The solution was then poured into an autoclave through a filtration paper. The autoclave was placed in a rotational oven at 448 K for 24 h. The preparation of boriented silicalite-1 crystals followed Yoon's method. 2 The round, silicalite-1 crystals were synthesized from a solution composition of 6 TEOS:0.9 TPAOH:620 H2O. The synthesis solution was prepared by adding TEOS to the solution containing TPAOH and H2O. The mixture was transformed into a clear solution after stirring in a sealed liner for 24 h at room temperature. The clear solution was filtered with filter paper and charged into a Teflon-lined autoclave. The hydrothermal reaction was carried out in a rotation oven at 423 K for 12 h.
Manual assembly of an a-oriented and a b-oriented silicalite-1 monolayer. Typically, the quartz plates for the a-oriented and b-oriented zeolite ZSM-5 films were pre-coated with polyethylenimine (PEI) or pretreated with H2O2 (35 wt % aqueous solution), respectively. Then ca. 30 mg of a-oriented and boriented silicalite-1 crystals were put at the surface of quartz plates followed by pressing and rubbing using a finger with soft latex gloves. Subsequently, the loosely attached top layers were removed by gently wiping using glass wool for ca. 30 s. The a-oriented monolayers were calcined at 823 K overnight to remove the organic layers before secondary growth. The b-oriented monolayers were calcined for 2 h at 423 K to enhance the strength of the attachment.

Preparation of a-oriented and b-oriented
Zeolite ZSM-5 films on quartz plates. The preparation of the a-oriented and b-oriented zeolite ZSM-5 films (thickness= ca. 800 nm) followed the same protocol reported by our group using the secondary growth method (see details in Supplementary Experimental Procedures). 3 The silicalite-1 monolayers supported on quartz plates were intergrown by the secondary growth method (SGM) at 448 K for 24 h. Typically, a solution comprising a required amount of Al2(SO4)3·18 H2O, 0.77 g H2SO4 (10 wt% aqueous solution) and 6.9 g of H2O was added to a sodium silicate solution made by mixing 0.977 g of sodium silicate and 7.31 g H2O as well as 0.32 g ethanol.
After another 30 min of stirring at room temperature, a clear solution (Si/Al= 125) was obtained and was transferred without filtration to a 20 mL Teflon-lined autoclave with the a-oriented and b-oriented crystals seeded substrates vertically placed at the bottom of the same autoclave. After the reaction, the autoclave was removed from the oven and quickly cooled to room temperature by immersing in cool H2O. As-synthesized zeolite membranes on the substrates were carefully removed from the autoclave and washed with copious amounts of deionized H2O and dried in air at 333 K.
After the synthesis, the films ( Figure 1) were first treated with 0.2 M NH4F solution for 2 h without stirring to remove the outermost amorphous silica layers that can block the channels of the membranes. Then the films were calcined at 823 K overnight with a ramp rate of 1.5 K/min under an airflow to remove the organic additive, i.e., ethanol. The H-form sample was obtained by three times repeating the ion exchange of the calcined Na-type sample with a 1 M NH4NO3 solution (Acros Organic, 99+%) at 333 K overnight followed by calcination at 823 K for 6 h with a ramp rate of 1.5 K/min. The films with hierarchical pores were prepared using the fluoride etching method. 4 Typically, the as-synthesized zeolite films were put into a 40 wt% NH4F aqueous solution and react at 333 K for 20 min under mechanical stirring. The etched products were thoroughly washed with distilled water after the fluoride treatment.

Section S2. Physicochemical characterization
The images of the zeolite monolayers and films were examined using a scanning electron microscope (SEM) from XL-30 (Philips) operating at an accelerating voltage of 5 kV. Before measurement, the surface of the materials was coated with a Pt layer of ca. 5 nm to avoid charging effects. X-ray diffraction (XRD) was used to confirm the orientation of the as-synthesized zeolite monolayers and films. The XRD patterns were collected using a Bruker D2 Phaser (2 nd Gen) instrument using a cobalt radiation source, Co kα = 1.789 Å. UV/Vis diffuse-reflectance spectroscopy (DRS) was applied to determine the absorption property of the zeolite ZSM-5 films after the oligomerization of furfuryl alcohol (FFA). The UV/Vis DRS spectra were collected using an AvaSpec 2048L spectrometer connected to a high-temperature UV/Vis optical fibre probe, which was used to collect spectra in reflection mode. Focused ion beam-scanning electron microscope (FIB-SEM) images were recorded on an FEI Helios NanoLab G3 UC (FEI Company) instrument. A protective layer of Pt was deposited on top of the region of interest before milling trenches on either side. Slices were milled perpendicular to the surface before SEM images were recorded (5 kV, 50 pA). Atomic force microscopy (AFM) micrographs were recorded using a Bruker Multimode 8 with a n-type silicon tips with a hard diamond-like coating and aluminium background (having a resonance frequency of 160 kHz and a force constant of 5 N/m). The data was post-processed in Gwyddion.

A. Details of experimental procedures
The protocol for the fabrication and loading of the reaction cell used in the experiments is shown in Figure S2. Single-molecule localization experiments were performed using a custom-made fluorescence microscope setup ( Figure S3). An inverted epi-fluorescence wide-field fluorescence microscope (Nikon Eclipse-TI) with a 100× oil immersion objective (NA=1.49) was used. A diode-pumped solid-state laser (Cobolt Jive-100, 560 nm) provided 54 mW to the sample. Fluorescence microscopy images were recorded as movies using an electron-multiplying charge-coupled device (EMCCD) camera (Andor iXon 897), after passing through a bandpass filter (Chroma, ZET560/10x), dichroic mirror (Chroma, ZT560rdc), long-pass filter (Chroma, ET575lp), quarter wave plate (Newport, 10RP44-1) and a 2.5× camera lens. Time-lapse videos were recorded with a resulting field of view of 32 × 32 μm 2 with 64 × 64 nm 2 per pixel (512 × 512 pixels), and a frame time of 30 ms. The reaction was performed in a Gene Frame cell (25 µL, 1 × 1 cm 2 , Thermo Scientific, see details in Figure S2). The cell consists of a Gene Frame, which forms a sealed liquid microscopy cell after adhesion to a coverslip and microscope slide on each side. The oligomerization of FFA was performed on activated zeolite ZSM-5 films loaded on the top of the cover glass in the gene frame cell. The zeolite materials were exposed to furfuryl alcohol (FFA, 99%, Sigma Aldrich), previously diluted in 1,4-dioxane, to achieve the desired catalytic activity. The optimal concentration (see details in Table S1) of FFA for high-resolution imaging was determined in a series of concentration-dependent measurements. The sample was equilibrated for at least two hours and exposed to maximum laser power for 5 minutes (10.9 mW through objective) for fluorescence photobleaching before the measurement.

B. Data Analysis
The recorded SMLM microscopy movies were analysed with the DoM plugin (Detection of Molecules, https://github.com/ekatrukha/DoM_Utrecht) for ImageJ. 23 The localization of fluorescent events was done by independent classification of each frame into emissive spots and background. A list of initial emitter positions was determined with a sub-diffraction limited spatial resolution by fitting a 2-D Gaussian using the Levenberg−Marquardt least-squares algorithm. For trajectory analysis, molecules were allowed to blink (i.e., the molecule is fluorescing intermittently) for maximally 3 consecutive frames (90 ms) and travel 6 pixels (384 nm) between consecutive localizations (tracking result in Supplementary Movies 1-4). Only trajectories with more than four localizations were considered to ensure sufficient displacements per trajectory for MSD analysis and to remove unphysical trajectories originated from incorrect localizations. Trajectory classification, analysis and plotting were done in MATLAB (The Mathworks) using DiffusionLab, a software developed in our group for the classification and motion analysis of single-molecule trajectories. Previous work of our group has demonstrated the possibility of using a machine learning approach to group a large set of trajectories into populations with the same motion behaviour. The classification of the trajectories was done following a hierarchical decision tree built from a training set. 24 Inspired by this work, we manually built a decision tree to group the trajectories in this work in a rational manner. Motion analysis was performed with mean squared displacement (MSD) analysis of individual trajectories or a set of trajectories, i.e., a population. The diffusion coefficient and localization error were obtained from a linear fit of the MSD curve. MSD analysis on individual trajectories was done including the first 25% of the delay times and at least three points. Only the first four points were used in the fit of the population MSD, because the number of trajectories contributing to the MSD was constant in this fit range.

C. Diffusion coefficient estimation
General notes for diffusion coefficient estimation. The free diffusion model in two dimensions were used for all motion types. We find that the population-averaged mean squared displacements are linear at short delay times (Figures 3c-d and 4b-c), indicating Brownian-type diffusion. One could expect that the channel geometry would give rise to one-dimensional diffusion. Nevertheless, a reason that we observed two-dimensional diffusion (Figure 2b) could be that the fluorophores travel through two or more different crystal domains with a different orientation in the observation plane. These domains are in the order of ~ 5 µm, therefore we don't expect this to occur very frequently and severely affecting the diffusion constant. We also fit the hybrid trajectories with a free diffusion model. A hybrid trajectory consists of mobile and immobile segments, due to transient trapping of the fluorophore. In twodimensional diffusion, the population-averaged MSD curve is linear indicating normal diffusion. Here, the slope of the MSD curve is a weighted average of the mobile and immobile segments. In threedimensional diffusion, the number of localizations per trajectory is heavily dependent on the length and number of immobile segments, because the fluorophore can diffuse out-of-focus during a mobile segment. Thus, only long trajectories with many immobile displacements contribute at long delay times.
This effect can give a plateau in the population-averaged MSD curve at long delay times, even though the fluorophores are not confined, which would give a similar effect. To prevent a bias in the measured slope of the MSD curve, the part of the MSD curve that considers displacements over longer delay times than the duration of the shortest trajectories must be excluded.
We determined the diffusion coefficient for individual trajectories and populations of trajectories with MSD analysis, as summarized in Table S2. To ensure sufficient displacements per trajectory for MSD analysis and remove artificial trajectories linked from random false localizations, we only considered trajectories with at least five localizations. Individual trajectories were analysed by fitting the mean squared displacement "MSD" as a function of lag time , 25 with the diffusion coefficient, and the localization error. This relation considers two-dimensional normal diffusion with a fixed localization error. The lag time is related to the number of frames over which the displacement spans and the frame time ∆ as = ∆ . The least squares fit was weighed by the number of squared displacements that contributed to the mean, and only the first 25% of the trajectory (or at least three data points) were taken for the fit. Due to statistical dependence of the MSD values, data points at longer t are strongly correlated and reduce the accuracy of the MSD fit. 26 The optimum number of data points to include in the fit is only known for trajectories without blinking, 27 therefore we used a fixed fraction of 25%.

Diffusion coefficients from individual trajectories and all trajectories in a single population.
Diffusion coefficients from individual trajectories were obtained by linear fitting of the MSD curves of each trajectory in a sample. The average motion of all trajectories in a single population via the "ensemble-averaged" (or population-averaged) MSD curve, which is constructed from the displacements of all trajectories. We fit equation (7) to the first four data points of the ensemble-averaged MSD curves. Only the first four data points were fit, because all trajectories had at least five localizations, from which we could extract displacements over up to at least four frames. At longer delay times, the MSD curve becomes irregular as not all trajectories in the population contribute to it anymore. These delay times were excluded from the fit.

Effective diffusion coefficient.
Comparison with the previously reported ensemble diffusion coefficients, e.g., obtained with PFG-NMR, requires the computation of an effective diffusion coefficient. 29 This is the weighted mean of the diffusion coefficient of each population following 29,30 eff = mobile mobile + hybrid hybrid with and the fraction of guest molecules and the diffusion in population X, respectively. As the diffusion coefficient of the immobile population is zero by definition, we leave this population out of the equation. We estimate weight as the average fraction of molecules per frame, i.e. the average localization density of a population. We compute the molecules per frame and include the frames between the first and last localization when the molecule is off. We used the standard error of the mean to compute the error in and the uncertainties were propagated following the well-known formulas. 31

Section S4. Time-dependent density functional theory simulations
The static and dynamic adsorption maxima (Table S3) for the neutral and protonated oligomers were calculated using the time-dependent density functional theory (TD-DFT) simulations. The ground state structures of a series of neutral and charged furfuryl alcohol oligomers were obtained using DFT method as implemented in version 6.1 of CP2K software. 5 The gas-phase geometries were optimized using GAPW method 6 with CAM-B3LYP functional 7,8 and 6-311++G** basis set. Plane-wave cut-off value of 600 Ry was used. The structures of selected oligomers derived from the literature and are shown in Figure S4. 9,10 For each of the oligomers, several possible conformers were generated and their stability both in a vacuum and within the zeolite (vide infra) was tested at CAM-B3LYP and BLYP-D3 level of theory, respectively. We find that while in a vacuum the "trans" type of oligomers is the most stable, the "cis" form is preferred in the zeolite (Table S4). Therefore, for all subsequent calculations we used "cis" type of oligomers.
Vertical electronic excitation energies were determined by using Time-Dependent DFT (TD DFT) as implemented in ADF package 11-14 at the CAM-B3LYP/TZV2P 15 level of theory, which has proven to be accurate for modelling of absorption bands of five-membered heterocyclic compounds 16 . The ability of CAM-B3LYP functional to correctly reproduce exited states of furfuryl alcohol and its derivates was further confirmed by comparing experimentally measured absorption maximum of furfuryl alcohol in water. The experimental value 217 nm 17 is very close to TDDFT computations (220 nm) providing that a polarisable continuum model (PCM) for water is used, therefore we conclude that CAM-B3LYP functional can be used to resolve main spectral features of UV/Vis absorption of the furfuryl alcohol oligomerization reaction products.
To account for a flexibility of the UV/Vis-active compounds and confinement arising from the zeolite network, we followed the same procedure as Hemelsoet et al. 18  The unit cell dimensions used for each simulation are summarized in Table S5. The system was initially allowed to equilibrate for 1 ps, which was followed by a 5 ps production run. The time-averaged spectra were obtained by taking 50 different snapshots (every 200 steps), on which separate TD DFT calculations were performed. During TD DFT simulations, the zeolite structure was not taken into account due to enormous computational cost this would require. Finally, the averaged spectrum from all snapshots was computed.

Supplementary videos 1-2:
wide-field movie of fluorescent molecules diffusing through the sinusoidal channels over a b-oriented zeolite film. Trajectories (lines) and localizations (circles) are overlaid in post-processing. Green circles mark localizations with a width (standard deviation of the Gaussian fit function) deviating less than 30% from the expected 128 nm width, while localizations with a larger deviation are marked red. Temporal resolution: 30 ms per frame. Playback at 15% real-time. In supplemental video 1, a zoom-in of a region with a high density of mobile and hybrid trajectories is shown. Supplemental video 2 displays the same time lapse without zoom.

Supplementary video 3:
wide-field movie of fluorescent molecules diffusing through the straight channels over an a-oriented zeolite film. Trajectories (lines) and localizations (circles) are overlaid in post-processing. Green circles mark localizations with a width (standard deviation of the Gaussian fit function) deviating less than 30% from the expected 128 nm width, while localizations with a larger deviation are marked red. Temporal resolution: 30 ms per frame. Playback at 15% real-time.

Supplementary video 4:
wide-field movie of fluorescent molecules diffusing through the sinusoidal channels over a b-oriented, hierarchical zeolite film. Trajectories (lines) and localizations (circles) are overlaid in post-processing.
Green circles mark localizations with a width (standard deviation of the Gaussian fit function) deviating less than 30% from the expected 128 nm width, while localizations with a larger deviation are marked red. Temporal resolution: 30 ms per frame. Playback at 15% real-time

Supplementary video 5:
wide-field movie of fluorescent molecules diffusing through the straight channels over an a-oriented, hierarchical zeolite film. Trajectories (lines) and localizations (circles) are overlaid in post-processing.
Green circles mark localizations with a width (standard deviation of the Gaussian fit function) deviating less than 30% from the expected 128 nm width, while localizations with a larger deviation are marked red. Temporal resolution: 30 ms per frame. Playback at 15% real-time. The zeolite films were grown from the same gel solution, so the difference in the a-oriented and boriented films would be solely in their orientations. The a-oriented zeolite films are with the sinusoidal and straight channels being perpendicular (z-axis) and parallel (xy-plane) to the substrates, respectively.
Similarly, the b-oriented zeolite films are with the straight and sinusoidal channels being perpendicular (z-axis) and parallel (xy-plane) to the substrates, respectively. Therefore, the channels in the z-axis are with the identical sizes and geometries, while those in the xy-plane are with the identical sizes and geometries. Due to the polarization of the laser and uniform orientation of the films, solely the molecules align in the channels in the xy-plane would be excited by the 560 nm laser. This allows us to exclusively investigate the diffusion behaviour of molecules in the channels with the same size and geometry in the xy-plane.  (Table S1). c) The reaction cell is closed with a glass slide adhered to the gene frame.  Table S2. Additionally, the stochastic single turnover dynamics were quantified.

Time-dependent quantification of the stochastic single turnover dynamics
To determine the catalytic behaviour of our zeolite thin films, we counted the number of new trajectories starting per unit of time, assuming that each start represents the formation of a new product molecule. 32 We found that the turnover rate increased in the initial stages of the experiment for both channel orientations and approached a constant value after approximately 2 hours of reaction (Figures 6a-b). A previous study of large ZSM-5 zeolite crystals showed a continuous increase of stochastic single turnovers as a function of time, which was attributed to the long diffusion path of ~20 µm in the large crystals. 32 However, our zeolite thin films have a short diffusion path (~800 nm), so FFA molecules can readily diffuse into the micropores of the zeolites and reach a steady-state concentration distribution within approximately 2 hours as shown in Figure S6c. Additionally, a sharp decrease in turnover rate was observed in the first 2 min of each measurement (not shown), indicating photobleaching of the molecules. To avoid this effect, 5 min of laser illumination was performed before each diffusion measurement.
The different reactivity in the straight and sinusoidal channels can affect the mobility of the molecular probes. As a result of a different reactivity, the fluorescent oligomers can have a different length, which directly affects the probe's mobility. We argue that the effect of the reactivity on the obtained diffusion constants is relatively small if we make two assumptions: a higher reactivity would result in a longer oligomer and this would lead to a lower diffusion constant. If the observed difference in probe diffusivity between the straight and sinusoidal pores were due to longer probes as a result of a higher reactivity, we expect a negative correlation between the reactivity and diffusivity. That says, a lower diffusion coefficient would be obtained for the straight channels. However, we find a positive correlation between the reactivity ( Figure S6) and diffusivity (Table S2) in the parent samples. Therefore, we conclude that the different reactivity is not a dominant factor in the probe diffusion.

Time-dependent density functional theory calculated absorption spectra of furfuryl alcohol oligomers
Time-dependent density functional theory (TD-DFT) calculations were performed to identify the UV-Vis absorption bands of furfuryl alcohol oligomers, and the results are summarized in Table S3 and Figures   S7-8. First, the absorption spectra of both protonated and neutral forms in the vacuum were considered.
The absorption maximum of the protonated species was always at longer wavelengths than the neutral form. As shown in Figure S7, the absorption bands of the neutral oligomers increased progressively in wavelength from 209 nm (monomers) to 475 nm (tetramers) and the bands of the pronated oligomers from 219 nm (monomers) to 574 nm (tetramers). From comparison with the experimental data ( Figure   S5a), it is evident that the absorption band around 460 nm corresponds to conjugated species, namely trimeric and tetrameric species. This is a guideline, because absorption bands of confined molecules can shift with respect to their unconfined equivalents as has been demonstrated by Hemelsoet et al. 18 The influence of the zeolitic framework on the oligomers trapped in the pores of the MFI zeolite was examined using ab-initio MD simulations. The results are compiled in Figure S8. protonated form. The only species that absorbs above 500 nm is the highly conjugated tetramer with a calculated absorption maximum at 544 and 614 nm for its neutral and protonated form. As we shown in Figure S8, the absorption maximum of the oligomers increases ca. 100 nm with the addition of an extra monomer. Therefore, we expected that the absorption for slightly larger oligomers, e.g., pentamer, would be with an absorption maxima > 644 and 714 for the neutral and protonated forms, respectively. The experimental results, however, show that the highest absorption observed for this reaction over zeolites is ca. 620 nm. This unambiguously suggest that the largest and dominant stable species formed during the furfuryl alcohol oligomerization that can be excited by a 560 nm laser, is the furfuryl alcohol tetramer.
However, as shown in Table S3, the adsorption peaks are further broadened by presence of different conformations. As the result, calculated standard deviation for the dynamic spectra is up to 32 nm.
Additionally, in our model we assume that oligomers are located in the straight zeolite channels. In a real system, they are formed in both straight and sinusoidal channels, which will result into further broadening of adsorption maxima. Thus, it cannot be excluded that using 560 nm laser we have excited also other species such as trimer. The tetramers, however, would be the most efficiently excited species as they show adsorption maximum close to 560 nm.   Tables S3-S5.

Determination of pixel jump and blinking gap
We use a tracking algorithm to identify single-molecule localizations that are likely the same molecule based on their vicinity in space and time. These localizations are grouped into what we call a "trajectory".
In the grouping process, two input values are important: the "pixel jump" and the "blinking gap" ( Figure   S9). The blinking gap accounts for the possibility that the fluorescence of a single molecule "blinks", i.e.
sometimes turns off for a few tens of ms. 33 The blinking gap is the maximum number of dark frames by which we allow a trajectory to be interrupted. Similarly, the pixel jump specifies the maximum spatial separation between two consecutive localizations in a trajectory.
Finding the optimal set of cut-off values is challenging, and even the most optimal set can limit the range of diffusion coefficients that can be reliably extracted from the experiments. If the pixel jump value is too small, the algorithm does not recognize that two consecutive localizations originate from the same molecule when they are further spaced apart from the pixel jump ( Figure S9a). The trajectory is falsely cut, and our motion analysis will measure a reduced average diffusion constant. Instead, if the pixel jump is too large, we will connect localizations of different molecules that are in each other's vicinity, which will increase the measured average diffusion constant in our analysis. Therefore, we aim at a pixel jump value that minimizes the probability of either scenario.
We have manually determined the optimal pixel jump and blinking gap by careful inspection of the trajectories overlaid with the recorded movies. Where the tracking algorithm cannot distinguish between localizations from different molecules, we humans sometimes can make an educated guess about the localizations that belong to the same molecule, for instance by considering the brightness of the localizations. By making such educated guesses, we select the largest value of the pixel jump and blinking gap that does not connect localizations of different molecules. We use these values in our analysis, which are given in Table S7. Figure S9. Examples of the grouping process by the tracking algorithm. a) The first example demonstrates the grouping process for a molecule that does not blink, and we only have to consider the pixel jump. When the pixel jump is shorter than the displacement between two subsequent localizations, such as between localization two and three, it is recognized by the algorithm. However, if the displacement is larger than the pixel jump, which happens between localization three and four, the displacement is not recognized. The trajectory is then falsely cut and we obtain two measured trajectories. b) In the second example, the molecule blinks and we also have to consider the blinking gap. Depending on its value, the algorithm still recognizes a displacement even though the molecule is off for one or more frames. The blinking gap is the maximum number of dark frames allowed and its value is one in this example. Thus, a pair of localizations with displacement shorter than the pixel jump will be recognized, if the number of dark frames between the localizations does not exceed one. For instance, between localizations one and three, the molecule blinks for one frame and the displacement is recognized. However, between localization three and six, the molecule is dark for two frames, which is longer than the blinking gap of one. The trajectory is falsely cut even though the displacement is smaller than the pixel jump, and we obtain two measured trajectories.

Modelling of the effect of the pixel jump on the measured diffusion constant
Depending on the pixel jump value, our tracking algorithm may tend to underestimate the diffusion constant of the fastest molecules in the sample. These molecules would show large displacements, which the tracking algorithm fails to recognize. When this happens frequently, the probability decreases that all localizations of a single molecule are grouped in the same trajectory. This probability is an important parameter, which we coin the "trajectory recognition success probability". A probability close to one means that the algorithm is able to recognize the molecules, while a probability close to zero shows that it often misses displacements and is not able to group all localizations of a molecule in the same trajectory. We want the algorithm to recognize the full trajectory, so a probability close to one is desired. In this section, we will model the trajectory recognition success probability to estimate whether we can reliably measure the diffusion constant of the fastest molecules in the sample. We do not consider the false grouping of different molecules into a trajectory, because we have manually selected values of the pixel jump and the blinking gap at which this should not occur.
To model the trajectory recognition success probability, we need to compute the probability that the displacements of every single step in the trajectory are recognized by the tracking algorithm. We start by modelling the recognition success probability of a single step in the trajectory, which allows us to compute the recognition success probability of the full trajectory. First, we consider a molecule that is continuously fluorescent and does not blink. This means that all steps in the trajectory span one frame.
We assume 2-dimensional Brownian motion with a displacement probability in x and y given by a normal distribution with a standard deviation of √2 , with the diffusion coefficient and the delay time. 34 We assume that the error in the localizations is normally distributed with a standard deviation of . The recognition success probability of a single step in the trajectory in the case of no blinking 35 is given by where ∆ is the pixel jump and MSD is the mean squared displacement. The subscript for 1 indicates that the displacement of this step spans one frame, which corresponds to a molecule that does not blink, and the MSD depends on the frame time ∆ as MSD = 4 ∆ + 4 2 .
The probability that a full trajectory is recognized by the tracking algorithm is given by with the length of the trajectory in frames. This situation is depicted in Figure S10a.
Blinking affects the probability that a single-molecule trajectory is properly identified. During periods when the molecule is off, it cannot be localized, and localizations are missing. These so-called "localization gaps" in the trajectory result in steps that span more than one frame. We can generalize equations (1) and (2) for steps that span frames with MSD = 4 ∆ + 4 2 .
To compute the trajectory recognition success probability, we multiply all values of weighted in the exponent by their occurrence in the trajectory. The probability that a full trajectory is recognized by the tracking algorithm becomes We do not know the blinking statistics of the molecules. Therefore, we estimate the average ℎ, , values from the statistical distribution of localization gaps in the experimental data. The length of these localization gaps is per definition limited by the blinking gap that was used in the tracking. Therefore, we can only approximate the true blinking statistical distribution from the experimental data. Localization gaps that are longer than the blinking gap are not measured. We cannot estimate these gaps from the experimental data and do not consider these in our analysis. This is not a problem for the diffusion constant estimation, because these long localization gaps do not affect the measured diffusion constant.
When a step in a trajectory spans over more frames than the blinking gap allows, the trajectory is falsely cut. The cutting is done irrespective of the magnitude of the displacement and therefore does not lead to an underestimation of the diffusion constant. Altogether, we can safely omit localization gaps larger than the blinking gap in our statistical analysis. We find that the probability that a full trajectory of a molecule is recognized by the tracking algorithm strongly depends on its diffusion constant and the length of the trajectory. We can plot this as a probability map if we assume a constant localization error, pixel jump, and frame time ( Figure S11a).
We find a very low recognition probability for fast molecules (D = 10 -11 -10 -12 m 2 s -1 ), near-perfect detection for slow molecules (D = 10 -14 -10 -13 m 2 s -1 ), and a transition region between D = 10 -12 -10 -13 m 2 s -1 . The transition region is important because it shows the fastest molecules that can be recognized by the tracking algorithm. We see that trajectories with a shorter length have a higher probability to be successfully recognized by the algorithm. The transition region moves to a lower diffusion constant when the pixel jump is decreased from six to two pixels ( Figure S11a-b). We can understand this intuitive result if we again realize that a shorter pixel jump limits the maximum displacement we can measure.
This leads to a lower probability to successfully recognize the step, therefore a decreased probability to recognize the full trajectory. In the end, this means that the diffusion constant of the fastest molecules we can accurately estimate is lower at a shorter pixel jump.
Blinking affects the trajectory recognition success probability as well, and we find a broadening of the transition region due to blinking (Figures 11a and 11c). This broadening is a result of many different combinations of localization gaps that can occur in a trajectory within the same statistical distribution.
We also find that the transition region moves to a slightly lower diffusion constant. The decreases for a larger localization gap: 1 > 2 > 3 > ⋯, because the molecule can have a larger displacement, as it has more time to travel between the detections ( Figure S11c). Therefore, it has a lower probability to have a displacement that is shorter than the pixel jump, which moves the transition region to a lower diffusion constant.

Validating the pixel jump
To quantitatively compare the experimentally obtained diffusion constants, we must be sure that we do not underestimate the diffusion constant of the fastest molecules. If we do not underestimate the diffusion constant of the molecules, we call the pixel jump "good"-otherwise it is "too short". We can investigate whether the pixel jump is good by comparison of the experimental trajectories with the trajectory recognition success probability map. Every trajectory is a point on the probability map, and from the distribution of trajectories, we can learn whether the pixel jump is good or too short. The problem is that we do not know the true diffusion constant and trajectory length of a measured trajectory; however, we can learn how the trajectory behaves near the transition region by looking at simulated trajectories with a known and .
First, we need to learn how the distribution of trajectories on the probability map looks like if the pixel jump is good. To this end, we simulate trajectories with a diffusion constant and statistical distribution of localization gaps representative for the experiments, and we call these the "true" trajectories. We assume two-dimensional Brownian motion as described above. The simulated trajectories were first converted into localizations. These localizations were linked back to new trajectories, which we call the "measured" trajectories. The linking was done with the same pixel jump and blinking gap as was used in the experiment. We simulate true trajectories with a length of 40 frames and diffusion constant of 10 -13 m 2 s -1 and overlay this with the probability map of ( Figure S12a). According to our analytical model, the simulation parameters result in trajectories with ≈ 0.9. We see that most measured trajectories span 40 frames, like the true trajectories, which confirms that the pixel jump is good. A few measured trajectories were shorter than the true trajectories, because of the ~0.1 probability that the true trajectory is not successfully recognized at ≈ 0.9.
Now that we know how the distribution of trajectories looks if the pixel jump is good, we simulate trajectories for which the pixel jump is too short. We simulate trajectories with a diffusion constant slightly higher (5•10 -13 m 2 s -1 ) and far higher than the transition region (10 -12 m 2 s -1 ), and overlay the resulting "measured" trajectories with their probability map ( Figure S12b-c). The pixel jump is too short according to our model, since the simulated parameters would result in < 0.1 for both simulations. Because the longest displacements are not included in the measured trajectories, we expect that these trajectories become shorter and have a lower diffusion constant than the true trajectories. We indeed find this trend in our simulation results. Overall, this leads to many trajectories in the blue-colored region of the probability maps (with < 0.2), which is in contrast to the simulation with a good pixel jump where we didn't find any trajectories in this region. We will use this characteristic to estimate whether the experimentally measured trajectories are linked by a pixel jump that is too short. We compare the experimentally obtained trajectories with their probability map of ( Figures S13a-d) and do not find strong indications that the pixel jump is too short. Indeed, all trajectories are in or below the transition region with a > 0.2. If, on the other hand, we take a smaller pixel jump value of 2 pixels, we find clear signatures that this affects our estimates of diffusion constants (Figures S13e-g). Here, a significant portion of the trajectories is in the < 0.2 domain, which indicates that the pixel jump is indeed too short. Altogether, this analysis demonstrates a proper selection of the pixel jump that does not introduce errors in the estimated diffusion constant. analysis fitting the first 25% or at least three points. The value of the y-axis for (e-h) is an order of magnitude lower than that for (a-d).

Validation of the diffusion heterogeneity observed by single-molecule tracking
As discussed in the main text, linear fits of the MSD curves of each trajectory in a sample produce a continuous range of diffusion coefficients produce a continuous range of diffusion coefficients spanning eight orders of magnitude (Figure 2a), which reflects the strong heterogeneity within theoretically identical straight zeolite channels. It should be noted that the imprecision in localization also contributes to the span of diffusion coefficients, particularly at small diffusion constants. Nevertheless, Figure S14 shows a distribution of "displacements", which is largely independent of the delay time. Closer inspection reveals a small shoulder at long displacements that increases with the delay time (inset of Figure S14), which is attributed to molecular diffusion. Therefore, these results collectively suggest the existence of We find a spread of more than two orders of magnitude in the diffusion coefficient, which is due to heterogeneity in the underlying motion and imprecision of the diffusion coefficient estimation. The distribution of diffusion coefficients is skewed and contains outliers-also to negative values. In these conditions, the median of the distribution is preferred to compare the populations. Because the distribution is not normal, the standard error of the median could not be computed. Although correlation in the displacements could lead to a large spread in the estimated diffusion coefficients 28 , our simulation shows that the variation is about 1 order of magnitude. Therefore, the up to 6 orders of magnitude variation of diffusion coefficients for trajectories in the same type of trajectories reflects the true heterogeneity of diffusivity in zeolite channels.

Classification of trajectories
The procedures for classification of trajectories are shown in Figure S15, with a detailed list of corresponding trajectory properties listed in Table S6. We build on the trajectory properties and hierarchical decision tree described in Hendriks and Meirer et al. 24 to classify the trajectories into mobile, hybrid and immobile populations. Inspired by this previous work in which trajectories were classified using a machine learning approach, 24 Table S6.    The error bars indicate the standard error.

Morphology analysis of zeolite films after the introduction of secondary porosity
Several approaches, including dealumination (i.e., steaming and acid leaching) and desilication (i.e., alkaline leaching) have been widely applied to produce meso-and macropores. 36,37 However, these approaches also change the chemical composition and the related spatial distribution of silicon and aluminium within zeolites and thereof their acidities. To exclusively study the effect of the introduction of meso-and macroporosity on the diffusion properties, we have applied a recently developed method, i.e., NH4F etching, to introduce secondary porosity by equally removing silicon and aluminium. 4 focused ion beam-scanning electron microscopy (FIB-SEM) analysis proved the successful introduction of secondary porosity to the zeolite films. Figure S20 shows that the intercrystallite domains grown from secondary growth are severally etched, demonstrating that these domains are more vulnerable to a fluoride attack. Similarly, Qin et al. reported that small intergrown domains on the surface of zeolite ZSM-5 crystals were removed. 38 Moreover, the formation of macropores was also observed at the surface of zeolite crystals, as shown in the high-resolution scanning electron microscopy (SEM) images in Figure   S20. A closer look has found small discrete domains with sizes of ~50 nm at around the crystal edges, most likely demonstrating the formation of mesopores in the zeolite crystals. To clearly illustrate the formation of mesopores, we also have attempted to study the morphology of the cross-sections using FIB-SEM. 39 As shown in Figures S20c-d 4 As the same etching approach was utilized, we postulated that the secondary pore network in the zeolite films also possess a similar arrangement being ordered and well-aligned.
To confirm the formation of mesopores and further highlight morphological changes to the film morphology, atomic force microscopy (AFM) measurements were performed on the b-oriented zeolite films prior and post etching, as shown in Figure S21. We selected b-oriented zeolite films as a representative as the material is known for its flat surface, which is suitable for AFM measurements.
The parent sample ( Figure S21a) shows typical zeolite crystal step edges of 1.3 ± 0.4 nm high, as calculated from the line profiles shown in Figure S21b. This value is equal to the height of the building unit of MFI, i.e., a pentasil chain. 41 Furthermore, a surface roughness, represented by a calculated RMS, of 2.4 and 3.0 nm can be found for the full zoom-in micrographs. If we only consider one crystal plane, a roughness of 0.3 -0.5 nm was found. This altogether suggests the formation of a perfect b-oriented zeolite film that is free of other pore networks.
The micrographs of the etched sample ( Figure S21c) confirm the successful introduction of secondary pore networks. Interestingly, the surface around the formed pores seems to be rougher and displays a granular structure as seen in the zoom-ins. The overall surface roughness of the samples increases to significant higher values (7.0 and 12.0 nm, excluding the cracks) and the crystal step edges are hard to be found, although still visible at certain locations, like in the micrograph with smaller Z-contrast (dashed blue box). This micrograph also highlights the large number of secondary pores found at the crystals.
From the displayed micrographs, all visible secondary pores were measured over the largest cross section and represented in the approximated pore size distribution in Figure S21d, indicating a median pore size of 112 nm. It should be noted that this is the pore size found at the surface, as the AFM is only to measure top-down and the cone-like tip shape prevents entering the pores if they are smaller than the tip apex.   Section S7. Supplementary Tables S1-S9 a. M-parent denotes "parent" zeolite films with preferential "M" orientation before the introduction of secondary pore networks. M-hierar denotes "hierarchical" zeolite films with preferential "M" orientation after the introduction of secondary pore networks.
b. The movement of individual fluorescent molecules within the sinusoidal and straight zeolite channels were tracked in the b-oriented and a-oriented zeolite films, respectively.