Catalytic conversion in nanoporous materials: Concentration oscillations and spatial correlations due to inhibited transport and intermolecular interactions

Traditional mean-field (MF) rate and reaction-diffusion equations of chemical kinetics apply in weakly interacting systems with locally well-stirred and randomized distributions of reactants and products. These formulations do not account for spatial correlations or the impact of such correlations on particle-number fluctuations.^1,2 In fact, intermolecular interactions generally induce non-trivial short-range ordering in fluids. However, extensive analysis of equilibrium systems has provided substantial insight into the associated thermodynamic pair correlations which aids assessment of their effect on reaction kinetics. Coincidentally, this type of short-range ordering is also reflected in the presence of concentration oscillations for fluids near walls and in confined environments.³ There is less appreciation of the feature that for systems with limited mobility, such as occurs in crowded reaction environments, distinct and sometimes strong non-equilibrium correlations can be induced by the presence of reaction.^1,2 Examples include catalytic surface reactions under high-pressure conditions,^4,5 and catalytic conversion in nanoporous materials with inhibited transport due to narrow pores.^2,6 The most extreme case for the latter where such non-equilibrium correlations should be strongest is where narrow pores impose a single-file diffusion (SFD) constraint, i.e., no passing of reactant and product molecules.^2,6 The lack of a general theoretical framework to precisely determine non-equilibrium correlations and their effect on reaction kinetics poses a major challenge for reliable beyond-MF assessment of the reaction yield.

In this contribution, we develop stochastic models for first-order catalytic conversion reactions A → B, where reactants A in a well-stirred external fluid diffuse into a decoupled array of narrow linear pores in a nanoporous material.^2,7–13 Conversion of A to product B can occur anywhere inside the pores, but significant yield relies on efficient removal of product from the pores to facilitate further entry of the reactant. We focus on the SFD regime where product removal is most inhibited. However, we emphasize that our model and theoretical formulation, as well as our basic observations, extend to the regime where passing within the pore is possible.² 

Investigation of such reaction systems from the early 1990’s was initially motivated by extensive studies of catalysis in zeolites, noting that a large subset of these materials do indeed consist of very narrow (∼1 nm) decoupled linear pores.⁶ Furthermore, experimental analysis for selected zeolite systems revealed clear indication of the presence of SFD and its influence on the reaction kinetics.^8,14,15 Subsequent studies have exploited new experimental techniques to characterize SFD in these systems.¹⁶ Additional interest in reaction systems subject to SFD was motivated by more recent studies of liquid-phase reactions utilizing catalytically functionalized mesoporous silica nanoparticles (MSN).^2,17,18 In general MSN particles can have coupling between pores. However, our particular interest relates to studies for MSN synthesized with hexagonal arrays of decoupled parallel linear pores with length ∼100 nm which traverse the entire nanoparticle, and where pores are not connected and are narrow.¹⁷ We should further emphasize that this synthesis procedure readily produces pores with nominal diameters of ∼2 nm, and adsorption of species on the pore walls under reaction conditions can lead to even narrower effective pore diameters. Significantly, such MSN systems have been shown to induce SFD.¹⁸ 

Coarse-grained spatially discrete stochastic modeling (described in more detail below) of catalytic conversion subject to SFD has typically been applied in order to efficiently treat the entire reaction-diffusion process.^2,7–13 Behavior of such models was precisely characterized by kinetic Monte Carlo (KMC) simulation. For previous simpler models, which did not incorporate non-trivial intermolecular interactions, it was recognized that a steady-state reaction occurs with reactant (product) concentration strongly decreasing (increasing) into the pore.^7–13 Furthermore, the total concentration was constant for these simple models, and thus the reaction-diffusion system was characterized by exactly counter-balancing gradients. This corresponds to a classic counter-diffusion mode.^2,13,19

With regard to analytic approaches to complement and provide further insight beyond KMC analysis, the limitations of MF-type treatments for even these simple systems have been recognized. In particular, such treatments were found to greatly overestimate diffusion fluxes in the presence of inhibited diffusion.^2,12,13 Such treatments do account for the non-trivial collective many-particle aspects of inhibited transport and its interplay with reaction. On the other hand, the so-called hydrodynamic treatments also fail in that they underestimate diffusion fluxes near pore openings.¹² Such treatments apply strictly only in the regime of small concentration gradients. As an aside, here the term “hydrodynamic” is used in a broad context of interacting particle systems,^20,21 and in particular diffusive systems, rather than just for convective fluid flow. Recent work on analytic treatments has shown the effectiveness of generalized hydrodynamic (GH) treatments of transport in the presence of strong mesoscale concentration gradients.¹³ This GH terminology is borrowed from early studies of convective fluid dynamics going beyond hydrodynamic treatments to describe transport on shorter time- and length-scales.²² Again our use is in a broader sense considering diffusive interacting particle systems.

Our goal here is to extend previous simple reaction models to include steric intermolecular interactions which are present in real systems. These interactions induce concentration oscillations in an external fluid approaching the catalytic nanoparticle, and would induce radial concentration oscillations in wider pores (although this feature is not included in or relevant for our modeling of narrow pores). Perhaps unexpectedly, we show that these interactions do induce oscillations along the pore axes near the pore openings, a key feature which must be incorporated in our modeling. This feature implies that it is necessary to extend the standard concept of counter-diffusion modes (applying just for constant total concentration), and also to adapt the previous GH formulation to this more complex scenario. In Sec. II, we present our model for catalytic conversion in nanopores with SFD which incorporates intermolecular interactions. Of particular significance is development of a strategy enabling explicit simulation of processes just inside the pore, avoiding computationally expensive simulation of the surrounding well-stirred equilibrated fluid. Sec. II also presents basic KMC simulation results for concentration profiles inside the pore. In Sec. III, we develop our analytic GH formulation, and demonstrate that it can accurately capture behavior seen in KMC simulation analysis. Thus, the GH formulations provide insight into the failure of simple MF-type formulations. Conclusions are provided in Sec. IV.

We first provide a detailed description of our model which is illustrated schematically in Fig. 1. In the spirit of classic lattice-gas descriptions of liquids,²³ each pore is divided into a linear array of cells labeled n = 1 − L whose centers correspond to discrete molecular positions. For prescription of adsorption and desorption at pore openings, it is actually convenient to extend this linear array within the pore to a three-dimensional simple-cubic array of cells in the surrounding fluid. The cell spacing, a ∼ 1 nm, is regarded as being slightly smaller than molecular dimensions, so that nearest-neighbor (NN) cells cannot be occupied. For convenience, we often set a = 1 below. This steric exclusion constraint suffices to induce all of the features (equilibrium spatial correlations, concentration oscillations, etc.) associated with more general molecular interactions.²³ The key ingredients of our model are as follows:

• Reactants A adsorb into the pore from the surrounding external fluid. This process is described by hopping at rate h from cells just outside the pore to an end cell n = 1 or L, provided the end cell and its NN cell within the pore are empty.

• Reactants A diffuse within the pore by hopping to NN empty cells at rate h (in either direction) provided that this creates no NN pairs of species.

• Reactants A convert to products B at rate k at any cell inside the pore.

• Products B diffuse within the pore by the same mechanism as for A. This prescription automatically imposes SFD, i.e., no passing of A and B in the pore.

• Reactants A and products B desorb from the pore to the surrounding external fluid, a process which is described by hopping from end sites of the pore at rate h to cells just outside the pore provided that such cells are empty, and also that all of the five NN cells of that target cell are also empty.

Thus, local diffusion within the pore in the direction along the pore axis is described by a single hop rate, h (and a corresponding low-concentration diffusion coefficient of D₀ = a²h). A central component of the analysis in Sec. III is to appropriately describe the corresponding chemical diffusion for finite concentrations in this multi-component system. Diffusion in the radial direction within pores is not relevant for the model. The exterior fluid is regarded as being in a well-stirred equilibrated state (corresponding to a lattice-gas with NN exclusion). We emphasize that this equilibrium assumption means that the associated diffusive or convective dynamics in the external fluid is not relevant for modeling. (As an aside, we note that one could regard this equilibrium state as being achieved by rapid effective hopping between neighboring cells subject to NN exclusion.) Another key feature of our model is that the exterior fluid has a large volume compared with the pores, so the desorbing product is quickly diluted and does not re-enter the pore. Thus, the external bulk reactant concentration, 〈A_b〉, matches the total external concentration, 〈X_b〉, and is a fixed constant. Finally, we emphasize that the equilibrium state of the external fluid is non-trivial with long-range ordering or crystallinity developing above a critical concentration 〈X_c〉 ≈ 0.209.²³ Consequently, we consider only the regime with short-range order for bulk concentrations 〈A_b〉 = 〈X_b〉 below 〈X_c〉.

Since in this model, reactants and products are “identical” in terms of interactions and diffusional dynamics, evolution of the total concentration corresponds to a pure diffusion problem for a single-component lattice-gas model with NN exclusion. The current study just focuses on steady-state behavior, so such evolution is not directly relevant. Nonetheless, we note that evolution is non-trivial even in the hydrodynamic regime of small concentration gradients given a non-trivial concentration-dependence of chemical diffusion in this model.²⁴ In the reactive steady-state of interest here, the total concentration matches that of an equilibrium model with NN exclusion. However, even this concentration distribution is non-trivial. The fluid + pore geometry induces concentration oscillations in the external fluid approaching the interface with the nanoporous material, and also a particularly complicated three-dimensional variation of the concentration near the pore opening. Furthermore, we shall see that there are also concentration oscillations within the pore along its axis within, but restricted to near the pore openings. All of these complex concentration variations will impact key adsorption and desorption rates at the pore openings, as discussed below.

Behavior of the above stochastic model can be assessed precisely by KMC simulation. The default treatment would simultaneously simulate behavior in both the pore interior and the external fluid. However, this approach is inefficient due to the large external fluid volume. Furthermore, it is unnecessary due to the assumed rapid equilibration of the external fluid. Thus, we are motivated to develop a strategy to enable explicit simulation of just the pore region while exactly accounting for the non-trivial coupling to the equilibrated external fluid. To this end, we first perform tailored simulations of the exterior fluid region to extract key adsorption and desorption parameters which will constitute the appropriate boundary conditions at pore openings for these stand-alone simulations of the pore region.

For adsorption, we first note that the concentration, 〈A₀〉, at cells just outside the pore, given that the end pore cell is empty, corresponds to the concentration of a fluid against the wall in a semi-infinite fluid system. Thus, we perform a simulation analysis of our lattice-gas model of the fluid with NN exclusion for a semi-infinite system, where the concentration only depends on the distance from the wall and exhibits strongly decaying oscillations away from the wall. Of most relevance, we find that the concentration, 〈A₀〉, is enhanced relative to the bulk concentration, 〈A_b〉. This enhancement is a natural consequence of the lower coordination of cells against the wall (with 5 neighbors which could possibly be occupied) relative to the coordination of cells in the bulk of the fluid (with 6 neighbors). See  Appendix A for further discussion and results for these concentration oscillations and enhancement at the wall, including a simple analytic estimate. This enhancement is quantified in Table I for a range of 〈A_b〉. Finally, we note that the adsorption rate at empty end cells of the pore (which also have empty NN cells within the pore) is given by R_ads = h〈A₀〉, and thus is not determined simply by the bulk concentration 〈A_b〉, but rather by 〈A₀〉.

For desorption, the presence of a particle at the end cell within the pore implies that the cell just outside the pore is empty. However, desorption requires that in addition all five cells adjacent to this cell are also empty. (The 2D analogue of these sites is denoted by * in Fig. 1.) Based on these observations, we perform additional tailored simulations of a semi-infinite fluid with one cell against a wall specified empty. These reveal a complicated three-dimensional variation of the concentration near the cell specified empty (in addition to the type of concentration oscillations approaching the wall away from this cell described above). See  Appendix B for further discussion. These tailored simulations allow determination of the conditional probability, Q₅, that these five additional cells are empty. Results for Q₅ are given in Table I. Then, it follows that the desorption rate from the filled end cell of the pore equals R_des = hQ₅.

As an aside, above we have described above the non-trivial and distinct concentration variations in the external fluid associated with both of our tailored simulations to extract adsorption and desorption parameters. Neither of these corresponds to the concentration variation in the external fluid under steady-state reaction conditions (which is just the equilibrium concentration of a lattice-gas with NN exclusion in a geometry corresponding to the fluid + pore system). However, we describe in  Appendix B how this latter concentration distribution can be reconstructed from the two distinct distributions extracted from our tailored simulations.

Below, we present KMC results of basic steady-state behavior. These and subsequent results are obtained from simulations just of the pore region with the appropriate non-trivial adsorption-desorption boundary conditions described in Sec. II B. However, we have confirmed in selected cases that results are consistent with large-scale simulations of the entire fluid + pore system. Fig. 2 shows typical steady-state concentration profiles in the pore for L = 30 with k/h = 0.001 and 〈A_b〉 = 0.2. Oscillations are apparent in both the total concentration and the reactant concentration near the pore openings. Thus, the steady-state does not correspond to a conventional counter-diffusion mode with constant total concentration and exactly counter-opposing gradients of A and B.¹⁹ However, we describe it as a quasi-counter diffusion mode since these conditions still apply away from the pore openings. See the supplementary material Fig. S1 for behavior with larger L where 〈A_n〉 ≈ 0 in the pore center.

With regard to total concentrations within the pore for L = 30 with k/h = 0.001 and 〈A_b〉 = 0.2, we specifically find that 〈X₁〉 = 0.321, 〈X₂〉 = 0.254, 〈X₃〉 = 0.279, 〈X₄〉 = 0.270, 〈X₅〉 = 0.273, etc., and a total concentration near the pore center of around 〈X_int〉 = 0.272. Clearly all of these values are substantially higher than in the bulk of the external fluid at 〈X_b〉 = 0.200, and also higher than the enhanced value of 〈X₀〉 ≈ 0.211 just outside the pore opening. This strong enhancement of concentration within the pore reflects the much lower coordination of cells within the pore (with 2 neighbors which could possibly be occupied) relative to the coordination of cells in the bulk of the fluid (with 6 neighbors which could be occupied). See  Appendix C for further discussion including a simple analytic estimate of this strong enhancement. The sudden transition from high-coordinated sites just outside the pore to lower coordinated sites within produces the concentration oscillations near pore openings as is evident in Fig. 1. We show in Sec. III that an accurate analytic description of this complicated behavior is possible within our GH formulation.

Deeper insight into reaction model behavior comes from an analytic formulation based on exact master equations for the stochastic process. Let 〈C_n〉 denote the probability that cell n in the pore is occupied by species C = A, B, or is empty E. It is also convenient to introduce the notation X = A or B for either type of species, so that 〈X_n〉 = 〈A_n〉 + 〈B_n〉 denotes the total concentration at cell n. Let 〈A_nE_n+1E_n+2〉 denote the probability that cell n is occupied by A and cells n + 1 and n + 2 are empty, etc. The NN exclusion constraint and conservation of probability impose various relations on these multisite probabilities.² The lowest-order evolution equations have the form

where

is the net diffusion flux of C = A or B from cell n to cell n + 1. Also ∇K_n = K_n − K_n−1 denotes a discrete derivative. Separate equations apply for the end cells, n = 1, 2 and L − 1, L, reflecting the non-trivial adsorption-desorption boundary conditions described above. See  Appendix D. The overall conversion rate of A to B is given by $Rtot=k∑n〈An〉$ simply reflecting the total amount of reactant inside the pore. Eq. (1) is not closed due to the appearance of triplet probabilities in J_C^n>n+1, but equations can be developed for such multisite probabilities generating a coupled hierarchy. See again  Appendix D.

Adding (1) for 〈A_n〉 and (1) for 〈B_n〉 leads to the pure diffusion equations

for the total concentration 〈X_n〉 = 〈A_n〉 + 〈B_n〉, for diffusion flux J_X^n>n+1 = h[〈X_nE_n+1E_n+2〉 − 〈E_n−1E_nX_n+1〉]. Again, separate equations are needed for end cells, n = 1, 2 and L − 1, L. In the steady-state, the spatial Markov property of 1D lattice models with NN interactions ensures the pair approximation factorization becomes exact, e.g.,

In obtaining the reduced expression after the last equality, we have also exploited NN exclusion. Using a similar relation for 〈E_n−1E_nX_n+1〉 together with the adsorption-desorption boundary conditions, one can solve exactly a coupled set of equations for 〈X_n〉 to recover the oscillations in the total concentration shown in Fig. 2. See  Appendix E. Such exact solution for steady-state 〈X_n〉 does not extend to the transient regime of pore filling, or to the individual reactant and product concentrations.

The fundamental challenge in solving the reaction-diffusion Equation (1) is to develop appropriate expressions for the diffusion fluxes, J_C^n>n+1. MF-type factorization approximations for probabilities of multi-cell configurations can fail dramatically. The site approximation neglects all spatial correlations and thus fails even to account for NN exclusion. Furthermore, it greatly overestimates diffusion fluxes for SFD, reactant penetration in the pore, and thus reactivity. The refined pair approximation accounts for NN correlations and thus excludes NN occupancy, but it still significantly overestimates diffusion fluxes and related quantities. See  Appendix D. Substantial additional insights into these shortcomings are provided below. An alternative hydrodynamic treatment applies for slowly varying concentration gradients, as mentioned previously. Thus, it is not geared to describe concentration oscillations occurring on the nanoscale near pore openings, but it is potentially relevant for description of longer mesoscale concentration variations deeper in the pore which do correspond to a classic counter-diffusion mode. The hydrodynamic diffusion fluxes satisfy J_C^n>n+1 = − D_tr∇〈C_n+1〉, where D_tr is the tracer diffusion coefficient for particles X.^2,12,19 However, for SFD, such D_tr are negligible, specifically decreasing to zero inversely with the pore length.^25–28 Consequently, this formulation greatly underestimates diffusion fluxes, reactant penetration, and thus reactivity for typical length pores.

Thus, another strategy is required to treat diffusive transport on the mesoscale, also accounting for concentration oscillations. A key ingredient which is motivated by generalized hydrodynamic (GH) treatments of fluids²² is to replace hydrodynamic transport coefficients with ones appropriate for a shorter mesoscale. In our case, these reflect distinct behavior near the pore openings where fluctuations in adsorption-desorption processes are prominent. Specifically, we replace D_tr with a spatially varying D_tr(n, n + 1) for each NN pair of cells which is enhanced near the pore openings (see below).¹³ In addition, to ensure the diffusion flux vanishes in the steady-state, we define fractional coverages 〈c_n〉 = 〈C_n〉/〈X_n〉 for C = A or B (and c = a or b) and adopt a specific GH form

Note that Eq. (5) automatically recovers the standard choice for conventional hydrodynamic counter-diffusion where J_C = − D_tr ∇〈C〉 in a continuum setting.^2,12,19

Next, we outline the determination of the GH D_tr(n, n + 1) from analysis of the so-called tracer counter-permeation (TCP).¹⁹ Here, a species labeled 1 enters a pore only from the left, and differently labeled species 2 (which is identical in terms of interactions and diffusional dynamics) enters only from the right. Otherwise adsorption and desorption are treated as for the above simulations incorporating non-trivial boundary conditions at the pore opening. The TCP simulations reach a steady-state with equal and opposite fluxes of magnitude J_TCP of 1 from left to right, and 2 in the opposite direction through the pore. See Fig. 3(a). Measuring the concentrations at different sites and equating the total flux with an expression of the type (5) allows extraction of the generalized tracer diffusion coefficients. Results are shown in Fig. 3(b) for L = 30 (and for larger L in the supplementary material Fig. S2). As expected, D_tr is naturally strongly reduced for higher total concentrations. Also, the GH D_tr(n, n + 1) decays to a non-zero plateau value, D_tr(plateau), in the pore center for sufficiently large L. Adapting previous studies which considered the overall tracer diffusivity for SFD in finite systems without NN exclusion^13,19,27,28 to account for NN exclusion in our model, we anticipate that

Here, 〈X_int〉 is the plateau value of the total concentration 〈X_n〉 in the pore center. For L = 30, we find that 〈X_int〉 = 0.272, 0.210, 0.134, and 0.059 for longer pores exceeds the external bulk fluid concentration 〈X_b〉 = 0.20, 0.15, 0.10, 0.05, respectively. This relation for D_tr(plateau) reasonably estimates precise KMC values reported in Fig. 3(c).

Finally, we remark that the above-mentioned overestimation of the diffusion fluxes by the site and pair approximations can be understood from the corresponding results for tracer diffusivity,

which far exceed D_tr(plateau) for typical L. Derivation of these results is indicated in  Appendix D.

Numerical solution of the GH reaction-diffusion Equation (1) is implemented incorporating the expression (2) for J_C^n>n+1 and our exact analytic solution for 〈X_n〉. The results almost exactly recover the individual reactant and product concentration profiles (including the concentration oscillations) obtained from KMC simulations shown in Fig. 2 for k/h = 0.001, L = 30, and 〈A_b〉 = 0.2. The degree of success of the GH theory for a range of k/h retaining 〈A_b〉 = 0.2 is shown in Fig. 4 focusing on the reactant profiles. Since 〈X_n〉 is recovered exactly, slight discrepancy in predicting reactant profiles is counterbalanced by a discrepancy of the same magnitude in prediction product profiles. To contrast the success of the GH theory, Fig. 4 also shows shortcomings of the pair approximation which predicts far too great a reactant penetration into the pore due to overestimation of the diffusion flux in the presence of SFD. See the supplementary material Fig. S3 for additional results. Since the total reaction rate, R_tot, for conversion of A to B simply reflects the total amount of reactant in the pore, success in predicting the reactant concentration profile automatically translates into success in predicting R_tot.

The above results indicate that our GH theory is well-suited to describe the regime of small k/h ≤ 0.001 where the reactant concentration exhibits slower mesoscale decay into the pore. For larger k/h where the reactant concentration decays more quickly on the nanoscale, the mesoscale GH treatment becomes somewhat less precise (although still reasonably accurate and certainly qualitatively correct). Actually for k/h ∼ 0.1, higher-order MF type approximations achieve comparable accuracy to the GH formulation. See the supplementary material Fig. S4.

To further illustrate the effectiveness of our GH formulation, one can also consider behavior for fixed k/h = 0.001, say, but varying the external fluid concentration 〈A_b〉. Naturally, this analysis necessarily incorporates the appropriate generalized tracer diffusion coefficients which depend strongly on 〈A_b〉, as shown in Fig. 3. Comparison of successful GH predictions with precise KMC results (and also generally inadequate pair approximation predictions) is shown in Fig. 5. Naturally, for lower concentrations, oscillations become significantly less prominent. Even the pair approximation becomes reliable for low enough concentrations where SFD constraints become less significant.

We have already provided one perspective on why MF-type approximations overestimate reactant penetration into the pore (and thus reactivity), specifically tying this feature to their overestimation of tracer diffusivity. Next, we provide an alternative perspective, and also a deeper understanding of the failure of the conventional MF type approximations. We emphasize that SFD in the presence of a reaction and also NN exclusion generates strong non-equilibrium spatial correlations. A direct consequence of these strong spatial correlations is the feature that the exact diffusion flux,

from (2) is far smaller than site or pair approximation predictions, and similarly for J_B^n>n+1. To restate this observation, these strong correlations imply that the triplet probabilities, 〈A_nE_n+1E_n+2〉 and 〈E_n−1E_nA_n+1〉, are far closer to each other than the site or pair approximation predictions.

In the site and pair approximations, neglecting oscillations in the total concentration, one has that

where G(x) = (1 − x)² in the site approximation and G(x) = (1 − 2x)/(1 − x) in the pair approximation. Thus, the large values of J_A^n>n+1 derive from the significant difference in the estimates of these triplet probabilities. This difference in turn derives from the significant difference between 〈A_n〉 and 〈A_n+1〉 near the pore openings given strong concentration variations in that region.

On the other hand, to understand exact behavior, it is useful to first note the exact identities

Here, we have used the feature that the site to the left (right) of A in the configuration A_nE_n+1E_n+2 (E_n−1E_nA_n+1) must be empty due to NN exclusion. Next, considering the central pair of cells n and n + 1 in the quartet configurations E_n−1A_nE_n+1E_n+2 and E_n−1E_nA_n+1E_n+2, we recognize that A is likely to hop back and forth between these two cells. This follows as the cells on each end of the quartet are specified empty ensuring that such motion is compatible with NN exclusion. This facile motion naturally tends to equalize these two probabilities. Precise results from KMC simulation analysis shown in Fig. 6 confirm this picture choosing a longer pore with L = 100 to clearly show behavior.

In summary, catalytic conversion subject to SFD produces strong non-equilibrium spatial correlations which occur in addition to correlations of thermodynamic origin due to intermolecular interactions (NN exclusion in our model). After accounting for equilibrium correlations and concentration oscillations in the total concentration, we show that a suitably refined GH treatment can capture both non-equilibrium and thermodynamic correlations. As a consequence, this formulation can reliably predict the mesoscale variation of the reactant concentration profile (as well as the concentration oscillations), and thus also predict the reaction yield. Our analytic formulation also provides deeper insight into the origin and nature of these correlations than is provided just from KMC simulation studies.

It should also be emphasized that our model is readily amenable to refinement and extension. One can relax the SFD constraint by allowing the exchange of A and B on second NN sites in the pore with rate p_exh, where p_ex reflects the passing propensity (and p_ex = 0 for SFD).² Passing reduces the strength on the non-equilibrium correlations that develop during reaction, so the GH formulation becomes even more accurate. Also, beyond treatment of just the initial stages of reaction, one can analyze the reaction yield for various specified fractions, f, of reactant converted to product in the external fluid, so now product can renter the pore. (Note that we assume a separation of time scales where the steady-state for a specific f is achieved on a short time scale compared to the overall reaction.) The overall reaction kinetics can be pieced together from a sequence of such simulations for increasing f.² 

See the supplementary material for a more comprehensive set of simulation results for both generalized tracer diffusivity and steady-state concentration profiles. Fig. S1: Steady-state concentration profiles for increasing pore lengths. Fig. S2: Generalized tracer diffusion coefficients, D_tr(n, n + 1), versus n for various pore lengths. Table SI gives values for D_tr(plateau). Fig. S3. Comparison of KMC, GH, and PA predictions for reactant profiles for various pore lengths. Fig. S4. Comparison of predictions from MF-type site, pair, and triplet approximations.

We acknowledge discussions with Igor Slowing and Marek Pruski motivating this study. We thank Tiago Oliveira for discussions on the theoretical formulation. This work was supported by the U.S. Department of Energy (USDOE), Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences, and Biosciences through the Ames Laboratory Chemical Physics program. The work was performed at Ames Laboratory which is operated for the USDOE by Iowa State University under Contract No. DE-AC02-07CH11358.

Determination of the adsorption rate for reactants into the pore in our model with NN exclusion requires analysis of the concentration variation approaching a planar wall in a semi-infinite lattice-gas model on a simple-cubic lattice with NN exclusion. Consistent with the notation of Sec. II, we let 〈X₀〉 denote the concentration in cells in the layer adjacent to the wall, 〈X₋₁〉 the concentration in cells in the next layer away from the wall, 〈X₋₂〉 the concentration in the next layer further away, etc., and 〈X_b〉 denotes the bulk concentration far from the wall.

Simple analytic estimation of this variation, and specifically of the (weakly) enhanced concentration adjacent to the wall, is possible utilizing a pair approximation. To this end, it is convenient to consider the semi-infinite equilibrated fluid as having arbitrary-range exchange dynamics described by a rate r, where exchange events are consistent with NN exclusion. In equilibrium, the corresponding flux of atoms from a cell adjacent to the wall to the bulk, J_w→b, and the reverse flux from the bulk to the wall, J_b→w, must balance. If P₇ is the probability of an empty cell in the bulk with all six NN cells also empty, then one has that

in the pair approximation. Let P₆ denote the probability that a cell against the wall, as well as all 5 of the NN cells, are empty. Then, one has that

in the pair approximation. Let us first assume that 〈X₋₁〉 ≈ 〈X_b〉, i.e., that concentration oscillations die out quickly away from the wall. Then, from the equality J_w→b = J_b→w, one obtains that 〈X₀〉 ≈ 0.2189 (versus the Monte Carlo simulation value 0.211) for 〈X_b〉 = 0.20. We also obtain 〈X₀〉 ≈ 0.1071 (versus the simulation value 0.106) for 〈X_b〉 = 0.10, etc. Thus, the pair approximation gives a quite reliable estimate of the (weak) concentration enhancement near the wall.

The above analysis can be refined to provide additional assessment of concentration oscillations away from the wall. The next level of analysis retains 〈X₋₁〉 as an independent variable, but assumes that 〈X₋₂〉 ≈ 〈X_b〉. Then, in addition to the equality J_w→b = J_b→w, one also balances fluxes between the layer of cells with concentration 〈X₋₁〉 and the bulk. Using the pair approximation, this yields two coupled equations for two unknowns, 〈X₀〉 and 〈X₋₁〉 in terms of 〈X_b〉. Their solution yields 〈X₀〉 ≈ 0.2192 and 〈X₋₁〉 = 0.1979, for 〈X_b〉 = 0.20. Thus, one predicts that 〈X₀〉 and 〈X₋₁〉 are 9% above and 1% below 〈X_b〉, respectively, versus simulation results which give values 6% above and 0.5% below 〈X_b〉, respectively. This pair approximation analysis also yields 〈X₀〉 ≈ 0.1072 and 〈X₋₁〉 = 0.0995 for 〈X_b〉 = 0.10, also mimicking the rapid decay seen in simulation studies. These results support the assumption in the simplest analysis that concentration oscillations decay quickly away from the wall.

The analytic treatment is readily further refined for an even more complete assessment of concentration oscillations. We have also performed a more complete Monte Carlo simulation analysis of the semi-infinite system with NN exclusion. However, 〈X₋₂〉, 〈X₋₃〉, etc., are very close to 〈X_b〉, so the above more limited analysis provides an essentially complete picture.

Our tailored simulations to extract adsorption and desorption parameters (described in Sec. II B) produce non-trivial and distinct concentration variations in the external fluid which might be regarded as a semi-infinite system. For the former, the concentration just depends on distance from the wall. For the latter there is a complicated three-dimensional variation with the strongest deviation from the bulk fluid concentration occurring around the cell specified empty just outside the pore. We argue that information from these tailored simulations provides boundary conditions at pore openings which allow simulation of the reaction model just in the pore region (which in turn recovers reaction behavior in the entire pore + external fluid system). From this perspective, one would also expect that information from the tailored simulations should allow recovery of the equilibrium concentration variations in the external fluid under steady-state conditions. We note that these equilibrium variations are distinct from those in tailored simulations for either adsorption or desorption parameters.

The tailored simulations for adsorption correspond to the situation where the end cell of the pore is empty, which occurs a fraction 〈E₁〉 = 1 − 〈X₁〉 of the time. Those corresponding to desorption correspond to the situation where this end cell in the pore is occupied, which occurs a fraction 〈X₁〉 of the time. Thus, we claim that the equilibrium distribution for the model is simply given by a corresponding weighted average of the distributions in the tailored simulations. This feature is illustrated schematically in Fig. 7, where we just show concentration variation along a 1D line of cells in the fluid which extend out from the pore opening. The ability to reconstruct the equilibrium distribution from the tailored simulations also reflects a spatial Markov field property of lattice-gas models with NN interactions²⁹ which applies not just for infinite systems, but also in more complex (e.g., pore + external fluid) geometries. We will elaborate on this feature in a separate paper dealing with more general models.

It is appropriate to provide further insight into the strong enhancement of total concentration in the center of the pore, 〈X_int〉, relative to that in the external bulk fluid, 〈X_b〉. The concentration in the center of long pores can be determined directly in terms of 〈X_b〉 by considering an infinite 3D lattice-gas model with NN exclusion suitably coupled to a 1D lattice-gas model with NN exclusion. Analogous to  Appendix A, this coupling is realized by direct exchange between the systems described by rate r, where exchange events are consistent with NN exclusion. In equilibrium, the corresponding flux of atoms from the 3D to the 1D system, J_3D→1D, and the reverse flux from the 1D to the 3D system, J_1D→3D, must balance. If P₇ denotes the probability of an empty cell in the 3D system with all neighbors empty as in (A1), then one has that

where a pair approximation estimate of P₇ is given in (A1). If P₃ denotes the probability of an empty cell in the 1D system with both neighbors empty, then one has that

For this 1D model, a pair approximation factorization is in fact exact, so the only approximation is in the factorization of P₇ in (C1). Then, from the equality J_1D→3D = J_3D→1D, one obtains 〈X_int〉 ≈ 0.3057 (versus the precise KMC value of 0.273) for 〈X_b〉 = 0.20. One also obtains 〈X_int〉 ≈ 0.1374 (versus the KMC value of 0.134) for 〈X_b〉 = 0.10, etc. Not surprisingly, one finds that the pair approximation is somewhat less accurate in predicting the strong enhancement of concentration in the pore interior relative to the bulk (at least for higher 〈X_b〉) compared to its success in predicting the weak enhancement near infinite walls in  Appendix A. There are other analytic strategies which could be employed, e.g., matching chemical potentials for the 1D and 3D systems, where the latter might be determined, e.g., from virial expansion. However, the pair approximation clearly captures the key feature of concentration enhancement inside the pore.

Precise direct assessment of concentration enhancement inside the pore can naturally also be achieved by Monte Carlo simulation of the coupled 3D and 1D systems. We have implemented such simulations and recover the previously reported values of 〈X_int〉 from KMC simulations of the reaction model.

In Sec. III A, we have described just the lowest-order equations in the coupled hierarchy of exact evolution equations for the stochastic reaction model, e.g.,

for 3 ≤n ≤ L − 2. As indicated in Sec. III A, separate equations are needed for cells at the end of the pore. For example, for n = 1, one has

and

where 〈E₀〉 = 1 − 〈A₀〉, 〈E₁A₂E₃〉 = 〈E₁A₂〉 = 〈A₂〉, and appropriate factorizations are implemented for probabilities of hopping involving the state of cells both inside and just outside the pore. An example of a next-highest-order equation in the hierarchy is

We have grouped terms for forward and reverse hopping events between pairs of sites corresponding to loss and gain of the configuration of interest. Since cells adjacent to A or B must be empty, including this feature means that grouped hopping terms include the same set of cells. For example, 〈A_nE_n+1E_n+2〉 = 〈E_n−1A_nE_n+1E_n+2〉 specifies the state of cells n − 1, n, n + 1, and n + 2, as does 〈E_n−1E_nA_n+1E_n+2〉.

Next, we comment further on MF-type factorization approximations which facilitate truncation of the hierarchy to yield a closed set of evolution equations. The site approximation ignoring all correlations sets

so, e.g., 〈A_nE_n+1E_n+1〉 ≈ 〈A_n〉〈E_n+1〉〈E_n+2〉, leading immediately to a closed set of equations for 〈A_n〉 and 〈 B_n〉. However, as noted previously, this approximation does not impose the basic constraint for models with NN exclusion that the concentration in any cell should be no higher than $12$⁠. The pair approximation sets

so, e.g., 〈A_nE_n+1E_n+2〉 ≈ 〈A_nE_n+1〉〈E_n+1E_n+2〉/〈E_n〉 = 〈A_n〉(1 − 〈X₁〉 − 〈X₂〉)/(1 − 〈X₁〉). This again leads to a closed set of equations for 〈A_n〉 and 〈B_n〉. Results from numerical analysis of these equations are shown in Figs. 4-6, and in the supplementary material. The triplet approximation sets

Thus, this approximation does not directly approximate any quantities (in the flux terms) in the lowest-order equations. However, in higher-order equations such as (D4), one must implement factorization, e.g., 〈A_n−1E_nE_n+1E_n+2〉 ≈ 〈A_n−1E_nE_n+1〉〈E_nE_n+1E_n+2〉/〈E_nE_n+1〉. This expression can be recast noting that 〈E_nE_n+1E_n+2〉 = 〈E_n+1E_n+2〉 − 〈A_nE_n+1E_n+2〉 − 〈B_nE_n+1E_n+2〉 and 〈E_n+1E_n+2〉 = 1 − 〈X_n+1〉 − 〈X_n+2〉. From numerical analysis of the equations for the triplet approximation, we find only minor improvement over the pair approximation. See the supplementary material Fig. S4. This further demonstrates the challenge of capturing strong non-equilibrium spatial correlations with MF-type approximations, and also highlights the success of the GH approach.

Finally, we discuss the evaluation of tracer diffusivity within the site and pair approximations. To this end, consider the generic form of the reaction-diffusion equations, and specifically the diffusion flux, away from the pore openings where 〈E_n〉 ≈ 〈E_int〉 = 1 − 〈X_int〉 is effectively constant. In the site approximation, factorizing 〈A_nE_n+1E_n+2〉 ≈ 〈E_int〉²〈A_n〉 and 〈E_n−1E_nA_n+1〉 ≈ 〈E_int〉²〈A_n+1〉 yields J_X^n>n+1 = − h〈E_int〉² ∇〈A_n+1〉. On the other hand, first utilizing exact identities and then factorizing corresponding to 〈A_nE_n+1E_n+2〉 = 〈A_n − E_n+2〉 ≈ 〈E_int〉〈A_n〉 and similarly for 〈E_n−1E_nA_n+1〉, yields

We adopt the last analysis which to some extent accounts for NN exclusion. Noting that this analysis applies for a standard counter-diffusion mode, it follows that D_tr(site) = h〈E_int〉 = h(1 − 〈X_int〉). In the pair approximation, factorizing 〈A_nE_n+1E_n+2〉 ≈ 〈E_nE_n+1〉〈A_n〉/〈E_n〉 ≈ (1 − 2〈X_int〉)〈A_n〉/(1 − 〈X_n〉) and similarly for 〈E_n−1E_nA_n+1〉 yields

Noting that this analysis applies for a standard counter-diffusion mode, it follows that D_tr(pair) = h(1 − 2〈X_int〉)/(1 − 〈X_n〉). Hence, these analyses provide a derivation of (7).

The (pure) diffusion equations, d/dt 〈X_n〉 = − ∇J_X^n>n+1, for the total concentration profile 〈X_n〉 within the pore are non-trivial due to the NN exclusion constraint. The non-trivial feature is the appearance of triplet probabilities in the expression for the diffusion flux, J_X^n>n+1 = h[〈X_nE_n+1E_n+2〉 − 〈E_n−1E_nX_n+1〉], for 3 ≤ n ≤ L − 2. As with reaction-diffusion equations, some modification is required for the end sites within the pore. For example, one has that

and

where 〈E₀〉 = 1 − 〈X₀〉, and appropriate factorizations are implemented for probabilities of hopping involving the state of cells both inside and just outside the pore.

Our interest in these equations is the behavior of the solutions in the equilibrium steady-state. We have argued in (4) that in the equilibrium state (but not for time evolution), the factorization of the pair approximation, e.g., 〈X_nE_n+1E_n+2〉 = 〈X_nE_n+1〉〈E_n+1E_n+1〉/〈E_n+1〉, becomes exact. This is a consequence of the Markov random field property of equilibrium lattice-gas models in any dimension with NN interactions.²⁹ It is applied here for the special case of a 1D lattice-gas model with NN exclusion. To clarify this issue, consider the conditional probabilities,

for cell n to be in state C given that cell n + 1 is in state D, cell n + 2 is in state F, etc. Then the spatial Markov property implies that 〈C_n|D_n+1F_n+1…〉 = 〈C_n|D_n+1〉, and in particular that 〈C_n|D_n+1F_n+1〉 = 〈C_n|D_n+1〉. The latter equality demonstrates that the factorization used in the pair approximation becomes exact.

Application of this factorization allows exact solution for steady-state 〈X_n〉 by solution of the resulting coupled set of equations given the values of 〈X₀〉 and P₅ in Sec. II B recover exactly the oscillations in total concentration within the pore, i.e., the concentration oscillation which would be seen in the coupled 1D pore + 3D extended fluid system.

As a final aside, we offer a simple test case for the validity of our strategy of capturing behavior in the pore for a coupled system with analysis just of the pore. Consider a coupled 1D pore + 1D extended fluid again with NN interactions. This just corresponds to an infinite 1D lattice-gas model with NN exclusion so the concentration should be constant, 〈X_b〉, everywhere in equilibrium. Refining the above equations for this 1D case (where 〈X₀〉 is replaced by 〈X_b〉/(1 − 〈X_b〉) and P₅ is replaced by (1 − 2〈X_b〉)/(1 − 〈X_b〉)), we find that the equations are consistent with a solution 〈X_n〉 = 〈X_b〉 for all n.