Dissociative adsorption onto a surface introduces dynamic correlations between neighboring sites not found in non-dissociative absorption. We study surface coverage dynamics where reversible dissociative adsorption of dimers occurs on a finite linear lattice. We derive analytic expressions for the equilibrium surface coverage as a function of the number of reactive sites, N, and the ratio of the adsorption and desorption rates. Using these results, we characterize the finite size effect on the equilibrium surface coverage. For comparable N’s, the finite size effect is significantly larger when N is even than when N is odd. Moreover, as N increases, the size effect decays more slowly in the even case than in the odd case. The finite-size effect becomes significant when adsorption and desorption rates are considerably different. These finite-size effects are related to the number of accessible configurations in a finite system where the odd-even dependence arises from the limited number of accessible configurations in the even case. We confirm our analytical results with kinetic Monte Carlo simulations. We also analyze the surface-diffusion case where adsorbed atoms can hop into neighboring sites. As expected, the odd-even dependence disappears because more configurations are accessible in the even case due to surface diffusion.
I. INTRODUCTION
Adsorption and desorption processes provide an essential mechanism for mass transport at fluid–solid interfaces.1,2 For example, heterogeneous catalysis3,4 relies on this mechanism to transport reactants from a gas onto a catalytic surface where chemical reactions occur and bring products back to the gas phase. Therefore, developing a correct description of adsorption and desorption processes is a crucial step in the computational modeling of gas–solid interfacial systems requiring appropriate modeling assumptions and careful analysis.5–8 While the inclusion of lateral interactions (also referred to as adsorbate–adsorbate interactions) is important for realistic modeling of the phenomenon, it makes the analytic investigation of the behavior of the resulting system intractable. In this paper, we consider a theoretical model of reversible dissociative adsorption based on Langmuir adsorption modeling and show that, even without lateral interactions, the phenomenon of reversible dissociative adsorption exhibits rich dynamics that requires detailed analysis.
The Langmuir adsorption model1,2,9 has served as the most influential theoretical model for reversible adsorption processes. Despite its simplicity, the model not only captures key molecular features but also gives analytical expressions for adsorption isotherms. One of the fundamental assumptions of the model is that there are no interactions between adsorbates on adjacent sites. This assumption implies another assumption, namely, that the occupancy of each site becomes uncorrelated in the infinite equilibrium system, from which one can derive the Langmuir isotherms for both non-dissociative (or one-site) and dissociative (or two-site) adsorption. However, we note that, for the dissociative adsorption case, the validity of the uncorrelated site occupancy assumption may not be guaranteed in a general situation (i.e., finite or nonequilibrium system) because adsorption and desorption events can lead to dynamic correlations between neighboring sites. These kinetically induced correlations present in the dissociative adsorption case make the analysis of system behaviors nontrivial. This is contrasted with the non-dissociative adsorption case, where the uncorrelated site occupancy assumption always holds and the system exhibits trivial exponential kinetics because the occupancy of each site can be modeled by an independent stochastic process under the assumption of no lateral interactions.
Theoretical investigations of kinetically induced correlations in a system undergoing two-site adsorption events date back to at least the 1939 study by Flory10 where intramolecular reactions on polymer chains were modeled by irreversible two-site adsorption events (i.e., without desorption events). This study has led to various analyses of emerging behavior in the random and cooperative sequential adsorption models.11 For the reversible dissociative adsorption case, a systematic analysis was performed in the context of deposition and evaporation of k-mers (k = 2, 3, …) on a linear lattice.12,13 Using an equivalent formulation based on a quantum-spin model, it was shown that the autocorrelation function for the number of adsorbed atoms exhibits a power-law decay in time (i.e., ∼t−1/2). To this end, in the context of the Goldstone theorem,14 a family of conservation laws was identified by dividing the lattice into k sublattices and considering the number of atoms adsorbed on each sublattice. While this approach provides an insightful explanation for the origin of the nontrivial power-law decay observed in the infinite-system limit, we note that the approach is not fully applicable to a finite system, particularly if the size of the lattice is not a multiple of k. Liu and Evans15 analyzed spatial correlations in one- or two-dimensional lattice systems undergoing reversible dissociative adsorption using a similar formulation with sublattices. They showed that the magnitude of nearest-neighbor and other short-range correlations decay like t−d/2 where d is the dimensionality of the lattice. This analysis for the dimer (O2) adsorption–desorption process was performed as part of a broader study of a CO oxidation model16 and strong spatial correlations appearing in some quasi-steady states during the evolution of surface coverage dynamics were used to explain why phenomenological kinetics (e.g., mean-field description) may fail to provide an adequate description of heterogeneous catalysis.17,18 We note that previous studies12,13,15 mostly considered the infinite-system limit and analyzed the power-law decay of dynamic correlations to demonstrate the intriguing nature of the reversible dissociative adsorption dynamics.
In this paper, we perform a systematic analysis of reversible dissociative adsorption occurring on a finite linear lattice. The main quantity of interest is the equilibrium surface coverage. We show that, contrary to the infinite system case, the effect of kinetically induced correlations can be seen in this static quantity because the occupancy of each site is not completely uncorrelated in a finite system. We also show that the finite system-size effect on the equilibrium surface coverage exhibits interesting behavior depending on whether the lattice has an odd or even number of sites and this can be explained by the number of accessible configurations that a finite system can have. Analyzing finite system-size effects and investigating their physical origins are crucial to understanding emerging behaviors in finite-sized physical systems (e.g., nano-engineered materials)19,20 or perform reliable simulation studies using finite-sized computational models (e.g., with periodic boundary conditions).21–24 In fact, considering that even high-quality single crystal surfaces have terraces that seldom exceed a size of a few hundred sites in one direction, we note that most real catalytic systems are nano-structured without any particular engineering.
We formulate a continuous-time Markov chain model25 for a linear lattice undergoing reversible dissociative adsorption and investigate the dynamics of this system using the following two methods. First, we develop an analytic approach based on the chemical master equation (CME). The CME is a set of first-order differential equations that describe the probabilistic time evolution of the system in state space.26,27 We derive analytical expressions for the equilibrium surface coverage as a function of system size and the ratio of the adsorption and desorption rates. Second, we perform lattice kinetic Monte Carlo (KMC) simulations28–30 to confirm the validity of our analytical results for the equilibrium surface coverage. KMC numerically solves the CME in the sense that the distribution of sample trajectories is the solution of the CME. Furthermore, we investigate the time evolution of the surface coverage dynamics to discuss how the equilibrium surface coverage is reached. In addition, using both CME and KMC approaches, we investigate the effect of surface diffusion, which is known to reduce spatial correlations.15,31
The rest of the paper is organized as follows. In Sec. II, we introduce our lattice system and reversible dissociative adsorption and derive the Langmuir isotherm. In Sec. III, we formulate a continuous-time Markov chain model and derive analytic results for the equilibrium surface coverage and the correlation coefficient for the occupancy of neighboring sites. Using these results as well as KMC simulations, we analyze finite system-size effects on equilibrium surface coverage. In Sec. IV, we analyze the effect of surface diffusion. In Sec. V, we conclude the paper with a summary and an outline for future work.
II. SYSTEM
We consider a theoretical model of reversible dissociative adsorption based on Langmuir adsorption modeling. The surface is represented as a linear lattice with N reactive sites, see Fig. 1(a). We assume that all reactive sites are identical and each site has two neighboring sites. In other words, the terminal sites are connected via a periodic boundary. One can consider this lattice system with a sufficiently large number of sites as an approximation to an infinite lattice. One may also view this linear lattice model with a small number of sites as a simple theoretical model for active sites adjacent to doped sites,32 see Fig. 1(b).
(a) An example of a linear lattice with N = 6 sites and a typical configuration with two occupied sites (marked with an X). (b) A linear lattice with a periodic boundary can be viewed as a simple theoretical model for active sites (in yellow) adjacent to doped sites (in black) on a surface of inert sites (in white). Two cases with N = 6 and 9 are shown.
(a) An example of a linear lattice with N = 6 sites and a typical configuration with two occupied sites (marked with an X). (b) A linear lattice with a periodic boundary can be viewed as a simple theoretical model for active sites (in yellow) adjacent to doped sites (in black) on a surface of inert sites (in white). Two cases with N = 6 and 9 are shown.
III. FINITE-SIZE EFFECTS: EVEN AND ODD N
In this section, we develop a continuous-time Markov chain description for the system described in Sec. II and investigate finite system-size effects on equilibrium surface coverage using both analytical and simulation approaches. The analytical approach is based on the master equation description, for which a complete characterization of all accessible configurations is a prerequisite. As briefly mentioned in the Introduction, the finite-size effect is related to the number of accessible configurations. We present a combinatorial argument to characterize all accessible configurations in Sec. III A. We then describe our master-equation-based approach in Sec. III B. We present analytic results for the equilibrium surface coverage and the correlation coefficient for the occupancy of neighboring sites in Secs. III C and III D, respectively. We finally present KMC simulation results in Sec. III E.
A. Accessible configurations
As described in Sec. II, we consider a periodic linear strip with N sites and label these sites from 1 to N. We assume that each site in the system is initially unoccupied. We define a configuration of the system by specifying whether each site is occupied (denoted by X) or unoccupied (denoted by ∅). We say that a configuration is accessible from another configuration if the former can be obtained from the latter via a sequence of reactions. For example, for N = 5, configuration {X∅∅X∅} is accessible from the initial configuration {∅∅∅∅∅} because the former is obtained from the latter via two adsorption events at sites 1-2 and 3-4 followed by a desorption event at sites 2-3. Using a combinatorial argument, we count the total number of accessible configurations, denoted by ntot. We also determine the number of accessible configurations with 2l occupied sites, denoted by , which will be used to derived in Sec. III C.
Since each site is either occupied or unoccupied, there are a total of 2N configurations. However, it is important to note that all configurations are not accessible to each other. This is because of properties that must continue to hold when the system undergoes a sequence of dissociative adsorption and associative adsorption events. We first note that the parity (i.e., whether odd or even) of the number of occupied sites does not change because the number of occupied sites increases (or decreases) by 2 via a dissociative adsorption (or associated desorption) event.
For the even case of N = 2m, there is an additional conserved quantity.12,13,15 To define this quantity, we consider the two alternating sublattices. The (+) sublattice only has sites with an odd site number, whereas the (−) sublattice only contains sites with an even site number. We denote the number of occupied sites in the (+) and (−) sublattices by N+ and N−, respectively. Since dissociative adsorption and associative desorption events occur at two neighboring sites, one of the sites belongs to the (+) sublattice and the other belongs to the (−) sublattice. As a result, the quantity N+ − N− is conserved when adsorption or desorption occurs. Hence, contrary to the odd case, all 2N−1 configurations with an even number of occupied sites are not accessible from the initial unoccupied configuration. For example, for N = 6, {XXX∅X∅} is not accessible from {∅∅∅∅∅∅} because the N+ − N− values of the former and latter configurations (2 and 0, respectively) are different.
B. Continuous-time Markov chain description
We construct a continuous-time Markov chain model for the system described in Sec. II by considering all accessible configurations and defining transition rates between each pair of configurations. If two configurations are obtained from each other by an adsorption or desorption event, the transition rates are set to ra and rd, respectively; otherwise, zero transition rates are set. To describe the time evolution of the continuous-time Markov chain model, one can use the chemical master equation (CME), which is a set of first-order differential equations whose solution gives the probability that the system is in a certain configuration at a certain time. However, the dimension of the CME, which is equal to ntot, grows quickly, see Eqs. (4) and (5). Moreover, even for small values of N, the values of ntot are rather large (e.g., ntot = 16 and 20 for N = 5 and 6, respectively), which makes it difficult to investigate the CME analytically. Instead of the standard approach which keeps track of all accessible configurations separately, we group configurations with the same characteristics into an aggregated state and write the CME for those aggregated states. Owing to the periodic boundary or the ring structure, configurations obtained via cyclic translation belong to the same aggregated state. We introduce a notation to denote an aggregated state containing the configuration α and all the other configurations that are reached from α via cyclic translation.
Transition diagrams for a five-site system. Panel (a) shows the standard approach, where 16 configurations accessible from the initial unoccupied configuration are considered separately. A double-sided arrow between two configurations indicates that one configuration is obtained from the other configuration via either adsorption or desorption and vice versa. The transition rate is ra (or rd) for the direction increasing (or decreasing) the number of occupied sites. Panel (b) shows our approach based on aggregated states. Note that the transition rates between aggregated states are multiples of ra and rd. See the main text for how the multiplicities are determined.
Transition diagrams for a five-site system. Panel (a) shows the standard approach, where 16 configurations accessible from the initial unoccupied configuration are considered separately. A double-sided arrow between two configurations indicates that one configuration is obtained from the other configuration via either adsorption or desorption and vice versa. The transition rate is ra (or rd) for the direction increasing (or decreasing) the number of occupied sites. Panel (b) shows our approach based on aggregated states. Note that the transition rates between aggregated states are multiples of ra and rd. See the main text for how the multiplicities are determined.
C. Equilibrium surface coverage
1. Analytic formulas
2. Derivation
3. Finite system-size effect
We now analyze the behavior of and . Figure 3 shows the curves of for several small values of N as a function of k. We first confirm that both the even and odd formulas give the same infinite-system limit, i.e., , which coincides with the Langmuir isotherm given in Eq. (3). This shows that, even without infinitely fast surface diffusion of X, each site in the infinite system becomes uncorrelated at equilibrium. In a finite system, however, reactive sites are not completely uncorrelated, causing finite system-size effects on the equilibrium surface coverage . We investigate this correlation in Sec. III D.
The equilibrium surface coverage of an N-site system is shown as a function of k = rd/ra for various values of N. The curves of are plotted for even values of N in panel (a) and for odd values of N in panel (b) using Eq. (11). To clearly show the finite system-size effect in , the curves with small values of N (solid lines) are compared with the infinite-limit case given in Eq. (3) (dashed lines).
The equilibrium surface coverage of an N-site system is shown as a function of k = rd/ra for various values of N. The curves of are plotted for even values of N in panel (a) and for odd values of N in panel (b) using Eq. (11). To clearly show the finite system-size effect in , the curves with small values of N (solid lines) are compared with the infinite-limit case given in Eq. (3) (dashed lines).
While both the even and odd cases converge to the same value of , they exhibit remarkably different convergence behaviors. Figure 3 shows that finite system-size effects are more significant and persist longer when N is even. In the odd case, we observe that for each value of N there is a range of k centered around k = 1 where the values are close to . The width of this region increases as N increases. In the even case, the discrepancy between and is significantly larger and convergence as N → ∞ is slower over the entire range of k.
Figure 4 shows a more detailed analysis of the finite system-size effect as measured by . As N increases, ɛ decreases like ɛ ∼ N−1 in the even case, whereas ɛ decreases exponentially [i.e., ɛ ∼ e−a(k)N] in the odd case. Interestingly, we find that the odd case, in Eq. (11) is equal to the Padé approximation33 of Eq. (3) around k = 1 of orders (m − 1, m), meaning that it is the best rational approximation around k = 1 up to a given order in the power series expansion. We also observe that the sign and magnitude of the finite system-size effect depend on the value of k. Figure 3 shows that if N is even then (or ) for k < 1 (or k > 1), while if N is odd then for all k ≠ 1. When the magnitudes of ra and rd are comparable (i.e., k is close to unity), finite system-size effects become less significant. In fact, when k is exactly equal to unity, both formulas match the infinite-limit value and thus there is no finite system size effect on .
Convergence behavior of the equilibrium surface coverage in the limit N → ∞. For even values [shown in panel (a)] and odd values [panel (b)] of N, the finite system-size effect is plotted. Panel (a): to show ɛ ∼ N−1 for the even case, the curves of ɛ are plotted in the log–log scale for various values of k and a straight line with slope −1 is also plotted for comparison. Panel (b): to show ɛ ∼ e−a(k)N for the odd case, the curves of ɛ are plotted in the semi-log scale for the same values of k as in panel (a).
Convergence behavior of the equilibrium surface coverage in the limit N → ∞. For even values [shown in panel (a)] and odd values [panel (b)] of N, the finite system-size effect is plotted. Panel (a): to show ɛ ∼ N−1 for the even case, the curves of ɛ are plotted in the log–log scale for various values of k and a straight line with slope −1 is also plotted for comparison. Panel (b): to show ɛ ∼ e−a(k)N for the odd case, the curves of ɛ are plotted in the semi-log scale for the same values of k as in panel (a).
4. Symmetry and limiting behaviors
Figure 3 also shows that the curve of has a reflection symmetry around . In other words, is satisfied, implying that the identical equilibrium surface coverage is obtained when switching the notions of X (occupied) and ∅ (unoccupied) and the values of ra and rd. However, this symmetry property does not hold for . In fact, if the fully occupied configuration is chosen as the initial configuration, the resulting equilibrium surface coverage is different from in Eq. (11) and given as .
D. Correlations between neighboring sites
Figure 5 shows the curves of as a function of k for small values of N. When N is an odd number, is an odd function in log k about k = 1. Except for N = 3, there is a range of k values centered around k = 1 where the values of are very close to zero, that becomes wider as N increases. As k → 0, whereas as k → ∞, . For an even number N > 2, as in the odd case, there is a range of k values around k = 1 where the magnitude of becomes smaller that increases as N increases. However, for the even case, is an even function in ln k about k = 1, is positive for all values of k, and has the minimum value at k = 1. As k → 0 or k → ∞, .
The correlation coefficient between two neighboring sites in an N-site system is shown as a function of k = rd/ra for various values of N. The curves of are plotted by colored solid lines for even values of N in panel (a) and for odd values of N in panel (b) using Eq. (19). The infinite-limit case (i.e., as N → ∞) is shown by black dashed lines.
The correlation coefficient between two neighboring sites in an N-site system is shown as a function of k = rd/ra for various values of N. The curves of are plotted by colored solid lines for even values of N in panel (a) and for odd values of N in panel (b) using Eq. (19). The infinite-limit case (i.e., as N → ∞) is shown by black dashed lines.
The behavior of the correlation coefficient for the even and odd cases explain why the finite system-size effect on the equilibrium surface coverage becomes more significant and persistent in the even case. In the odd case, correlations between neighboring sites are negligible in equilibrium in a neighborhood of k = 1, in contrast to the even case where they are not. Furthermore, the range of k with negligible correlation increases as N increases. We revisit this relation of the finite system-size effect and correlations between neighboring sites in Sec. IV, where we consider surface diffusion.
E. Time-transient behavior of θN(t)
By performing KMC simulations, we numerically validate our analytic results for the equilibrium surface coverage given in Eq. (11), and also observe the time-transient behavior of θN(t). For the setup of KMC simulations, see the supplementary material.
In Fig. 6, we show the KMC and CME results for a large value of k = 50. Panels (a) and (b) show the time profiles of θN(t) for small values of N up to 8. Since the transition matrix R can be explicitly given for these N values, θN(t) can be also obtained by numerically solving the CME (7). We first confirm the agreement between the KMC and CME results, which cross-validates both approaches. Due to the large value of k, we observe significant system-size effects on the long-time limit of θN(t) (i.e., ). In addition, these effects are more severe when N is even. For example, the result with N = 8 has larger system-size effects than N = 5. We also investigate how these effects develop as time increases. Early in the simulation θN(t) grows rapidly and its curves with different values of N coincide. As later times, however, curves with smaller values of N start to reach their equilibrium values and diverge from curves corresponding to larger values of N, reflecting a lower long-time limit value for smaller N. This behavior appears in both even and odd cases. Panels (c) and (d) show the time profiles of θN(t) for larger values of N, 9 ≤ N ≤ 16 as obtained by KMC. As expected, those curves converge to the values predicted by Eq. (11), which confirms the validity of these analytic results. The characteristic behaviors appearing in panels (a) and (b) are also observed. In particular, remarkably slow convergence of θN(t) to θ∞(t) is observed for even values of N. In addition, for smaller values of N, is smaller as a result of the earlier rollover of the θN(t) curves compared to larger values of N for both even and odd cases.
For a large value of k = 50 (ra = 1 and rd = 50), time-dependent surface coverages θN(t) obtained by KMC simulations (colored lines) are compared with the CME results (black solid lines): for 2 ≤ N ≤ 8, the even case in panel (a) and the odd case in panel (b); for 9 ≤ N ≤ 16, the even case in panel (c) and the odd case in panel (d). The curves of θN(t) computed by the CME are plotted in panels (a) and (b). For larger systems shown in panels (c) and (d), where CME results are not available, we plot the values of computed by Eq. (11). The infinite-system equilibrium coverage is also shown by the dashed line. Error bars for KMC simulations are not shown for visual clarity. The magnitude of error bars is comparable to that of fluctuations appearing in each curve.
For a large value of k = 50 (ra = 1 and rd = 50), time-dependent surface coverages θN(t) obtained by KMC simulations (colored lines) are compared with the CME results (black solid lines): for 2 ≤ N ≤ 8, the even case in panel (a) and the odd case in panel (b); for 9 ≤ N ≤ 16, the even case in panel (c) and the odd case in panel (d). The curves of θN(t) computed by the CME are plotted in panels (a) and (b). For larger systems shown in panels (c) and (d), where CME results are not available, we plot the values of computed by Eq. (11). The infinite-system equilibrium coverage is also shown by the dashed line. Error bars for KMC simulations are not shown for visual clarity. The magnitude of error bars is comparable to that of fluctuations appearing in each curve.
IV. EFFECT OF SURFACE DIFFUSION
Including surface diffusion is expected to reduce correlations between neighboring sites that are caused by reversible dissociative adsorption. Hence, while each surface diffusion event itself does not change the instantaneous value of surface coverage, the surface coverage dynamics is modified by surface diffusion. In this section, we investigate how surface diffusion affects the finite system-size effect on surface coverage dynamics.
Before discussing analytic results, we first present KMC simulation results to emphasize different behaviors of and θN(t) when surface diffusion is considered. Figure 7 shows the time profiles of θN(t) for small values of N up to 8. The value of k is set to 50 using ra = 1 and rd = 50 whereas the rate for surface diffusion is set to rdiff = 1. In contrast with the no-diffusion case shown in Fig. 6 the positions of the equilibrium surface coverage are in order (i.e., ). In other words, the significant finite system-size effect in the even case that appears in the no-diffusion case is absent. We also observe that, for each value of N, θN(t) reaches its equilibrium value faster due to surface diffusion.
For a lattice system undergoing reversible dissociative adsorption as well as surface diffusion, time-dependent surface coverages θN(t) obtained by KMC simulations (colored lines) are shown for small N values, 2 ≤ N ≤ 8. The value of k is set to 50 using ra = 1 and rd = 50, whereas the rate for surface diffusion is set to rdiff = 1. The black solid lines indicate the positions of the equilibrium surface coverage estimated by the analytic formulas in Eq. (23) for each N. The infinite-system equilibrium coverage is also shown by the dashed line. Note that the simulation results without surface diffusion are shown in Fig. 6.
For a lattice system undergoing reversible dissociative adsorption as well as surface diffusion, time-dependent surface coverages θN(t) obtained by KMC simulations (colored lines) are shown for small N values, 2 ≤ N ≤ 8. The value of k is set to 50 using ra = 1 and rd = 50, whereas the rate for surface diffusion is set to rdiff = 1. The black solid lines indicate the positions of the equilibrium surface coverage estimated by the analytic formulas in Eq. (23) for each N. The infinite-system equilibrium coverage is also shown by the dashed line. Note that the simulation results without surface diffusion are shown in Fig. 6.
In order to further investigate these behaviors, for each value of N (3 ≤ N ≤ 8), we compute the time profiles of θN(t) for different values of rdiff and compare them with the no-diffusion case (i.e., rdiff = 0) in Fig. 8. For each odd value of N, the equilibrium is the same for all values of rdiff, including zero. The main difference due to the rdiff value is that θN(t) reaches faster as rdiff increases. For the even case, the same observations are made for all nonzero values of rdiff. However, the no-diffusion case with rdiff = 0 is singular in the sense that its equilibrium value is different from that obtained from all nonzero values of rdiff.
For each system size N (3 ≤ N ≤ 8), the time profiles of the surface coverage θN(t) obtained from KMC simulations with k = 50 (ra = 1, rd = 50) and different values of the surface-diffusion rate rdiff are plotted. The infinite-system equilibrium coverage is also shown by the dashed line. Note that surface diffusion cannot be considered for N = 2 and thus this case is omitted.
For each system size N (3 ≤ N ≤ 8), the time profiles of the surface coverage θN(t) obtained from KMC simulations with k = 50 (ra = 1, rd = 50) and different values of the surface-diffusion rate rdiff are plotted. The infinite-system equilibrium coverage is also shown by the dashed line. Note that surface diffusion cannot be considered for N = 2 and thus this case is omitted.
Figure 9 shows the curves of and vs k for small even values N when surface diffusion is considered. We note that the odd case is exactly the same as the no-diffusion case shown in Figs. 3 and 5. Contrary to the no-diffusion case, for each even value of N, there is a range of k values where is much closer to . This explains why significant system-size effects observed in Fig. 5 for the even case do not appear here. In fact, for both even and odd cases, it is observed that decreases exponentially, i.e., ɛ ∼ e−a(k)N; the convergence plot is similar to Fig. 4(b) (see the supplementary material). We also notice that, in the plot of the correlation coefficient in Fig. 9(b), for each even value of N > 2 there is a corresponding range of k where is much closer to zero. This demonstrates the close relation between the finite system-size effect on the equilibrium surface coverage and correlations between neighboring sites.
In panel (a), the equilibrium surface coverage with surface diffusion is shown as a function of k for small even numbers N using Eq. (23). The infinite-system limit is also shown. In panel (b), the correlation coefficient with surface diffusion is shown for small even numbers N using Eqs. (19), (23), and (24). The infinite-system limit (i.e., as N → ∞) is also shown. The odd case is omitted because it is exactly same as the no-diffusion case, see Figs. 3(b) and 5(b).
In panel (a), the equilibrium surface coverage with surface diffusion is shown as a function of k for small even numbers N using Eq. (23). The infinite-system limit is also shown. In panel (b), the correlation coefficient with surface diffusion is shown for small even numbers N using Eqs. (19), (23), and (24). The infinite-system limit (i.e., as N → ∞) is also shown. The odd case is omitted because it is exactly same as the no-diffusion case, see Figs. 3(b) and 5(b).
V. CONCLUSION
We have considered the surface coverage dynamics where reversible dissociative adsorption occurs on an initially unoccupied linear lattice. Unlike the molecular (or non-dissociative) adsorption case, this system exhibits finite system-size effects caused by dynamic correlations between neighboring sites. We investigated this finite size effect on the equilibrium surface coverage and relate it to non-vanishing static site correlations introduced by reversible dissociative adsorption. We also investigated the effects of surface diffusion of adsorbed atoms, which reduces site correlations.
We modeled the equilibrium surface coverage and time-transient surface coverage θN(t) of a finite lattice with N reactive sites using the chemical master equation (CME) and kinetic Monte Carlo (KMC). We derived analytical expressions for and verified them numerically for the case without surface diffusion, see Eq. (11), and the case with surface diffusion, see Eq. (23). Without surface diffusion, finite system-size effects are significant for even N when the ratio k = rd/ra is much larger or smaller than unity. By comparing to the case with surface diffusion, the behavior observed for the even case in Eq. (11) was explained by an insufficient number of accessible configurations. We further related this behavior with the persistent positive correlation coefficient for even N, and demonstrated the close relation between the finite system-size effect on the equilibrium surface coverage and correlations between neighboring sites.
We draw the reader’s attention to the following points. First, in our study, the equilibrium was defined as the steady state that the system attains with a given initial configuration as opposed to the one defined via a grand canonical distribution of configurations. Our analysis relies on the fact that not every pair of configurations is mutually accessible via reversible dissociative adsorption and, as a result, configurations are partitioned into classes such that only configurations in the same class are accessible to each other. In the sense that the steady-state of the system depends on the initial state, the system is not ergodic. Second, as discussed in Ref. 15, the conservation of the quantity N+ − N− still holds in a 2D square lattice. Hence, a similar system-size effect is expected when a square lattice has an even number of sites in each direction. Third, in a real system, the range of k = rd/ra can be much wider than the range [10−2, 102] considered in our study. While the considered range roughly corresponds to a range of [−0.1, 0.1] eV for adsorption free energies at room temperature, a model with strong binding can easily have a value beyond this range.34 Hence, more significant finite-size effects can be expected. Although this system-size effect may not be significant when surface diffusion or other surface reactions are introduced, our study implies that caution should be exercised when lattice KMC modeling is used for a surface lattice system undergoing reversible dissociative adsorption of dimers.
This study has the following possible future directions. First, one can investigate non-periodic systems or two-dimensional lattice systems. Alternatively, since our model can be considered as a simple theoretical model for active sites adjacent to doped sites,32 one can also consider a system consisting of several strips where strip k has Nk sites and Nk follows, for example, a Poisson distribution. Second, one can also consider dissociative adsorption of heteronuclear diatomic molecules, e.g., NO.35 Third, as mentioned in Introduction, our findings will be useful for the development of a multiscale simulation method for a fluid–solid interfacial system, where KMC is coupled with a mesoscopic continuum method, for example, fluctuating hydrodynamics.36,37
SUPPLEMENTARY MATERIAL
KMC simulation setup; Convergence behavior of θ̄N in the presence of surface diffusion.
ACKNOWLEDGMENTS
C.K. thanks Dr. François Blanchette and Dr. Lei Yue (both at UC Merced) for helpful discussions on the combinatorial expressions for the number of accessible configurations. This work was supported in part by the National Science Foundation under Grant Nos. CHE-2213368 and DMS-1840265. This work was also supported in part by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Applied Mathematics Program under Contract No. DE-AC02-05CH11231.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
Enrique Mercado: Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (supporting); Software (lead); Visualization (lead); Writing – original draft (lead). Hyun Tae Jung: Formal analysis (supporting); Investigation (supporting); Software (supporting); Validation (lead); Visualization (lead). Changho Kim: Conceptualization (equal); Data curation (supporting); Formal analysis (lead); Funding acquisition (lead); Investigation (lead); Methodology (lead); Project administration (lead); Resources (lead); Software (lead); Supervision (lead); Validation (supporting); Visualization (supporting); Writing – original draft (lead); Writing – review & editing (lead). Alejandro L. Garcia: Conceptualization (equal); Formal analysis (supporting); Investigation (supporting); Writing – original draft (supporting); Writing – review & editing (supporting). Andy J. Nonaka: Funding acquisition (lead); Project administration (supporting); Software (supporting); Supervision (supporting); Writing – original draft (supporting); Writing – review & editing (supporting). John B. Bell: Conceptualization (equal); Funding acquisition (supporting); Investigation (supporting); Writing – original draft (supporting); Writing – review & editing (supporting).
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.