4272
Langmuir 1996, 12, 4272-4280
Simulation of Model Heterogeneous Surfaces in the Presence of Correlation Alessandra Adrover, Massimiliano Giona,* and Manuela Giustiniani Centro Interuniversitario sui Sistemi Disordinati e sui Frattali nell’Ingegneria Chimica, c/o Dipartimento di Ingegneria Chimica, Universita´ di Roma “La Sapienza” via Eudossiana 18, 00184 Roma, Italy Received January 16, 1996. In Final Form: May 14, 1996X We propose a general method to generate correlated energy landscapes for model heterogeneous surfaces. It is based on the representation of an energy landscape as a superposition of exponentially correlated Gaussian processes. This method is the extension of a previous work of Giona and Adrover to reconstruct binary lattice models of porous structures to a random energy field with values between 0 and ∞. The application to dual bond-site models is also developed.
1. Introduction The spatial distribution of adsorption energies influences both equilibrium and kinetic properties. The correlation properties of the energy landscape play a significant role especially in connection with surface diffusion.1,2 In the description of heterogeneous energy landscapes, two main archetypes are usually adopted:3 random topographies (in which the adsorption energies in two different sites are uncorrelated with each other) and patch models (in which the energies within a given patch are highly correlated). A way to introduce correlations in the generation of model heterogeneous surfaces has been recently analyzed by Riccardo et al.4 In the article of Riccardo et al. attention is focused on the correlation between nearest neighboring sites. The energy landscape is described by means of the probability density function f(Ei,Ei+1) associated with the site energy Ei and with the nearest neighboring site energy Ei+1. In the model of Riccardo et al. the presence of correlations between nearest neighboring sites induces correlations between two arbitrary sites which can be calculated in terms of the joint distribution of nearest neighboring sites. The correlation degree between two (or more) adjacent sites is particularly important in the study of dimer and polymer adsorption. The generation of an energy landscape with a prescribed correlation between nearest neighboring sites is achieved through a Monte Carlo relaxation algorithm, and the analysis is also extended to the dual site-bond model.5-8 Another interesting approach in the simulation of correlated energy landscapes and models of porous media has been proposed by Mayagoitia et al.9 * To whom correspondence should be addressed at Universita´ di Roma “La Sapienza”. E-mail:
[email protected]. Author is also affiliated with Dipartimento di Ingegneria Chimica, Universita´ di Cagliari, Piazza d’Armi, 09123 Cagliari, Italy. X Abstract published in Advance ACS Abstracts, July 15, 1996. (1) Riccardo, J. L.; Pereyra, V.; Rezzano, J. L.; Rodriguez Saa´, D. A.; Zgrablich, G. Surf. Sci. 1988, 204, 289. (2) Riccardo, J. L.; Chade, M. A.; Pereyra, V. D.; Zgrablich, G. Langmuir 1992, 8, 1518. (3) Rudzinski, W.; Everett, D. H. Adsorption of Gases on Heterogeneous Surfaces; Academic Press: New York, 1992. (4) Riccardo, J. L.; Steele, W. A.; Ramirez Cuesta, A. J.; Zgrablich, G. Langmuir, submitted for publication. (5) Mayagoitia, V.; Rojas, F.; Riccardo, J. L.; Pereyra, V. D.; Zgrablich, G. Phys. Rev. B 1990, 41, 7150. (6) Mayagoitia, V.; Rojas, F.; Pereyra, V. D.; Zgrablich, G. Surf. Sci. 1989, 221, 395. (7) Faccio, R. J.; Vidales, A. M.; Zgrablich, G.; Zhdanov, V. P. Langmuir 1993, 9, 2499. (8) Riccardo, J. L.; Pereyra, V.; Zgrablich, G. Langmuir 1993, 9, 2730.
S0743-7463(96)00043-1 CCC: $12.00
In the present article, we analyze the generation of a model heterogeneous surface by focusing attention on the entire spatial behavior of the energy correlation function and not only on nearest neighboring correlations. The approach followed is similar to that developed in the characterization of porous media to reconstruct lattice models with assigned statistical features,10-12 and the mathematical apparatus follows the analysis developed by Giona and Adrover13 for porous structures. The main difference between the reconstruction of porous media and the generation of model heterogeneous surfaces with prescribed probability density function and correlation properties is that in the first case the state of each site is discrete (either 0 if the site belongs to the porous matrix or 1 if the site belongs to the pore network), while in the latter the adsorption energies are defined on a continuum (between 0 and ∞). In the present analysis, a model heterogeneous surface is characterized by means of a probability density function (pdf) g(E) for the adsorption energies and by the secondorder correlation function C2E(x). This article is organized as follows. First, we analyze the generation of a model heterogeneous surface in the presence of correlation by considering site energies exclusively. Generation is developed by means of a superposition of exponentially decaying Gaussian processes. Some examples of application of this method are presented and the influence of correlation on adsorption in the presence of surface diffusion is briefly discussed. Finally, we develop the application of this method to the dual site-bond model in order to generate a distribution of site and bond energies with prescribed correlation properties. 2. Statement of the Problem The aim of this article is to provide a general method for the simulation of model heterogeneous surfaces with a given spatial energy distribution and a given energy density function g(E). Let us characterize the spatial energy distribution by means of the normalized spatial correlation function C2E(x) (9) Mayagoitia, V.; Rojas, F.; Kornhauser, I.; Perez-Aguilar, H. Proceedings of the Second ISSHAC, Zakopane, Poland, September 4-10, 1996, p 95. (10) Joshi, M. Y. A class of stochastic models for porous media. Ph.D. Dissertation, University of Kansas, Lawrence, KS, 1974. (11) Quiblier, J. A. J. Colloid Interface Sci. 1984, 98, 84. (12) Adler, P. M.; Jacquin, C. G.; Quiblier, J. A. Int. J. Multiphase Flow 1990, 16, 691. (13) Giona, M.; Adrover, A. AIChE J., in press.
© 1996 American Chemical Society
Simulation of Model Heterogeneous Surfaces
C2E(x) )
〈(E(y) - 〈E〉)(E(y+x) - 〈E〉)〉 〈E2〉 - 〈E〉2
Langmuir, Vol. 12, No. 17, 1996 4273
(1)
where x is the lattice coordinate, E(x) the adsorption energy field, and 〈E〉 the average energy
∫0∞Eg(E) dE
〈E〉 )
(2)
The average in eq 1 is obtained by summing over all the lattice sites labeled with the coordinate y. Under the assumption of isotropy, C2E(x) is solely a function of x ) |x|, i.e. of the distance between sites y and y + x. We assume, unless otherwise stated, that the energy distribution is isotropic. The statistical characterization of the spatial energy distribution is limited to the second-order correlation function C2E(x), so that the lattice model to be generated must have a given density function g(E) and a given correlation function C2E(x). A simple way to obtain this is to consider a continuous correlated random process Y(x)
Y(x) )
a(r′)ξ(r′+x) ∑ r′
(3)
a(r′)a(r′+x)/∑a (r′) ∑ r′ r′ 2
E(x) ) G(Y(x))
(5)
depending on the (Gaussian) distribution function FY of Y and on the distribution function FE(E) ) ∫E0 g(e) de of E. For each point x, the corresponding value of the adsorption energy E(x) can be obtained by solving the probability balance equation
∫-∞Y exp(-y2/2) dy ) FE(E(x)) )
1 x2π
∫0Eg() d
(6)
The nonlinear filter G is therefore given by G(Y) ) F-1 E (FY(Y)). For example, by choosing an exponential form for the energy density function,
g(E) ) β exp(-βE), 〈E〉 ) 1/β
(7)
∫-∞∞dy1 ∫-∞∞dy2 ×
[
]
1 ln[1 - FY(Y(x))] β
〈E2〉 - 〈E〉2
p(y1,y2) (9)
where
p(y1,y2) )
[
2
]
2
y1 + y2 - 2C2Y y1y2 1 exp (10) 2 1/2 2π(1 - C2Y) 2(1 - C22Y) is a bivariate Gaussian density function for the twodimensional Gaussian random variable (y1, y2) possessing correlation C2Y. By choosing the exponential energy density function eq 7, the functional relationship between C2E(x) and C2Y(x) attains the form
C2E )
∫-∞∞dy1 ∫-∞∞dy2 ×
[1 + ln(1 - FY(y1))][1 + ln(1 - FY(y2))]p(y1,y2) (11) and is independent of the value of the parameter β (see Figure 1). In many cases, the interaction of gases with solid surfaces is represented by energy density functions possessing several local maxima. Figure 1 shows the behavior of C2E(x) vs C2Y(x) for an energy density function given by the superposition of two truncated Gaussian pdf’s
g(E) ) 1
( (
x2πσ2
A1 exp -
)
(E - R1)2 2σ2
(
+ A2 exp -
))
(E - R2)2 2σ2
and
E ∈ [0, ∞)
the nonlinear filter G(Y(x)) attains the form
G(Y(x)) ) -
C2E )
(4)
The transformation from the Y process to the E(x) process10-12 representing the model heterogeneous surface is given by the nonlinear filter G
FY(Y(x)) )
only further condition to be imposed on E(x) is to admit the assigned spatial correlation function C2E(x). The correlation function C2E(x) is related to the corresponding correlation function C2Y(x) through the relation
(G(y1) - 〈E〉)(G(y2) - 〈E〉)
where ξ(x) is a normalized Gaussian random variable (with zero mean and unit variance), a(r) is the kernel of the linear filter, and the summation is extended over the values r′ belonging to the lattice. Since Y is a linear superposition of Gaussian variables, Y remains Gaussian with zero mean but exhibits a normalized correlation function14 C2Y(x) given by
C2Y(x) )
Figure 1. C2E(x) vs C2Y(x) for an exponential (a) and a twopeak (b) energy density function g(E). The dotted line is the curve C2E(x) ) C2Y(x).
(8)
Equation 8 ensures statistically that the E(x) process admits the assigned density function g(E), so that the (14) Doob, J. L. Stochastic Processes; John Wiley: New York, 1953.
where Ai ) [1 + erf(Ri/σx2) (i ) 1, 2), with R1 ) 4, R2 ) 9 au and σ ) 1. By applying the probability balance eq 6, it is possible to obtain the expression for the nonlinear filter G(Y) for the assigned energy density function. Therefore, for a given value of C2E(x), the corresponding value of C2Y(x) can be evaluated through eqs 9 and 10.
4274 Langmuir, Vol. 12, No. 17, 1996
Adrover et al.
The C2E-C2Y curve can be regarded as the calibration curve that relates the value of the correlation of the E(x) process possessing the pdf g(E) to the corresponding correlation value of the Gaussian Y(x) process. From the knowledge of C2Y(x), the coefficients of the linear filter {a(r)} are then calculated starting from eq 4, e.g. by means of optimization methods. However, we shall see in the next section that a suitable decomposition of Y(x) in terms of exponentially correlated basis processes considerably simplifies the evaluation of the coefficients a(r′).
Y(x) )
Y(x) )
∫a(r′) ξ(r′+x) dr′
and
C2Y(x) )
∫a(r′) a(r′+x) dr′/∫a2(r′) dr′
(12)
Unless otherwise specified, we will always consider the continuous formalism. 3.1. One-Dimensional Case. Let us first consider the one-dimensional case. As discussed in section 2, the generation of a model heterogeneous surface reduces to the generation of a one-dimensional Gaussian process with specified correlation function obtained from C2E(x) through eqs 9 and 10. Let x ∈(-∞, ∞) and let us define as basis processes {Y(x,λ)} the Gaussian processes defined as
Y(x,λ) )
∫0
∞
a(x′,λ)ξλ(x+x′) dx′
and
a(x,λ) ) x2λ exp(-λx)
(13)
Nloc
xRiν(x,λi) ∑ i)1
pˆ (λ) Y(x,λ) dλ +
Nloc
C2Y(x) )
∫0∞πˆ (λ)e-λx dλ + ∑Rie-λ x i
C2Y(x,λ) ) 〈Y(x′,λ) Y(x′+x,λ)〉 ) e-λx
(14)
i.e. it exhibits an exponential decay in x. In this case (d )1) 1/λ represents a correlation length (expressed in lattice units). A generic correlated process Y(x) can be regarded as the superposition of a continuous and a discrete spectrum of basis processes {Y(x,λ)}
(16)
i)1
and the weight function π(λ), encompassing both the continuous and the discrete part, is given by Nloc
π(λ) ) πˆ (λ) +
Riδ(λ - λi), ∑ i)1
πˆ (λ) ) pˆ 2(λ)
(17)
Therefore, under the condition that C2Y g 0, it follows that π(λ) is the inverse Laplace transform of the correlation function C2Y. It is important to observe that C2Y(x) is defined for real x, while the Laplace inversion requires a Laplace transform on the complex plane. We therefore need an analytic (p) (x) of C2Y(x) valid for all complex x continuation C2Y (whose restriction to real x coincides with C2Y(x)) such that
π(λ) )
1 2πi
+i∞ (p) C2Y(x)eλx dx ∫xx-i∞ o
o
(18)
The analytic continuation can be achieved by considering rational approximations (e.g. Pade´ approximants) or by means of other methods (e.g. orthogonal polynomial expansion). This topic is discussed in detail by Giona and Adrover.13 3.2. d-Dimensional Isotropic Structures. For ddimensional isotropic model surfaces (d > 1) it is convenient to define the basis processes {Y(x,λ)} in a slightly different way than in section 3.1, by convoluting a Gaussian uncorrelated process with a Gaussian kernel,
Y(x,λ) )
∫E a(u,λ) ξλ(u+x) du ) d/4 (4λπ) ∫E e-2λu ξλ(u+x) du d
2
d
(19)
where Ed is the Euclidean d-dimensional space, Ed ) {x| - ∞ < xi < ∞ (i ) 1, ..., d)}. In the d-dimensional case, 2/xλ represents a correlation length (expressed in lattice units). The correlation function for the basis process given by eq 19 is Gaussian: 2
where ξλ(x) is a Gaussian uncorrelated process with zero mean and unit variance. The correlation function C2Y(x,λ) of Y(x,λ) is given by
(15)
where Nloc is the number of discrete (localized) exponential contributions, {Y(x,λi)}, and the basis processes {Y(x,λ)} are uncorrelated to each other. Starting from its definition, eq 15, the process Y(x) admits the correlation function C2Y(x)
3. Closed-Form Solution with Linear Filters In the previous section we have seen that the problem of generating a model hereogeneous surface (with a given energy density function and a given spatial energy correlation function) reduces to the solution of the nonlinear functional eq 4. The approach followed to solve this equation is customarily based on optimization methods. In this section, we show that this problem can be recast in the form of a linear functional equation, analogous to a Laplace transform, which can be solved in closed form in many cases. This approach is based on the choice of a suitable set of basis processes with specified decay in the correlation function. The choice of the basis processes depends on the domain of definition. For this reason, different cases are discussed separately. In the theoretical development it is useful to recast the previous equations in integral form by taking the continuous limit of eqs 3 and 4. The continuous counterparts of these equations are
∫0
∞
C2Y(x,λ) ) e-λx
(20)
As discussed above, by considering the linear superposition of basis processes, eq 15, and putting z ) x2 and C2Y(x) ) C ˜ 2Y(z), it follows that
C ˜ 2Y(z) )
∫0 πˆ (λ)e ∞
Nloc -λz
dλ +
Rie-λ z ) ∫0 π(λ)e-λz dλ ∑ i)1 i
∞
(21)
Therefore, as in the one-dimensional case, eq 21 expresses the property that the weight function π(λ) can be evaluated
Simulation of Model Heterogeneous Surfaces
Figure 2. Comparison of C2E(x) ) exp(-λx) and the correlation function of the one dimensional E(x) process generated by considering Y(x) ) Y(x,λ), eq 13, for two different values of λ. The energy distribution function is exponential, eq 7.
Langmuir, Vol. 12, No. 17, 1996 4275
Figure 4. Comparison of C2E(z)x2) ) exp(-λz) and the correlation function of the two-dimensional energy landscape generated by considering Y(x) ) Y(x,λ), eq 19, for two different values of λ. The energy distribution function is exponential, eq 7.
Figure 5. Contour plot of the two-dimensional energy field E(x) with an exponential density function g(E), eq 7, (β ) 0.1) and a Gaussian spatial correlation function CE(x) ) exp(-λx2). (a) λ ) 0.1; (b) λ ) 0.01. Simulations were performed on a 200 × 200 square lattice.
Figure 3. One-dimensional energy field E(x) with an exponential density function g(E), eq 7, (β ) 0.1) and an exponential spatial correlation function CE(x) ) exp(-λx): (a, top) λ ) 0.1; (b, bottom) λ ) 0.01.
as the inverse Laplace transform of the analytic continu(p) ation on the complex plane C ˜ 2Y (z) of C ˜ 2Y(z). 4. Examples and Applications Let us consider one-dimensional model surfaces. As a consequence of the method adopted, it is quite easy to generate exponentially correlated processes, C2E(x) ) exp(-λx). When the energies are distributed exponentially, figure 1 shows that for positive correlations (C2E g 0) we have C2E = C2Y. Therefore an acceptable approximation for C2E(x) ) exp(-λx) is to choose Y(x) ) Y(x, λ), i.e. a basis process with the decay constant equal to λ. Figure 2 compares the assumed exponential decay in C2E(x) and the correlation function of the one-dimensional E process generated by means of the basis process Y(x, λ) for two different values of λ. The energy landscape for two different values of λ is shown in Figure 3. In a similar way, it is quite easy to generate twodimensional model surfaces which show Gaussian correlation properties, eq 20.
Figure 4 compares the assumed exponential decay in C2E(z)x2) and the correlation function of the twodimensional energy landscape generated by means of the basis process Y(x,λ), eq 19, for two different values of λ. The corresponding energy contour plots are shown in Figure 5. As the correlation length Lc ) 2/xλ increases, the transition from random to patchwise topographies becomes evident. 4.1. Closed-Form Solutions. Equations 16-18 and 21 state that the weight function π(λ) is the inverse Laplace transform of the correlation function C2Y. The entire apparatus of Laplace transform theory can therefore be applied to obtain π(λ) in closed form in many situations. As an example, let us consider a one-dimensional process Y(x) with a power-law decay in the normalized correlation function
C2Y(x) )
aυ (a + x)υ
(22)
where a is a positive constant and ν > 0. Such a process can be generated in a straightforward way by means of eq 15 because the inverse Laplace transform of C2Y(x) is analytic, and the weight function π(λ) attains the form
π(λ) ) πˆ (λ) )
aνλν-1 exp(-aλ) Γ(ν)
(23)
where Γ(υ) is the Gamma function of argument υ. Figure 6 compares eq 22 and the correlation function of the one-dimensional process Y(x) obtained by implementing eq 15 with eq 23. The E(x) process itself exhibits
4276 Langmuir, Vol. 12, No. 17, 1996
Adrover et al.
where P ˆ k(y) are the modified orthonormal Legendre polynomials of order k on [0, 1], P ˆ k(y) ) (2k + 1)1/2Pk(2y - 1), where Pk(x) are the Legendre polynomials, orthogonal in x ∈[-1, 1]. In this way, the problem is recast into the calculation of the expansion coefficients ck of a generalized Fourier series.
Figure 6. Comparison of the correlation function C2Y(x) of the stochastic one-dimensional process Y(x) generated by means of eqs 15 and 23 and the predicted power law decay, eq 22 (a ) 0.5, υ ) 1 ). Dots represent the behavior of the correlation function C2E(x) of the energy field E(x) obtained from Y(x) by means of the nonlinear filter, eq 8, in the case of an exponential g(E), eq 7. The simulation was performed on a one-dimensional lattice with 104 sites.
Figure 7. Comparison of the energy density function g(E) (dots) of the one-dimensional energy field E(x) obtained by applying the nonlinear filter, eq 8, to the Y(x) process and the expected exponential decay, eq 7 (continuous line).
a power-law behavior of the spatial correlation function C2E(x) (see Figure 6) because, for the assigned exponential decay of g(E), the relationship between C2E and C2Y, eq 11, can be approximated by C2E(x) = C2Y(x) for C2Y g 0, independently of the value of the parameter β (see Figure 1). Figure 7 compares the energy probability density function of the E(x) process obtained by appying the nonlinear filter, eq 8, to Y(x) and the expected exponential decay, eq 7. This example clearly shows that it is possible to generate correlated Gaussian processes with long tail correlation effects, avoiding any numerical problems. 4.2. Numerical Analysis. In all the cases for which closed-form expressions for the weight function π(λ) are not available, eqs 16-18 should be solved numerically. In order to apply eq 18, we need the analytic continuation of the correlation function on the complex plane. Orthogonal polynomial expansion can be useful to avoid the numerical problems associated with the Laplace inversion for the evaluation of πˆ (λ). This approach is discussed in detail by Giona and Adrover13 by considering some examples from the reconstruction of porous solids. It is convenient to define the new variable y ) 1/(z + 1) so that the positive real line z ∈[0, ∞) is mapped into y ∈[0, 1], and ˜ *2Y(y) ) C ˜ 2Y(1/y - 1). In this way C ˜ *2Y(y) can be to put C easily approximated by means of an orthonormal polynomial expansion N
C ˜ *2Y(y) )
∑ ckPˆ k(y)
n)0
(24)
5. Influence of Correlation on Adsorption Isotherms The effect of spatial correlations on adsorption isotherms has been extensively studied. For a review see Rudzinski and Everett,3 Jaroniec and Madey,15 and Steele.16 The influence of correlation on the expression of the virial expansion coefficients has been studied by Ripa and Zgrablich,17 and the functional form of the resulting adsorption isotherms on carbonaceous materials by Rudzinski et al.18 Lattice theories of adsorption on nonuniform surfaces have been proposed by Tovbin.19 The presence of correlations in the energy landscape of a solid surface influences the molecular processes evolving on it. A systematic study of the various cases (adsorption, surface diffusion, chemical reaction) has been recently developed by Zgrablich et al.20 (see also refs 21 and 22). In this section we briefly consider the effects of spatial energy correlation on adsorption isotherms by analyzing a simple model of adsorption in the presence of surface diffusion by a classical macroscopic mean-field approach.23 In the one-dimensional model considered, we assume an exponential distribution of adsorption energies and an exponential decay for the spatial correlation function. The spatial correlation properties of the energy landscape strongly influence the resulting adsorption equilibria as a consequence of the surface motion of the adparticles. The simulated conditions can be described in terms of adsorption, desorption, and diffusion contributions. In each site i (i ) 1, .., N ), characterized by the adsorption energy Ei and the corresponding coverage θi, the local adsorption and desorption rates are given by
rads,i ) ka(Ei)P(1 - θi)
(25)
rdes,i ) kd(Ei)θi where P is the gas-phase pressure and the adsorptiondesorption kinetic rates at a given temperature attain the form
ka(Ei) ) Ka ) constant
(26)
kd(Ei) ) Kd exp(-γEi) Surface diffusion effects are taken into account by introducing a hopping process between nearest neighboring sites driven by the energy field and by the thermic contribution to the two-dimensional motion within the adsorbed phase.24 (15) Jaroniec, M.; Madey, R. Physical Adsorption on Heterogeneous Surfaces; Elsevier: Amsterdam, 1988. (16) Steele, W. A. The Interaction of Gases with Solid Surfaces; Pergamon Press: New York, 1974. (17) Ripa, P.; Zgrablich, G. J. Phys. Chem. 1975, 79, 2118. (18) Rudzinski, W.; Nieszporek, K.; Cases, J. M.; Michot, L. I.; Villieras, F. Langmuir 1996, 12, 170. (19) Tovbin, Yu. K. Russ. J. Phys. Chem. 1990, 64, 461. (20) Zgrablich, G.; Mayagoitia, V.; Rojas, F.; Bulnes, F.; Gonzalez, A. P.; Nazzarro, M.; Pereyra, V.; Ramirez-Pastor, A. J.; Riccardo, J. L.; Sapag, K. Langmuir 1996, 12, 129. (21) Bulnes, F.; Riccardo, J. L.; Zgrablich, G.; Pereyra V. Surf. Sci. 1992, 260, 304. (22) Gonzales, A. P.; Pereyra, V. D.; Riccardo, J. L.; Zgrablich, G. J. Phys. Condens. Matter 1994, 6, 1. (23) Nicholson, D.; Parsonage, N. G. Computer Simulation and the Statistical Mechanics of Adsorption; Academic Press: London, 1982.
Simulation of Model Heterogeneous Surfaces
Langmuir, Vol. 12, No. 17, 1996 4277
We therefore confine ourselves to pointing out the crucial effect of the spatial distribution of adsorption energies and not only of their distribution function.
6. Generation of Correlated Dual Site-Bond Models
Figure 8. Coverage θt vs pressure P for different values of the parameter λ ) 5 × 10-3, 5 × 10-2, 5 × 10-1, characterizing the exponential decay of the energy correlation function.
Adparticles move randomly over the surface with hopping probabilities
hifi(1 ) exp[(Ei(1 - Ei)/kBT] 1 + exp[(Ei+1 - Ei)/kBT] + exp[(Ei-1 - Ei)/kBT]
(27)
The local hopping rates between nearest neighboring sites are therefore given by
rifj ) hohifjθi(1 - θj)
(28)
where j ) i ( 1. The factor ho is related to the characteristic adparticle mobility and can be tuned by assuming as a model parameter the ratio between the characteristic times for diffusion τd and for adsorption τa, i.e. φa ) τd/τa ) Ka/ho. The spatially discretized local equilibrium conditions
The expansion in exponentially correlated basis processes Y(x,λ) can be applied in the generation of correlated dual site-bond (DSB) models. The formulation of DSB models has been extensively discussed by Mayagoitia et al.5,6 and Riccardo et al.8 According to the analysis of Mayagoitia et al.5 and of Giona,25 site and bond energies can be considered as a dual representation of a single physical quantity, the adsorption energy field. In this section we focus on one-dimensional model surfaces. The analysis developed can be extended without major problems to higher dimensional models. Following Giona,25 we consider the generation of a DSB model by applying the case A method, in which the site energies are generated first and the bond energies are assigned subsequently by taking into account the local exclusion principle, EB e min{ES1,ES2}, where ES1 and ES2 are the energies of the sites connected by the bond. If we consider the energies Ei as a single process susceptible to dual representation in terms of site and bond energies, the correlation function C2E(x) can be defined in terms of the correlation functions C2ES, C2EB, and CESEB
1 C2E(2k) ) [C2ES(k) + C2EB(k)] 2 C2E(2k - 1) ) CESEB(k)
k ) 0, 1, ... (30) k ) 1, 2, ...
(31)
rads,i - rdes,i + ri+1fi + ri-1fi - rifi+1 - rifi-1 ) 0 (29) must be solved in order to obtain the equilibrium coverage N θi at a fixed value of the pressure P. θt ) ∑i)1 Figure 8 shows the adsorption isotherms obtained by changing the characteristic correlation parameter λ of the adsorption energies. The simulations were obtained for N ) 104 sites (one-dimensional lattice), β ) 0.1, γ ) 0.05, φa ) 10-4, and Emax/kBT =7.0 for three different values of λ in the low-pressure region. The increase in the correlation (lower values of λ ) generally implies higher coverages. In general, adsorption of monomers in the absence of surface diffusion is sensitive to surface topography only if molecules interact with each other. The presence of nonuniform energy fields induces correlation in particle motion on the surface. As a consequence, the correlation properties of the energy landscape become significant on the shape of the resulting isotherm. It should be observed that the results shown in Figure 8 depend on the choice of the functional form of the hopping rates, eq 27. It is clearly shown in Figure 8 that the spatial distribution of the adsorption energies is by no means irrelevant with respect to the nonlocalized equilibrium properties of an adsorbed phase. It would be interesting to dwell on the physical implications offered by such a phenomenon, but this lies outside the scope of this article. (24) Kapoor, A.; Yang, R. T.; Wong, C. Catal. Rev.sSci. Eng. 1989, 31, 129.
where the correlation functions C2A(x) (A ) ES, EB) and CESEB(x) are defined as
C2A(x) )
CESEB(x) )
〈(A(x+y) - 〈A〉)(A(y) - 〈A〉)〉 〈A2〉 - 〈A〉2
〈(ES(x+y) - 〈ES〉)(EB(y) - 〈EB〉)〉
x〈E2S〉 - 〈ES〉2x〈E2B〉 - 〈EB〉2
(32)
(33)
Let us express both site and bond energies as a function of a single Gaussian correlated stochastic process Y(x) possessing the correlation function C2Y(x). By following the same approach developed in section 2 (see eq 5 and the subsequent analysis), let G*S, G* B, be the nonlinear filters for bond and site energies,
G*S(Y) ) F*S-1(FY(Y))
G*B(Y) ) F*B-1(FY(Y)) (34)
-1 are the inverse functions of the where F*S-1 and F* B effective site and bond energy distributions F*S(E) and F*B(E). Let us consider the sites located at even values of the lattice coordinate x and the bonds located at odd values. The relation between C2E(x) and C2Y(x) is therefore given, for even x, by
(25) Giona, M. Unpublished.
4278 Langmuir, Vol. 12, No. 17, 1996
Adrover et al.
f(Y,ES1,ES2) ) B1 + (min{B2,ES1,ES2} - B1)KB(Y) (40) Let p(ES1,y,ES2) denote the probability density function for the three-dimensional random variable (ES1,y,ES2), where ES1 ) G* S(yi), y ) yi+1, ES2 ) G* S(yi+2), and yi ) Y(i). By definition, the effective bond energy distribution function F*B(EB) is given by
F*B(EB) )
∫∫∫f(y,E
S1ES2)eEB
Figure 9. C2E vs C2Y for uniform effective bond and site energy distributions. The dotted line is the curve C2E ) C2Y.
C2E(x) )
[
{∫
1 2
∞
dy1 -∞
∫
dy2 × -∞
[
∫BS dES1 ∫BS dES2 ∫-∞H (E ,B )p(ES1,y,ES2) dy + S E j (ES1,ES2) dES2 + 2∫E dES1∫S p ∫SE dES1 ∫SE pj (ES1,ES2) dES2 + B S H (E ,E ) p(ES1,y,ES2) dy 2∫E dES1 ∫E dES2 ∫-∞ 2
]
(G*S(y1) - 〈ES〉)(G*S(y2) - 〈ES〉)
∫-∞∞dy1 ∫-∞∞dy2
After some steps, we obtain
F*B(EB) )
∞
〈E2S〉 - 〈ES〉2
p(ES1,y,ES2) dES1 dES2 dy (41)
2
2
B
2
p(y1,y2) +
〈E2B〉 - 〈EB〉2
(35)
1
B
S1
1
1
2
p(y1,y2)
2
B
B
] }
(G*B(y1) - 〈EB〉)(G*B(y2) - 〈EB〉)
B
2
2
B
B
B
C2E(x) )
∫-∞dy1 ∫-∞dy2 ×
HB(EB,) ) K-1 B
∞
[x
]
(G*S(y1) - 〈ES〉)(G*B(y2) - 〈EB〉) 〈E2S〉 - 〈ES〉2x〈E2B〉 - 〈EB〉2
p(y1,y2) (36)
where p(y1,y2) is given by eq 10 and depends on C2Y(x). For example, in the case of the uniform effective pdf p*S(E), p*B(E), eqs 35 and 36 coincide and the curve C2E(x) vs C2Y(x) is independent of the value of Si and Bi (i ) 1, 2). Figure 9 shows the behavior of C2E(x) vs C2Y(x) in the case of uniform effective pdf’s. As can be seen, C2E is very close to C2Y. This situation occurs also for exponentially distributed site and bond energies (see eq 11 and Figure 1). Therefore, it is a fairly reasonable assumption that
C2E = C2Y
EB ) f(Y,ES1,ES2)
(38)
depending on the energies ES1, ES2 of the two sites connected by the bond. Let KB(Y) be a generic continuous monotonically increasing function mapping (-∞,∞) onto [0, 1]; then
{
f(Y,ES1,ES2) ) B1 + (B2 - B1)KB(Y) if ES1, ES2 g B2 B1 + (ES1 - B1)KB(Y) if ES1 e B2, ES2 e ES1 (39) B1 + (ES2 - B1)KB(Y) if ES2 e B2, ES1 g ES2 Equation 39 can be rewritten in compact form as
(
)
EB - B1 - B1
(43)
Note that in the two triple integrals on the right-hand side of eq 42 containing HB(EB,), the value of ( ) B2 in the first, ) ES1 in the second) is always greater than or equal to B2. By taking the derivative with respect to EB, it follows that
p*B(EB) ) dH (E ,B )
∫BS dES1 ∫BS dES2 p(ES1HB(EB,B2),ES2) BdEBB 2 B S 2∫E dES1 ∫E dES2 p(ES1,HB(EB,ES1),ES2) × 2
2
2
+
2
2
B
2
S1
dHB(EB,ES1) (44) dEB
(37)
Equation 37 simplifies the generation of correlated DSB energy landscapes. As discussed above, let us consider the generation of DSB one-dimensional energy landscapes by assuming the case A assignment. The local exclusion principle implies that the bond energies EB of a generic bond are assigned according to a nonlinear relation
(42)
where p j (ES1,ES2) ) ∫∞-∞p(ES1,y,ES2) dy and HB(EB,) is given by
and, for odd x, by ∞
S1
S1
The pdf p(ES1,y,ES2) can be written in terms of the pdf pY(yi,yi+1,yi+2) for the three-dimensional Gaussian random variable (yi,yi+1,yi+2) as
p(ES1,y,ES2) ) pY(G*S-1(ES1),y,G*S-1(ES2)) × dG*S-1(ES1) dG*S-1(ES2) | (45) dES1 dES2
|
where the pdf pY(yi,yi+1,yi+2) ) pY(y) is given by the usual expression for a multivariate Gaussian random variable
pY(y) )
1 1 exp - ytC-1y 2 (2π)3/2 det(C)
[
]
(46)
and the covariance matrix C for the three-dimensional normalized Gaussian random variable y attains the form
( )
1 c 1 c2 C ) c1 1 c 1 c2 c1 1
(47)
Simulation of Model Heterogeneous Surfaces
Langmuir, Vol. 12, No. 17, 1996 4279
Figure 10. p*B(E) vs E (B1 ) 1.4, B2 ) 3.4 au) for a uniform site energy distribution (S1 ) 2, S2 ) 4 au) and for the nonlinear filter HB given by eq 49. The dots are the numerical results and the lines the theoretical predictions, eq 44: (a) µ ) 1; (b) µ ) 4. λ ) 0.05, and the number of sites in the simulation is N ) 2 × 105.
Figure 11. C2E(x) vs x in the case µ ) 1, λ ) 0.05. The dots are the numerical results, and the line is the theoretical energy correlation function C2E ) exp(-λx).
where c1 ) C2Y(1) and c2 ) C2Y(2), since yi, yi+1 and yi+1, yi+2 are associated with nearest neighboring sites and yi, yi+2 with next nearest neighboring sites and the correlation function C2Y(x) is expressed in lattice units. The determinant of the covariance matrix is det(C) ) 1 - 2c21 c22 + 2c2c21. By performing a coordinate change, eq 45 can be rewritten as
p*B(EB) )
∫y(B∞ )dyi ∫y(B∞ )dyi+2 pY(yi,HB(EB,B2),yi+2) × 2
2
dHB(EB,B2) +2 dEB
∫y(Ey(B ))dyi × 2
B
dHB(EB,G*S(yi)) dyi+2 pY(yi,HB(EB,G*S(yi)),yi+2) (48) yi dEB
∫
∞
where y(EB) ) G*S-1(EB) and y(B2) ) G*S-1(B2). To consider an example, let us take the case of a nonlinear filter HB(EB,) given by
[(
HB(EB,) ) tan π
) ]
EB - B1 - B1
1/µ
-
π 2
(49)
where µ is a positive constant. Correspondingly,
KB(y) )
(
Figure 12. ψ(E) vs E for the overlapping degree I ) 0.75 (S1 ) 2, S2 ) 4, B1 ) 1.5, B2 ) 3.5 au) for different values of λ: (a) λ ) 0.1; (b) λ ) 0.5; (c) λ ) 1.37; (d) λ ) 2.0.
)
arctan(y) + π/2 π
µ
(50)
Let us further consider an exponentially correlated energy landscape C2E(x) ) exp(-λx), where λ is related to the correlation length Lc )1/λ (expressed in lattice units). Figure 10 shows the behavior of the effective bond energy pdf p*B(E) obtained in the case of a uniform site energy pdf (S1 ) 2, S2 ) 4 au) for λ ) 0.05 and for two different values of µ obtained from the numerical simulation of the DSB model compared with the theoretical prediction eq 44. The comparison between the correlation energy function obtained from the simulations and the prescribed one is shown in Figure 11. Equation 44 can of course be applied to solve the inverse problem, i.e. to find the function HB, which furnishes a given effective bond energy pdf with prescribed correlation properties. This problem is analogous to that examined by Giona25 and will be discussed in detail elsewhere. However, some general properties can be readily obtained from the analysis of eq 42. Solving the inverse problem means finding a suitable function HB, mapping [0, 1] onto (-∞,∞), which is the solution of eq 42 or eq 44. This can be achieved if the function
∫ES dES1 ∫SEdES2 pj (ES1,ES2) E E j (ES1,ES2) 2∫S dES1 ∫S dES2 p
ψ(E) ) F*B(E) - 2
2
1
S1
1
1
attains positive values. The function ψ(E) depends on the effective bond energy distribution, on the site energy distribution (in the case A assignment, local and effective site energy distributions coincide), and on the correlation properties between next nearest neighboring sites. In the uncorrelated case, p j (ES1,ES2) ) pS(ES1)pS(ES2) 2 and eq 51 reduces to ψ(E) ) F* B(E) - 2FS(E) + FS(E), 25 whose properties have been discussed by Giona. In the correlated case, let us consider the simple example of uniform effective site and bond energy distributions, in the case S1 ) 2, S2 ) 4 au and for exponentially correlated energy profiles C2E(x) ) exp(-λx). In the uncorrelated case, there is a threshold value of the overlapping degree Ic ) 0.75, such that above Ic it is not possible to generate a uniform effective bond energy distribution by applying the case A assignment procedure. The presence of correlations tends to remove this limitation. Figure 12 shows the behavior of ψ(E) in the interval [S1, B2] for Ic ) 0.75 and for different values of λ. As can be seen, with the decreasing of λ (i.e. increasing the correlation length L ) 1/λ) there is a critical value λc such that, for λ ) λc, ψ(E) is tangent to the E-axis. For values of λ smaller than λc, ψ(E) is strictly positive in [S1, B2] (by definition, ψ(E) is intrinsically non-negative in the interval [B1, S1] since FS(E) ) 0 in this region). Therefore, it is possible to define a function HB (and its inverse KB) solving eq 44 or equivalently eq 42, which, applied in the local assignment process of bond energies, gives an effective bond energy distribution F*B(E). The behavior of the critical correlation length Lc ) 1/λc vs the overlapping degree is shown in Figure 13.
4280 Langmuir, Vol. 12, No. 17, 1996
Figure 13. Lc ) 1/λc vs I in the case of a uniform effective site energy distribution (S1 ) 2, S2 ) 4 au) and of a uniform effective bond energy distribution F*B(E) ) (E - B1)/(B2 - B1). The line is a visual aid.
7. Concluding Remark In this article we present a general method to generate correlated energy fields in model surfaces. This method is the extension of a similar approach applied in the reconstruction of porous media. The concept of expanding a given stochastic process in terms of exponentially correlated Gaussian processes enables us to frame the determination of the expansion coefficients in the form of a linear problem similar to a Laplace transform inversion. The presence of correlation in the energy landscapes considerably modifies the behavior of experimentally measurable quantities (e.g. adsorption isotherms), especially in the case of mobile adsorption, in which surfacediffusion phenomena play a significant role. As the aim of this article is mainly methodological, the physical phenomenology in adsorption associated with
Adrover et al.
correlated energy patterns will be discussed in future papers. The analysis is also extended to the DSB model. Despite the formal difficulties intrinsic to the dual description of site and bond energies, the conceptual steps in the generation of DSB model surfaces are fairly simple and consist in reconstructing the correlation function of the Y(x) process, which can in most cases be considered approximately equal to C2E(x), and relating the effective bond pdf to the local pdf in the assignment of the bond energies. This relation, eq 44, depends on the correlation properties of the energy field. Moreover, we discuss in detail the effect of correlations in the case A assignment. In the absence of correlation, there is a threshold value of the overlapping degree above which it is not possible to generate a given uniform effective bond pdf. The presence of correlations in the energy field removes this limitation: above the critical threshold xc ) 0.75, there is a value of the correlation length Lc ) 1/λc above which the generation of the uniform bond pdf becomes feasible. It is important to observe that eq 44 can be used, by applying the same approach discussed by Giona,25 to solve the inverse problem, i.e. to find the nonlinear filter HB or its inverse KB, defined in eq 39, which enables us to generate a given effective bond energy pdf possessing the assigned correlation properties. The details of the solution of this inverse problem will be discussed elsewhere. Acknowledgment. One of the authors (A.A.) thanks J. Riccardo for useful discussions, and all the authors express their gratitude to W. Rudzinsky for organizing the ISSHAC Conference in Poland, during which most of this work was mapped out. LA960043M
