Surface Coverage Dynamics for Reversible Dissociative Adsorption on Finite Linear Lattices

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 catalysis 3,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][6][7][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 model 1,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 Flory 10 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 Evans 15 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 (O 2 ) adsorption-desorption process was performed as part of a broader study of a CO oxidation model 16 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 studies 12,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][22][23][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 model 25 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) simulations [28][29][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 Section II, we introduce our lattice system and reversible dissociative adsorption and derive the Langmuir isotherm.In Sections 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 Section IV, we analyze the effect of surface diffusion.In Section 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 Figure 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 Figure 1(b).
As a reaction model, we consider reversible dissociative adsorption of diatomic gas molecules X 2 between the gas phase and the lattice: ( A molecule X 2 in the gas phase dissociates into atoms that are adsorbed onto two neighboring empty sites (denoted as ∅∅) and vice versa.The rates of the forward reaction (dissociative adsorption) and the reverse reaction (associative desorption) are denoted as r a and r d , respectively.r a and r d have units of inverse time and are taken as constant.Note that the dependence on the partial pressure of X 2 is included in r a and no lateral interactions are assumed.In addition we will consider surface diffusion of the adsorbed atom X in Section IV, see Eq. (21).The main quantity of interest in this paper is the surface coverage of the system.The surface coverage θ N (t) is defined as the ratio of the mean number of occupied sites at time t to the number of reactive sites, N. We assume that the lattice is initially unoccupied, i.e. θ N (0) = 0. We define the equilibrium surface coverage θN as the long-time (or steady-state) limit of θ N (t), i.e. θN = lim t→∞ θ N (t).
Before closing this section, we derive the Langmuir isotherm, which is the equilibrium surface coverage in the infinite system, i.e. θ∞ = lim N→∞ θN .We let brackets denote the probability that a certain n-site cluster (n = 1, 2, 3) is found in the lattice system.Specifically, [X] and [∅] denote the probabilities that a site is occupied by X and unoccupied, respectively, whereas and [∅∅] denote the probabilities that a nearest-neighbor pair of sites is in states of XX, X∅, ∅X, and ∅∅, respectively.By the hierarchical rate equations, the time evolution of [X] is given as 15 Although Eq. ( 2) is exact, it cannot be solved as an initial value problem because the equation is not closed.The two-site cluster quantities, [∅∅] and [XX], are not determined by the onesite cluster quantities, [X] and [∅] = 1 − [X].However, assuming an infinite equilibrium system, one can derive the steady-state value of [X], which is θ∞ , from Eq. ( 2).Since we assume there are no lateral interactions, the occupancy of each site becomes uncorrelated in the thermodynamic equilibrium 9 and thus one has By combining these with the condition that d dt [X] becomes zero in the steady state, one finally obtains Note that Eq. ( 3) cannot be used to derive any finite-system or nonequilibrium results, such as θ N (t), θN , or θ ∞ (t), because the uncorrelated site occupancy assumption is not guaranteed to hold and

III. FINITE-SIZE EFFECTS: EVEN & ODD N
In this section, we develop a continuous-time Markov chain description for the system described in Section II and investigate finite system-size effects on equilibrium surface coverage θN 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 Section III A. We then describe our master-equation-based approach in Section III B. We present analytic results for the equilibrium surface coverage and the correlation coefficient for the occupancy of neighboring sites in Sections III C and III D, respectively.We finally present KMC simulation results in Section III E.

A. Accessible Configurations
As described in Section 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 n tot .We also determine the number of accessible configurations with 2l occupied sites, denoted by ñl , which will be used to derived θN in Section III C.
Since each site is either occupied or unoccupied, there are a total of 2 N 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 ) is given by N 2l where m j = m! j!(m− j)! denotes a binomial coefficient.The total number of accessible configurations, n tot = 2 N−1 , is then given by the sum of ñl : 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 2 N−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.
Combinatorial expressions for ñl and n tot in the even case are obtained as follows.Since we only consider configurations accessible from the initial unoccupied configuration, those configurations have N + − N − = 0. Thus, if a configuration has 2l occupied sites (i.e.N + + N − = 2l), we know N + = N − = l.We note that there are m l ways to arrange N + = l atoms in the (+) sublattice with m sites and the same expression holds for the (-) sublattice.Hence, we obtain Note that the conservative quantity, N + − N − , was identified using sublattices in previous studies 12,13,15 .However, the main focus in these works was to explain the power-law decay of timecorrelation functions observed in the infinite system limit and combinatorial arguments were not developed.

B. Continuous-Time Markov Chain Description
We construct a continuous-time Markov chain model for the system described in Section 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 r a and r d , 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 n tot , grows quickly, see Eqs. ( 4) and ( 5).Moreover, even for small values of N, the values of n tot are rather large (e.g., n tot = 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 C(α) to denote an aggregated state containing the configuration α and all the other configurations that are reached from α via cyclic translation.
We use a 5-site system as an example to explain how to construct the transition diagram and For each state i, we introduce l i so that each configuration in the state has 2l i X atoms.For the current example, we have (l 1 , l 2 , l 3 , l 4 ) = (0, 1, 1, 2).The transition rate from state i to state i ′ is given as a multiple of r a or r d when adsorption (for l i ′ = l i + 1) or desorption (for l i ′ = l i − 1) occurs, respectively.The multiplicity factor is determined by counting how many configurations in state i ′ can be obtained from each configuration in state i.Thus, the transition rates for 1 → 2, 2 → 4, and 3 → 4 are 5r a , 2r a , and r a , respectively, whereas the transition rates for where p i denotes the probability that the system is in state i.Equivalently, we have a matrix form: Once we have written down a CME (7) for an N-site system, its solution can be expressed using a matrix exponential: p(t) = e tR p(0).While obtaining an analytic expression of p(t) is not a trivial task even for small N, one can determine p(t) accurately using a numerical method (e.g.Runge-Kutta).The equilibrium probability distribution p * is then given as the long-time limit of p(t), i.e. p * = lim t→∞ p(t).Alternatively, the equilibrium probability distribution p * can be obtained as the unique invariant probability distribution satisfying Rp * = 0 and ∑ i p * i = 1.Note that p * is the (right) eigenvector associated with zero eigenvalue satisfying the probability normalization condition; uniqueness is guaranteed because the system is finite and irreducible as formulated 25 .
Unlike computing the matrix exponential e Rt , determining p * analytically is a feasible task for small N.For N = 5, we obtain In Appendix A, we provide the transition matrix R and the equilibrium distribution p * for a 6-site system, where six aggregated states are used.

Analytic Formulas
Using the solution p(t), one can compute the surface coverage θ N (t) by a weighted sum of the components of p(t), where each weight is given as the respective surface coverage that the corresponding configurations represent (i.e.2l i /N).Since p * = lim t→∞ p(t), the equilibrium coverage θN can be obtained as θN = lim t→∞ θ N (t).Alternatively, we can directly compute θN by a weighted sum of the components of p * .We thus have For N = 5, using the equilibrium distribution of p * from Eq. ( 8), we obtain since (l 1 , l 2 , l 3 , l 4 ) = (0, 1, 1, 2).Using a similar procedure, we obtain analytic expressions of θN for 2 ≤ N ≤ 8, see Appendix B. These expressions satisfy the following general formulas, depending on whether N is even or odd, In Section III E, we numerically confirm these formulas by performing KMC simulations for larger values of N.

Derivation
Before investigating various behaviors of θ2m and θ2m+1 in Eq. ( 11), we note that these analytic results can be actually derived by observing that detailed balance is satisfied in the continuoustime Markov chain system.For any pair of configurations α (with 2l atoms) and β (with 2l + 2 atoms) that are connected by a certain pair of adsorption and desorption events, their equilibrium probabilities q * α and q * β satisfy r a q * α = r d q * β , and equivalently, q * β /q * α = 1/k.This implies that the equilibrium probability of a configuration depends on the number of adsorbed atoms and is proportional to k −l if there are 2l atoms in the configuration.Hence, the equilibrium probability of the aggregated state i can be written as where n i is the number of configurations in the state i and c(k) is the normalization constant for Using Eqs. ( 9) and ( 12), we obtain Note that we rewrote the summation by using index l for possible values of l i (l = 0, 1, • • • , m for both N = 2m and N = 2m + 1) and ñl = ∑ i n i δ l,l i (i.e. total number of configurations with 2l occupied sites) with the Kronecker delta δ l,l ′ .Similarly, we rewrite the normalization condition as  11).To clearly show the finite system-size effect in θN , the curves with small values of N (solid lines) are compared with the infinite-limit case θ∞ given in Eq. (3) (dashed lines).

Finite System-Size Effect
We now analyze the behavior of θ2m and θ2m+1 .Figure 3 shows the curves of θN 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. lim m→∞ θ2m = lim m→∞ θ2m+1 = θ∞ , 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 θN .We investigate this correlation in Section III D.
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 Figure 4 shows a more detailed analysis of the finite system-size effect as measured by ε = | θN − θ∞ |.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, θ2m+1 in Eq. ( 11) is equal to the Padé approximation 33 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 θN > θ∞ (or θN < θ∞ ) for k < 1 (or k > 1), while if N is odd then θN < θ∞ for all k ̸ = 1.When the magnitudes of r a and r d 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 θ∞ = 1 2 and thus there is no finite system size effect on θN . 2 ).In other words, 1 − θ2m ( 1 k ) = θ2m (k) is satisfied, implying that the identical equilibrium surface coverage is obtained when switching the notions of X (occupied) and ∅ (unoccupied) and the values of r a and r d .However, this symmetry property does not hold for θ2m+1 .In fact, if the fully occupied configuration is chosen as the initial configuration, the resulting equilibrium surface coverage θ * 2m+1 is different from θ2m+1 in Eq. ( 11) and given as θ * 2m+1 (k) = 1 − θ2m+1 ( 1 k ).We also notice that θ2m exhibits correct limiting behaviors for both k → 0 and k → ∞,

Symmetry and Limiting Behaviors
whereas θ2m+1 shows the correct limiting behavior only for k → ∞; however, there is a significant system-size effect for k → 0, This system-size effect reflects the fact that for odd N there is always at least one empty site.Likewise, if the fully occupied configuration is chosen as the initial configuration, a correct limiting behavior is expected for k → 0 but a significant system-size effect is expected for k → ∞ because there is always at least one occupied site: Disparate behaviors shown in the odd and even cases can be related to different accessibility of configurations.As mentioned in Section III A, in the odd case, the set of all configurations accessible from the unoccupied configuration is disjoint from the set of all configurations accessible from the fully occupied configuration.In the even case, the unoccupied and full occupied configurations belong to the same set characterized by N + − N − = 0.

D. Correlations between Neighboring Sites
Using the equilibrium probability distribution of accessible configurations, we quantitatively investigate correlations between neighboring sites.To this end, we first define a random variable Z n for the occupancy of the nth site (n = 1, • • • , N), that is, Z n = 1 if the nth site is occupied and 0 otherwise.We then define the correlation coefficient ρN between Z 1 and Z 2 .Note that any pair of two neighboring sites gives the same result of ρN .As shown in Appendix C, we express ρN in terms of [X] and [XX]: where [X] = θN is given in Eq. ( 11) and Figure 5 shows the curves of ρN as a function of k for small values of N. When N is an odd number, ρN 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 ρN are very close to zero, that becomes wider as 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 ρN becomes smaller that increases as N increases.However, for the even case, ρN is an even function in ln k about k = 1, is positive for all values of k, and has the minimum value 1 N−1 at k = 1.As k → 0 or k → ∞, ρN → 2 N .The behavior of the correlation coefficient ρN 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 Section ??, 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 θN 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 Figure 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.θN ).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 θN = lim t→∞ θ N (t) 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 θN 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, θ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.

IV. EFFECT OF SURFACE DIFFUSION
We now consider the case where the lattice system undergoes not only reversible dissociative adsorption but also surface diffusion.In other words, we allow an absorbed X to hop into a neighboring site if the site is unoccupied as described by X∅ where the rate is denoted by r diff .
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 θN 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 r a = 1 and r d = 50 whereas the rate for surface diffusion is set to r diff = 1.In contrast with the no-diffusion case shown in Figure 6 the positions of the equilibrium surface coverage θN 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 θN faster due to surface diffusion.
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 r diff and compare them with the no-diffusion case (i.e.r diff = 0) in Figure 8.For each odd value of N, the equilibrium θN is the same for all values of r diff , including zero.The main difference due to the r diff value is that θ N (t) reaches θN faster as r diff increases.For the even case, the same observations are made for all nonzero values of r diff .
However, the no-diffusion case with r diff = 0 is singular in the sense that its equilibrium value is different from that obtained from all nonzero values of r diff .
The singular behavior of the no-diffusion case for an even value of N results from an insufficient number of accessible configurations, as discussed earlier.For an odd value of N, the total number of configurations accessible from the initially unoccupied state is equal to n tot = 2 N−1 whether surface diffusion is included or not.On the contrary, for an even value of N, if surface diffusion is not included, some configurations become inaccessible due to the conserved quantity N + − N − = 0 and n tot = N N/2 , see Eq. ( 5).If surface diffusion is included, however, all configurations with an even number of occupied sites become accessible.Hence, for both N = 2m and N = 2m + 1, we Based on this observation, we derive analytic expressions for θN for the surface-diffusion case.
Under the assumption that detailed balance also holds in the presence of surface diffusion, Eq. ( 15) is valid.By substituting Eq. ( 22) to Eq. ( 15), we obtain Similarly, we obtain analytic expressions of the correlation coefficient ρN using Eqs.( 19), (23)   and See Appendix C for the derivation of Eq. ( 24).We note that, in the odd case, the analytic expressions of θN and ρN are exactly the same as the no-diffusion case, see Eqs. ( 11) and ( 20), and this is why the singular behavior does not appear in the odd case in Figure 8.
Figure 9 shows the curves of θN and ρN versus 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 Figures 3 and 5. Contrary to the no-diffusion case, for each even value of N, there is a range of k values where θN is much closer to θ∞ .This explains why significant system-size effects observed in Figure 5 for the even case do not appear here.In fact, for both even and odd cases, it is observed that ε = θN − θ∞ decreases exponentially, i.e. ε ∼ e −a(k)N ; the convergence plot FIG. 9.In panel (a), the equilibrium surface coverage θN 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 ρN with surface diffusion is shown for small even numbers N using Eqs.( 19), ( 23), and ( 24).The infinite-system limit (i.e.ρ∞ → 0 as N → ∞) is also shown.The odd case is omitted because it is exactly same as the no-diffusion case, see Figure 3 is similar to Figure 4(b) (see the Supplementary Material).We also notice that, in the plot of the correlation coefficient ρN in Figure 9(b), for each even value of N > 2 there is a corresponding range of k where ρN 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.

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 θN 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 θN 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 = r d /r a is much larger or smaller than unity.By comparing to the case with surface diffusion, the behavior observed for the even case θ2m in Eq. ( 11) was explained by an insufficient number of accessible configurations.We further related this behavior with the persistent positive correlation coefficient ρN 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 = r d /r a can be much wider than the range [10 −2 , 10 2 ] 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 N k sites and N k 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 .

FIG. 1 .
FIG. 1.(a) An example of a linear lattice with N = 6 sites and a typical configuration with 2 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.

(
or decreases) by 2 via a dissociative adsorption (or associated desorption) event.For the odd case of N = 2m + 1, this parity can completely characterize the two groups of configurations within which all configurations are accessible from each other.The 2 N−1 configurations with an even number of occupied sites are accessible from the completely unoccupied configuration {∅∅ • • • ∅}, whereas the other 2 N−1 configurations with an odd number of occupied sites are accessible from the completely occupied configuration {XX • • • X}.Hence, we will consider only the former group of configurations in the continuous-time Markov chain description for the odd case.One can easily see that the number of configurations with 2l occupied sites

FIG. 2 .
FIG. 2. Transition diagrams for a 5-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 r a (or r d ) 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 r a and r d .See the main text for how the multiplicities are determined.

FIG. 3 .
FIG. 3. The equilibrium surface coverage θN of an N-site system is shown as a function of k = r d /r a for various values of N. The curves of θN 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 θN , the curves with small values

FIG. 4 .
FIG. 4. Convergence behavior of the equilibrium surface coverage θN in the limit N → ∞.For even values (shown in panel (a)) and odd values (panel (b)) of N, the finite system-size effect ε(N) = | θN − θ∞ | 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).

Figure 3
Figure 3 also shows that the curve of θ2m has a reflection symmetry around (k = 1, θ = 12 ).In other words, 1 − θ2m ( 1 k ) = θ2m (k) is satisfied, implying that the identical equilibrium surface coverage is obtained when switching the notions of X (occupied) and ∅ (unoccupied) and the

FIG. 5 .
FIG. 5.The correlation coefficient ρN between two neighboring sites in an N-site system is shown as a function of k = r d /r a for various values of N. The curves of ρN 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.ρN → 0 as N → ∞) is shown by black dashed lines.

FIG. 6 .
FIG. 6.For a large value of k = 50 (r a = 1 and r d = 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 θN 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.

FIG. 7 .Figure 6
FIG. 7.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 r a = 1 and r d = 50, whereas the rate for surface diffusion is set to r diff = 1.The black solid lines indicate the positions of the equilibrium surface coverage θN 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 Figure6.

FIG. 8 .
FIG. 8.For each system size N (3 ≤ N ≤ 8), the time profiles of the surface coverage θ N (t) obtained from KMC simulations with k = 50 (r a = 1, r d = 50) and different values of the surface-diffusion rate r diff 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.