- Email: [email protected]
Application of the Replica Ornstein–Zernike Equations to Study Submonolayer Adsorption on Energetically Heterogeneous Surfaces
Journal of Colloid and Interface Science 219, 184 –189 (1999) Article ID jcis.1999.6459, available online at http://www.idealibrary.com on
Application of the Replica Ornstein–Zernike Equations to Study Submonolayer Adsorption on Energetically Heterogeneous Surfaces W. Rz˙ysko,* ,† O. Pizio,* S. Sokołowski,† and Z. Sokołowska‡ ,1 *Instituto de Quı´mica de la UNAM, Coyoaca´n 04510, Me´xico D.F., Mexico; †Department for the Modelling of Physico-Chemical Processes, Maria CurieSkłodowska University, Lublin 20031, Poland; and ‡Institute of Agrophysics, Polish Academy of Sciences, Dos´wiadczalna 5, 20290 Lublin, Poland Received April 13, 1999; accepted July 23, 1999
particles of remaining components can move in the field created by the quenched molecules. The theories of quenched– annealed systems have been applied to interpret the adsorption of gases in microporous systems. By the construction, the adsorbent obtained by the quenching process is highly geometrically and energetically nonuniform. Ideas similar to those used to explain gas adsorption in a heterogeneous microporous adsorbent can be also applied to describe adsorption on geometrically flat but energetically highly heterogeneous surfaces. Thus, the principal aim of this article is to model submonolayer adsorption on an energetically heterogeneous surface using replica methodology (6, 7). The energetic surface heterogeneity follows from the spatial distribution of adsorbing centers. In the present study the surface energetic centers are obtained by quenching some “particles” on a flat surface. In order to test the proposed approach we also perform Monte Carlo simulation.
We apply the replica methodology and the replica Ornstein– Zernike equations to calculate adsorption isotherms on an energetically heterogeneous surface. The surface is modeled by quenching a two-dimensional “fluid of centers.” After the quenching, we study gas adsorption using the replica Ornstein–Zernike equation with the hypernetted chain closure. The results of theoretical predictions are compared with the data obtained from the Grand Canonical Monte Carlo simulation, and the obtained agreement is satisfactory. © 1999 Academic Press Key Words: adsorption; energetic heterogeneity; replica Ornstein–Zernike equation; computer simulations.
I. INTRODUCTION
The importance of the role that energetic heterogeneity of solid surfaces plays in determining submonolayer adsorption isotherms has been stressed in several papers; see, e.g., (1–3). Early attempts to characterize gas adsorption on energetically heterogeneous surfaces were based on the application of the discrete models which mapped the energetic heterogeneity of an adsorbent into a finite set of individual surfaces, each characterized by a value of the adsorbate–adsorbent energy. By increasing the number of such homeostatic surfaces up to infinity, the next class of models, namely, the so-called patchwise models, have been proposed. According to those models, the total adsorption is treated as an average value of “local adsorptions” on individual, energetically homeostatic patches, averaged with a distribution of the number of patches vs their energy, i.e., with the so-called energy distribution function. Obviously, in the case of real heterogeneous adsorbents, the energetic centers are more or less randomly distributed, and, in order to be able to reproduce different energetically heterogeneous surfaces, different degrees of spatial correlation should be taken into account (1, 2). Recently, there has been much progress in theoretical studies of the so-called quenched–annealed systems (4 –11), i.e., the systems in which the particles of one component of a mixture are immobilized and form a matrix (or a network), while the 1
To whom correspondence should be addressed.
0021-9797/99 $30.00 Copyright © 1999 by Academic Press All rights of reproduction in any form reserved.
II. THEORY AND SIMULATION
To begin with, we describe first the model and next the computer simulation procedure. We investigate monolayer adsorption of a gas on a planar surface with energetic sites (centers) distributed on it. The distribution of the centers can be modeled similarly to a matrix formation in the case of the quenched–annealed systems; i.e., it can be obtained by quenching a “fluid of centers.” In order to do that, we consider a two-dimensional fluid consisting of particles of species 0 at the density r 0. At equilibrium, the microscopic state of that fluid is characterized by the pair correlation function h 00 (r). After quenching a two-dimensional structure is formed, and its microscopic state is still characterized by the distribution function h 00 (r) (6, 7). Now, we imagine that a gas of species 1 is adsorbed on an already prepared surface and forms a submonolayer. The gas has the chemical potential m 1 (or the activity a 1); its particles interact among themselves via the potential u 11 (r) and with the centers—via the potential u 10 (r). Depending on the density of the centers and on the shape of the potential u 10 (r), the resulting adsorbing potential field, V( x, y) ( x and y are the Cartesian coordinates on the surface), being the sum (12) of interactions
184
185
REPLICA ORNSTEIN–ZERNLIKE EQUATIONS
of a test gas particle with the quenched centers, may vary quite rapidly from point to point on the surface. In other words, our model mimics a situation of submonolayer adsorption on an energetically heterogeneous solid surface. Obviously, the shape of the field V( x, y) depends on the surface preparation. In our mental experiment it has been done by quenching a “fluid.” However, other methods, such as a two-dimensional diffusion limited aggregation, can be also applied. The usual theory describing systems with quenched components is based on the replica Ornstein–Zernike (ROZ) approach. According to this approach, the total pair correlation and the direct function of the species 0, i.e., h 00 (r) and c 00 (r), are obtained solving the Ornstein–Zernike (OZ) equation (4–11) h 00 2 c 00 5 r 0 c 00 # h 00 ,
b u ij~r! 1 ln@h ij~r! 1 1# 5 h ij~r! 2 c ij~r!, [1]
(the r-dependencies here and below are omitted for the sake of brevity and V denotes convolution), complemented by a closure relation. One can use the Percus–Yevick closure for that purpose; c 00 ~r! 5 $exp@2b u 00 ~r!# 2 1%$1 1 h 00 ~r! 2 c 00 ~r!%,
in the usual way (6, 7) as the sum of all graphs in the corresponding total correlation functions, h ij , with no nodal points. The correlation functions h b11 and c b11 , i.e., the socalled blocked parts of the total and direct gas– gas correlation functions, are the sums of all graphs in h 11 and c 11 , respectively, such that all paths between the two roots pass through at least one r 0 field point (4 –11). In many applications the equations given above have been solved using the Percus–Yevick approximation. However, in the case of systems with attractive forces a more sophisticated hypernetted chain (HNC) approximation is more appropriate (9). This approximation reads
[2]
where b 5 1/kT, T is the temperature, k is the Boltzmann constant, and u 00 (r) is the potential between a pair of particles of “0” type. The submonolayer adsorbate is considered at a two-dimensional adsorbate density r 1. Everywhere in what follows we choose the diameter of particle 1 as a length unit, i.e., s 1 5 1. The equilibrium between an adsorbed submonolayer phase and its bulk (three-dimensional) counterpart occurs at a constant chemical potential; in an experiment its value can be changed, e.g., by changing the bulk gas pressure. However, in the theoretical procedure we need to pick up the value for the two-dimensional density, r 1, first and calculate the chemical potential afterward (13). The ROZ equations are applied to describe the gas– center and gas– gas correlations in a quenched system of centers with the structure given by h 00 (r). The ROZ equations are the (4 –11) h 10 2 c 10 5 r 0 c 10 # h 00 1 r 1 c 11 # h 10 h 11 2 c 11 5 r 0 c 10 # h 10 1 r 1 c c11 # h 11 1 r 1 c b11 # h c11
for ~i, j! 5 ~1, 1! and ~1, 0!, c b11 ~r! 5 exp@h b11 ~r! 2 c b11 ~r!# 2 1 2 h b11 ~r! 1 c b11 ~r!. [4] The ROZ equation [3], in conjunction with Eq. [4], can be solved numerically by direct iterations. To determine the equilibrium between the adsorbed phase and its bulk (three-dimensional) counterpart we need to calculate the chemical potential, bm 1(r 1). The chemical potential of an adsorbate can be presented as the sum of an ideal and an excess term (13),
bm 1 ~ r 1, r 0 ! 5 bm 1,id ~ r 1! 1 bm 1,ex ~ r 1, r 0 !,
[5]
where the ideal gas contribution is bm 1,id ( r 1 ) 5 ln(r 1). The excess term is represented as a sum of two contributions,
bm 1,ex~ r 1 , r 0!
E E r1
5 bm 1,ex ~ r 1 5 0, r 0 ! 2
dr9
drc c11 ~r, r 9!,
[6]
0
where bm 1,ex ( r 1 5 0, r 0 ) is the chemical potential of species 1 at infinite dilution. The adsorption isotherm, i.e., the condition of equilibrium between the adsorbed monolayer and the bulk gas of density r 1,b , has the form
bm 1,b ~ r 1,b ! 5 ln~ r 1 ! 1 bm 1,ex ~ r 1 5 0, r 0 !
h c11 2 c c11 5 r 1 c c11 # h c11
E E r1
2
h 11 5 h c11 1 h b11
dr9
drc c11 ~r, r 9!.
[7]
0
c 11 5 c c11 1 c b11 .
[3]
Here, the direct correlation functions between gas particles and between gas and centers, c ij , (ij 5 11 or ij 5 10) are defined
In a typical case of adsorption of gases at low temperatures bm 1,b ( r 1,b ) can be accurately approximated by b p, where p is the bulk gas pressure.
RZ˙YSKO ET AL.
186
As we have already noted, the part bm 1,ex ( r 1 5 0, r 0 ) is the contribution to the chemical potential in the limit of zero adsorption, r 1 3 0. Formally, bm 1,ex ( r 1 5 0, r 0 ) 5 2ln ^exp[2 b V ( x, y)]& 5 ln K H, where K H is a constant that may be identified with the usual Henry constant (12), and ^. . .& denotes averaging over the entire surface and over all replicas. The average ^exp[2 b V( x, y)]& is connected with an insertion of a single particle of species 1 into the system containing quenched centers; i.e., it is equal to the average of the exponent of the mean force potential between a single particle of type 1 and the quenched centers,
^exp@2b V~ x, y!#& 5 r 0
E
drg 010 ~r!,
[8]
0 0 where g 10 (r) is the radial distribution function, g 10 (r) 5 h 10 (r; r 1 5 0) 1 1, evaluated in the limit of r 1 3 0. We have applied the theory quoted above to the simple system in which the interaction potentials between pairs of 00 and 11 particles are just hard-sphere type with diameters s 0 5 0.5 s 1 and s 1 5 1. The interactions between an adsorbate particle and a given center have been assumed in the form
H
e9 r , s0 u 10 ~r!/kT 5 2e exp@2r/a#/r r . s , 0
[9]
where e9 denotes the energy parameter. The potential [9] has Yukawa type tail, appropriate to model attraction between an adsorbate particle and an adsorbing center; however, an infinite repulsion branch of the Yukawa potential has been replaced by a finite barrier, allowing the adsorbate molecules “to enter” the center. In fact, monolayers built up on real substrates are at some distance from the geometrical surface plane of the underlying solid; i.e., the adsorbate particles when mapped on the surface plane “overlap” the adsorbing centers, and the aim of the applied model is to mimic this feature. Obviously, the proposed model is rather simple and the application of the theory to real adsorbing surfaces would require application of more sophisticated forms of u 10 (r). The system described above has been also investigated via grand canonical Monte Carlo simulations. For the purpose of our study, it is necessary first to generate a multiple set of surfaces (replicas) that differ by the distribution of particles of species 0. We assume that the configuration of an array of these particles (sites) corresponds to an equilibrium configuration of hard spheres of diameter s 0. This system has been generated and equilibrated by using canonical Monte Carlo simulation. Next, a grand canonical Monte Carlo simulation (GCMC) has been performed for species 1. The particles of species 1 are hard spheres of the unit diameter, and they interact with the centers via the potential [9]. The GCMC simulation experiment has been performed for a fixed configuration of the centers “0”
under the assumption that the gas–matrix interaction potential, u 10 (r), permits penetration of adsorbate species into the “interior” of the centers. Thus, we attempted to create (and annihilate) the particles in the entire simulation box. The GCMC simulations have been performed in a commonly used manner. One attempts to displace particles and to create and annihilate them in the entire box, keeping the overall acceptance ratio around 30%. In the majority of runs, equilibration has been reached after at most 10 6 Monte Carlo steps; the productive part of a run, during which the averages have been collected, involved approximately 10 7 steps. The input parameter for the GCMC run is the chemical potential of the adsorbate, m 1 /kT. All the quantities obtained in a simulation run for a given matrix configuration have been saved. We have performed five independent GCMC simulations (for different configurations of the centers) in a square box with the edge l 5 100 s 1 , and, as a final step, we have averaged the results. Actually, the adsorption isotherms for each of matrix configurations have not differed much. Even at high values of the chemical potential the deviation of adsorbed density from configuration to configuration has been of the order of 2–3%. III. RESULTS AND DISCUSSION
Our model calculations have been carried out assuming the parameters of the potential [9] equal to e 9/kT 5 2 e /kT 5 2 and a 5 4. The density of the centers was r 0 s 12 5 0.2 (i.e., r 0 s 02 5 0.05). The concentration of the active centers is thus similar in the case of some adsorbents, such as silica gel or alumina (14). We should note here that in our approach the attractive energy is negative, in contrast to the usual theories of adsorption on heterogeneous surfaces (1, 2), in which the attractive adsorbing energy is taken to be positive. The most energetic parts of the surface generated in the Monte Carlo simulations exhibited the energies of the order of 215kT, but there also were some parts of the surface with repulsive energy of the order of 2kT. Thus, the model surface considered by us is characterized by a reasonably broad spectrum of energies. The system considered by us is also characterized by strongly nonadditive interactions. However, the calculations performed (15, 16) for fluids in bulk random porous matrices have indicated that the ROZ theory leads to accurate description of the systems with nonadditive interactions. Each distribution of adsorbing centers generates the field V( x, y). At a given point x, y on the surface and for a given distribution of the centers (i.e., for a given replica), the potential V( x, y) is calculated by summing the interactions of a test adsorbate particle with all energetic centers, i.e., with all particles of the species 0. An example of the function V( x, y) is shown in Fig. 1. Note that only 1/64 of the entire surface used as the simulation box is displayed here. We see that the adsorbing potential field varies rapidly along the surface. Obviously, the change of the density r 0 and/or the parameters of the potential [9] may lead to a quite different field V( x, y). In
REPLICA ORNSTEIN–ZERNLIKE EQUATIONS
FIG. 1. An example of the function V( x, y) for a selected replica. Only 1/64 of the entire simulation box is shown here.
particular, an increase of the density r 0 would increase the average value of V( x, y). We also stress that the surface generated in our computer experiment surface does not correspond to the so-called “random surface” (1, 2), because the adsorbing centers are not randomly distributed. The correlations between a pair of the adsorbing centers gives the pair correlation function g 00 (r). In usual approaches to describe adsorption on energetically heterogeneous surfaces, workers usually have used the socalled energy distribution function, x(e), defined in such manner that x ( e )d e gives the fraction of the surface with the energy e. Obviously, the knowledge of V( x, y) allows us to calculate the function x(e). Assuming a small energy grid, De, we count the number of elements DxDy with a given value of the energy e. Figure 2 compares the function x(e), obtained for one replica, with that averaged over five replicas of the system. Note that this individual replica has been selected using the criterion of the largest deviations between its energy distribution function and the average function x(e). All energy distribution functions exhibit a Gaussian-like shape with a maximum at about e /kT ' 25. For the most energetic parts of the surface the energy e /kT is close to 215. The results presented in Fig. 2 demonstrate that all replicas are characterized by similar functions x(e); i.e., according to standard theories of adsorption on heterogeneous surfaces, they are energetically almost identical. Thus, although the plot of V( x, y), presented in Fig. 1, may not be representative for the entire surface, only small changes in x(e) evaluated for particular replicas should be observed, providing that the sampled surface is big enough. The dependence between the function x(e) and the parameters of the potential [9], as well as between the density r 0 and the pair correlation function h 00 , is quite complicated. We can say that an increase of r 0 would lead to a shift of the entire energy distribution function x(e) toward more attractive energies e, but not necessarily to wider energy distribution. In fact, one can expect that the surface would exhibit maximum heterogeneity at some intermediate density r 0. Details of the distribution function x(e) would depend on the pair correlation function h 00 . Figure 3 compares the adsorbate–adsorbate and adsorbate– center correlation functions, g ij (r) 5 h ij (r) 1 1, for (i, j) 5
187
(1, 1) and (1, 0) evaluated from the solution of the ROZ equation and from simulation. These functions illustrate the probabilities of finding of an adsorbate particle at the distance r from another adsorbate particle (the functions g 11 (r)) and from an adsorbing center (the functions g 10 (r)). Solid lines and diamonds are for adsorbate–adsorbate correlations, (i, j) 5 (1, 1), whereas dashed lines and crosses denote adsorbate– adsorbing center correlations, (i, j) 5 (1, 0). The latter functions have been plotted for r . s 0 only, because the assumed intracore repulsive barrier of the adsorbate–adsorbing center potential [9] is high enough and the correlation functions g 10 (r) are close to zero for r , s 0 , especially at low adsorbate densities. Figure 3, part a, shows the results at a low two-dimensional adsorbate density (i.e., at a low value of adsorption), r 1 5 0.205. The adsorbate–adsorbing center correlation function exhibits a single maximum, indicating that the adsorbed particles are located in the nearest vicinity of the adsorbing centers. Because of the attractive adsorbate– center potential, the first peak of the adsorbate– center correlation function, g 10 (r), is significantly higher that the first peak of the adsorbate–adsorbate correlation function, g 11 (r). The inset in Fig. 3a shows the center– center total correlation function, g 00 (r) 5 h 00 (r) 1 1. Obviously, because of low concentration of the centers, g 00 (r) is close to unity everywhere. The curves shown in Fig. 3, part b, have been evaluated at r 1 5 0.469. The contact value of g 10 (r) is now lower than that in Fig. 2a; however, the second peak in g 10 (r) is well developed. The presence of this peak indicates the formation of a secondary “coordination shell” around adsorbing centers. This
FIG. 2. The energy distribution function x(e), evaluated as described in the text. The solid line denotes the function averaged over replicas, whereas the dashed line is the function for one replica, exhibiting the largest deviations from the average.
RZ˙YSKO ET AL.
188
FIG. 3. Fluid–fluid (solid lines and diamonds) and fluid– center (dashed lines and crosses) distribution functions obtained from theory (lines) and from Monte Carlo simulation (points). The amount of adsorbed fluid is r 1 5 0.205 (part a), 0.469 (part b), and 0.720 (part c). The inset in part a shows the center– center distribution function.
secondary shell is located at a distance about 1.1s 1 from the first coordination shell around the adsorbing center, i.e., at a distance slightly larger than the diameter of the adsorbate particles. In comparison with the situation presented in Fig. 3a, the adsorbate–adsorbate correlations are also more pronounced. Further increase of the adsorption up to r 1 5 0.720, cf. Fig. 3c, leads to increase of the adsorbate– center and adsorbate–adsorbate correlations—the distribution functions plotted in Fig. 3c exhibit three very visible maxima. In all cases the agreement between simulated and theoretical results is good. Small deviations are observed only at the highest twodimensional densities, i.e., at the highest values of adsorption. Figure 4 compares the adsorption isotherm obtained from simulation and the theory. We have also plotted here the submonolayer isotherm on a model homogeneous surface, i.e., on a surface without a deposited matrix of centers. The latter isotherm has been evaluated solving the two-dimensional OZ equation hom hom hom h hom 11 2 c 11 5 r 1 c 11 # h 11 .
[10]
To make these calculations consistent with those for a heterogeneous surface, the two-dimensional OZ equation has been supplemented by the HNC closure for the (1, 1) correlation functions, Eq. [4]. In comparison with the isotherm for the considered heterogeneous surface, the isotherm for a homogeneous surface is shifted toward higher values of the chemical potential. Obviously, energetic heterogeneity, i.e., the existence of highly energetic centers on the surface, leads to strong increase of adsorption, especially at values of the chemical potential. Again, we stress a quite satisfactory agreement between the isotherms predicted theoretically and obtained from simulations. Note that the isotherm would have been identical also in the case in which the initial surface (without frozen centers) had extorted an additional (constant) potential, e 0. The imposed constant value of the adsorption energy e 0 would only change the value of the Henry constant. The model under consideration is very simple. However, it permits several interesting extensions. One of them may be thought of in terms of the sophistication of adsorbate–adsorbate interaction. If the attractive forces, like in the Lennard–
189
REPLICA ORNSTEIN–ZERNLIKE EQUATIONS
simulation studies. We hope to address some of these issues in our future work. ACKNOWLEDGMENTS Z.S. thanks KBN for financial support under grant 3 P06B 029 14. The work of S.S. has been supported by KBN under grant 3T09 A 8216.
REFERENCES
FIG. 4. A comparison of theoretical (line) and simulated (points) adsorption isotherms, evaluated as a function of the chemical potential (shifted by the logarithm of the Henry constant). Dotted line shows the adsorption isotherm on a homogeneous surface, evaluated solving the two-dimensional OZ equation [10] with the HNC closure.
Jones potential, were included, a study of phase transitions in sub- and monolayer adsorbed phase on energetically heterogeneous surfaces could be attempted. One can also consider more a complicated structure of matrices obtained, for example, by freezing a polydisperse mixture. In such case the theory proposed in (17) would be useful. There is also the possibility of going beyond monolayer systems by using an inhomogeneous ROZ equation (18, 19), or approaches based on the Born– Green–Yvon equation that uses a concept of the averaged density (20). Moreover, one can consider different mechanisms of formation of the energetically heterogeneous surface, with a possibility of arrangement of the adsorbing sites in fractal structures (21), obtained, for example, via diffusion limited aggregation process. Evidently, the adsorption of fluid mixtures on energetically heterogeneous surfaces would be of interest for the study of a competition between different species. Theoretical progress in this area, either in the framework of the ROZ equations or of more difficult inhomogeneous ROZ equations and related approaches, however, requires extensive
1. Jaroniec, M., and Madey, J., “Physical Adsorption on Heterogeneous Surfaces.” Elsevier, Amsterdam, 1988. 2. Rudzinski, W., and Everett, D., “Adsorption of Gases on Heterogeneous Surfaces.” Academic Press, New York, 1992. 3. Mamleev, M. S., Zolotareev, P. P., and Glagyshev, P. P., “Heterogeneity of Adsorbents. Phenomenological Models.” Nauka, Alma-Ata, 1989. [In Russian] 4. Madden, W. G., and Glandt, E. D., J. Stat. Phys. 51, 537 (1988). 5. Madden, W. G., J. Chem. Phys. 96, 5422 (1988). 6. Given, J. A., and Stell, G., Physica A 209, 495 (1994). 7. Given, J. A., and Stell, G., J. Chem. Phys. 97, 4573 (1992). 8. Lomba, E., Given, J. A., Stell, G., Weis, J. J., and Levesque, D., Phys. Rev. E 48, 233 (1993). 9. Meroni, A., Levesque, D., and Weis, J. J., J. Chem. Phys. 105, 1101 (1996). 10. Vega, C., Kaminsky, R. D., and Monson, P. A., J. Chem. Phys. 99, 3003 (1993). 11. Page, K. S., and Monson, P. A., Phys. Rev. E 54, R29 (1996). 12. Nicholson, D., and Parsonage, N. G., “Computer Simulation and the Statistical Thermodynamics of Adsorption.” Pergamon, London, 1982. 13. Ford, D. M., and Glandt, E. D., J. Chem. Phys. 100, 2391 (1994). 14. Oscik, J., “Adsorption.” PWS/Ellis Horwood, Chichester, 1982. 15. Ibarra-Bracamontes, L., Pizio, O., and Sokołowski, S., Rev. Mexicana Fisica 48, 1858 (1999). 16. Ibarra-Bracamontes, L., Pizio, O., Sokołowski, S., and Trokhymchuk, A., Molec. Phys. 96, 1341 (1999). 17. Ilnytskij, J., Pizio, O., Patrykiejew, A., and Sokołowski, S., J. Phys. Chem. B 103, 868 (1999). 18. Pizio, O., and Sokołowski, S., Phys. Rev. E 56, R63 (1997). 19. Kovalenko, A., Henderson, D., Pizio, O., and Sokołowski, S., Phys. Rev. E 57, 1824 (1998). 20. Perez, L., Sokołowski, S., and Pizio, O., J. Chem. Phys. 109, 1147 (1998). 21. Jensen, P., Barabasi, A-L., Larralde, H., Havlin, S., and Stanley, H. E., Fractals 4, 231 (1996).
Application of the Replica Ornstein–Zernike Equations to Study Submonolayer Adsorption on Energetically Heterogeneous Surfaces W. Rz˙ysko,* ,† O. Pizio,* S. Sokołowski,† and Z. Sokołowska‡ ,1 *Instituto de Quı´mica de la UNAM, Coyoaca´n 04510, Me´xico D.F., Mexico; †Department for the Modelling of Physico-Chemical Processes, Maria CurieSkłodowska University, Lublin 20031, Poland; and ‡Institute of Agrophysics, Polish Academy of Sciences, Dos´wiadczalna 5, 20290 Lublin, Poland Received April 13, 1999; accepted July 23, 1999
particles of remaining components can move in the field created by the quenched molecules. The theories of quenched– annealed systems have been applied to interpret the adsorption of gases in microporous systems. By the construction, the adsorbent obtained by the quenching process is highly geometrically and energetically nonuniform. Ideas similar to those used to explain gas adsorption in a heterogeneous microporous adsorbent can be also applied to describe adsorption on geometrically flat but energetically highly heterogeneous surfaces. Thus, the principal aim of this article is to model submonolayer adsorption on an energetically heterogeneous surface using replica methodology (6, 7). The energetic surface heterogeneity follows from the spatial distribution of adsorbing centers. In the present study the surface energetic centers are obtained by quenching some “particles” on a flat surface. In order to test the proposed approach we also perform Monte Carlo simulation.
We apply the replica methodology and the replica Ornstein– Zernike equations to calculate adsorption isotherms on an energetically heterogeneous surface. The surface is modeled by quenching a two-dimensional “fluid of centers.” After the quenching, we study gas adsorption using the replica Ornstein–Zernike equation with the hypernetted chain closure. The results of theoretical predictions are compared with the data obtained from the Grand Canonical Monte Carlo simulation, and the obtained agreement is satisfactory. © 1999 Academic Press Key Words: adsorption; energetic heterogeneity; replica Ornstein–Zernike equation; computer simulations.
I. INTRODUCTION
The importance of the role that energetic heterogeneity of solid surfaces plays in determining submonolayer adsorption isotherms has been stressed in several papers; see, e.g., (1–3). Early attempts to characterize gas adsorption on energetically heterogeneous surfaces were based on the application of the discrete models which mapped the energetic heterogeneity of an adsorbent into a finite set of individual surfaces, each characterized by a value of the adsorbate–adsorbent energy. By increasing the number of such homeostatic surfaces up to infinity, the next class of models, namely, the so-called patchwise models, have been proposed. According to those models, the total adsorption is treated as an average value of “local adsorptions” on individual, energetically homeostatic patches, averaged with a distribution of the number of patches vs their energy, i.e., with the so-called energy distribution function. Obviously, in the case of real heterogeneous adsorbents, the energetic centers are more or less randomly distributed, and, in order to be able to reproduce different energetically heterogeneous surfaces, different degrees of spatial correlation should be taken into account (1, 2). Recently, there has been much progress in theoretical studies of the so-called quenched–annealed systems (4 –11), i.e., the systems in which the particles of one component of a mixture are immobilized and form a matrix (or a network), while the 1
To whom correspondence should be addressed.
0021-9797/99 $30.00 Copyright © 1999 by Academic Press All rights of reproduction in any form reserved.
II. THEORY AND SIMULATION
To begin with, we describe first the model and next the computer simulation procedure. We investigate monolayer adsorption of a gas on a planar surface with energetic sites (centers) distributed on it. The distribution of the centers can be modeled similarly to a matrix formation in the case of the quenched–annealed systems; i.e., it can be obtained by quenching a “fluid of centers.” In order to do that, we consider a two-dimensional fluid consisting of particles of species 0 at the density r 0. At equilibrium, the microscopic state of that fluid is characterized by the pair correlation function h 00 (r). After quenching a two-dimensional structure is formed, and its microscopic state is still characterized by the distribution function h 00 (r) (6, 7). Now, we imagine that a gas of species 1 is adsorbed on an already prepared surface and forms a submonolayer. The gas has the chemical potential m 1 (or the activity a 1); its particles interact among themselves via the potential u 11 (r) and with the centers—via the potential u 10 (r). Depending on the density of the centers and on the shape of the potential u 10 (r), the resulting adsorbing potential field, V( x, y) ( x and y are the Cartesian coordinates on the surface), being the sum (12) of interactions
184
185
REPLICA ORNSTEIN–ZERNLIKE EQUATIONS
of a test gas particle with the quenched centers, may vary quite rapidly from point to point on the surface. In other words, our model mimics a situation of submonolayer adsorption on an energetically heterogeneous solid surface. Obviously, the shape of the field V( x, y) depends on the surface preparation. In our mental experiment it has been done by quenching a “fluid.” However, other methods, such as a two-dimensional diffusion limited aggregation, can be also applied. The usual theory describing systems with quenched components is based on the replica Ornstein–Zernike (ROZ) approach. According to this approach, the total pair correlation and the direct function of the species 0, i.e., h 00 (r) and c 00 (r), are obtained solving the Ornstein–Zernike (OZ) equation (4–11) h 00 2 c 00 5 r 0 c 00 # h 00 ,
b u ij~r! 1 ln@h ij~r! 1 1# 5 h ij~r! 2 c ij~r!, [1]
(the r-dependencies here and below are omitted for the sake of brevity and V denotes convolution), complemented by a closure relation. One can use the Percus–Yevick closure for that purpose; c 00 ~r! 5 $exp@2b u 00 ~r!# 2 1%$1 1 h 00 ~r! 2 c 00 ~r!%,
in the usual way (6, 7) as the sum of all graphs in the corresponding total correlation functions, h ij , with no nodal points. The correlation functions h b11 and c b11 , i.e., the socalled blocked parts of the total and direct gas– gas correlation functions, are the sums of all graphs in h 11 and c 11 , respectively, such that all paths between the two roots pass through at least one r 0 field point (4 –11). In many applications the equations given above have been solved using the Percus–Yevick approximation. However, in the case of systems with attractive forces a more sophisticated hypernetted chain (HNC) approximation is more appropriate (9). This approximation reads
[2]
where b 5 1/kT, T is the temperature, k is the Boltzmann constant, and u 00 (r) is the potential between a pair of particles of “0” type. The submonolayer adsorbate is considered at a two-dimensional adsorbate density r 1. Everywhere in what follows we choose the diameter of particle 1 as a length unit, i.e., s 1 5 1. The equilibrium between an adsorbed submonolayer phase and its bulk (three-dimensional) counterpart occurs at a constant chemical potential; in an experiment its value can be changed, e.g., by changing the bulk gas pressure. However, in the theoretical procedure we need to pick up the value for the two-dimensional density, r 1, first and calculate the chemical potential afterward (13). The ROZ equations are applied to describe the gas– center and gas– gas correlations in a quenched system of centers with the structure given by h 00 (r). The ROZ equations are the (4 –11) h 10 2 c 10 5 r 0 c 10 # h 00 1 r 1 c 11 # h 10 h 11 2 c 11 5 r 0 c 10 # h 10 1 r 1 c c11 # h 11 1 r 1 c b11 # h c11
for ~i, j! 5 ~1, 1! and ~1, 0!, c b11 ~r! 5 exp@h b11 ~r! 2 c b11 ~r!# 2 1 2 h b11 ~r! 1 c b11 ~r!. [4] The ROZ equation [3], in conjunction with Eq. [4], can be solved numerically by direct iterations. To determine the equilibrium between the adsorbed phase and its bulk (three-dimensional) counterpart we need to calculate the chemical potential, bm 1(r 1). The chemical potential of an adsorbate can be presented as the sum of an ideal and an excess term (13),
bm 1 ~ r 1, r 0 ! 5 bm 1,id ~ r 1! 1 bm 1,ex ~ r 1, r 0 !,
[5]
where the ideal gas contribution is bm 1,id ( r 1 ) 5 ln(r 1). The excess term is represented as a sum of two contributions,
bm 1,ex~ r 1 , r 0!
E E r1
5 bm 1,ex ~ r 1 5 0, r 0 ! 2
dr9
drc c11 ~r, r 9!,
[6]
0
where bm 1,ex ( r 1 5 0, r 0 ) is the chemical potential of species 1 at infinite dilution. The adsorption isotherm, i.e., the condition of equilibrium between the adsorbed monolayer and the bulk gas of density r 1,b , has the form
bm 1,b ~ r 1,b ! 5 ln~ r 1 ! 1 bm 1,ex ~ r 1 5 0, r 0 !
h c11 2 c c11 5 r 1 c c11 # h c11
E E r1
2
h 11 5 h c11 1 h b11
dr9
drc c11 ~r, r 9!.
[7]
0
c 11 5 c c11 1 c b11 .
[3]
Here, the direct correlation functions between gas particles and between gas and centers, c ij , (ij 5 11 or ij 5 10) are defined
In a typical case of adsorption of gases at low temperatures bm 1,b ( r 1,b ) can be accurately approximated by b p, where p is the bulk gas pressure.
RZ˙YSKO ET AL.
186
As we have already noted, the part bm 1,ex ( r 1 5 0, r 0 ) is the contribution to the chemical potential in the limit of zero adsorption, r 1 3 0. Formally, bm 1,ex ( r 1 5 0, r 0 ) 5 2ln ^exp[2 b V ( x, y)]& 5 ln K H, where K H is a constant that may be identified with the usual Henry constant (12), and ^. . .& denotes averaging over the entire surface and over all replicas. The average ^exp[2 b V( x, y)]& is connected with an insertion of a single particle of species 1 into the system containing quenched centers; i.e., it is equal to the average of the exponent of the mean force potential between a single particle of type 1 and the quenched centers,
^exp@2b V~ x, y!#& 5 r 0
E
drg 010 ~r!,
[8]
0 0 where g 10 (r) is the radial distribution function, g 10 (r) 5 h 10 (r; r 1 5 0) 1 1, evaluated in the limit of r 1 3 0. We have applied the theory quoted above to the simple system in which the interaction potentials between pairs of 00 and 11 particles are just hard-sphere type with diameters s 0 5 0.5 s 1 and s 1 5 1. The interactions between an adsorbate particle and a given center have been assumed in the form
H
e9 r , s0 u 10 ~r!/kT 5 2e exp@2r/a#/r r . s , 0
[9]
where e9 denotes the energy parameter. The potential [9] has Yukawa type tail, appropriate to model attraction between an adsorbate particle and an adsorbing center; however, an infinite repulsion branch of the Yukawa potential has been replaced by a finite barrier, allowing the adsorbate molecules “to enter” the center. In fact, monolayers built up on real substrates are at some distance from the geometrical surface plane of the underlying solid; i.e., the adsorbate particles when mapped on the surface plane “overlap” the adsorbing centers, and the aim of the applied model is to mimic this feature. Obviously, the proposed model is rather simple and the application of the theory to real adsorbing surfaces would require application of more sophisticated forms of u 10 (r). The system described above has been also investigated via grand canonical Monte Carlo simulations. For the purpose of our study, it is necessary first to generate a multiple set of surfaces (replicas) that differ by the distribution of particles of species 0. We assume that the configuration of an array of these particles (sites) corresponds to an equilibrium configuration of hard spheres of diameter s 0. This system has been generated and equilibrated by using canonical Monte Carlo simulation. Next, a grand canonical Monte Carlo simulation (GCMC) has been performed for species 1. The particles of species 1 are hard spheres of the unit diameter, and they interact with the centers via the potential [9]. The GCMC simulation experiment has been performed for a fixed configuration of the centers “0”
under the assumption that the gas–matrix interaction potential, u 10 (r), permits penetration of adsorbate species into the “interior” of the centers. Thus, we attempted to create (and annihilate) the particles in the entire simulation box. The GCMC simulations have been performed in a commonly used manner. One attempts to displace particles and to create and annihilate them in the entire box, keeping the overall acceptance ratio around 30%. In the majority of runs, equilibration has been reached after at most 10 6 Monte Carlo steps; the productive part of a run, during which the averages have been collected, involved approximately 10 7 steps. The input parameter for the GCMC run is the chemical potential of the adsorbate, m 1 /kT. All the quantities obtained in a simulation run for a given matrix configuration have been saved. We have performed five independent GCMC simulations (for different configurations of the centers) in a square box with the edge l 5 100 s 1 , and, as a final step, we have averaged the results. Actually, the adsorption isotherms for each of matrix configurations have not differed much. Even at high values of the chemical potential the deviation of adsorbed density from configuration to configuration has been of the order of 2–3%. III. RESULTS AND DISCUSSION
Our model calculations have been carried out assuming the parameters of the potential [9] equal to e 9/kT 5 2 e /kT 5 2 and a 5 4. The density of the centers was r 0 s 12 5 0.2 (i.e., r 0 s 02 5 0.05). The concentration of the active centers is thus similar in the case of some adsorbents, such as silica gel or alumina (14). We should note here that in our approach the attractive energy is negative, in contrast to the usual theories of adsorption on heterogeneous surfaces (1, 2), in which the attractive adsorbing energy is taken to be positive. The most energetic parts of the surface generated in the Monte Carlo simulations exhibited the energies of the order of 215kT, but there also were some parts of the surface with repulsive energy of the order of 2kT. Thus, the model surface considered by us is characterized by a reasonably broad spectrum of energies. The system considered by us is also characterized by strongly nonadditive interactions. However, the calculations performed (15, 16) for fluids in bulk random porous matrices have indicated that the ROZ theory leads to accurate description of the systems with nonadditive interactions. Each distribution of adsorbing centers generates the field V( x, y). At a given point x, y on the surface and for a given distribution of the centers (i.e., for a given replica), the potential V( x, y) is calculated by summing the interactions of a test adsorbate particle with all energetic centers, i.e., with all particles of the species 0. An example of the function V( x, y) is shown in Fig. 1. Note that only 1/64 of the entire surface used as the simulation box is displayed here. We see that the adsorbing potential field varies rapidly along the surface. Obviously, the change of the density r 0 and/or the parameters of the potential [9] may lead to a quite different field V( x, y). In
REPLICA ORNSTEIN–ZERNLIKE EQUATIONS
FIG. 1. An example of the function V( x, y) for a selected replica. Only 1/64 of the entire simulation box is shown here.
particular, an increase of the density r 0 would increase the average value of V( x, y). We also stress that the surface generated in our computer experiment surface does not correspond to the so-called “random surface” (1, 2), because the adsorbing centers are not randomly distributed. The correlations between a pair of the adsorbing centers gives the pair correlation function g 00 (r). In usual approaches to describe adsorption on energetically heterogeneous surfaces, workers usually have used the socalled energy distribution function, x(e), defined in such manner that x ( e )d e gives the fraction of the surface with the energy e. Obviously, the knowledge of V( x, y) allows us to calculate the function x(e). Assuming a small energy grid, De, we count the number of elements DxDy with a given value of the energy e. Figure 2 compares the function x(e), obtained for one replica, with that averaged over five replicas of the system. Note that this individual replica has been selected using the criterion of the largest deviations between its energy distribution function and the average function x(e). All energy distribution functions exhibit a Gaussian-like shape with a maximum at about e /kT ' 25. For the most energetic parts of the surface the energy e /kT is close to 215. The results presented in Fig. 2 demonstrate that all replicas are characterized by similar functions x(e); i.e., according to standard theories of adsorption on heterogeneous surfaces, they are energetically almost identical. Thus, although the plot of V( x, y), presented in Fig. 1, may not be representative for the entire surface, only small changes in x(e) evaluated for particular replicas should be observed, providing that the sampled surface is big enough. The dependence between the function x(e) and the parameters of the potential [9], as well as between the density r 0 and the pair correlation function h 00 , is quite complicated. We can say that an increase of r 0 would lead to a shift of the entire energy distribution function x(e) toward more attractive energies e, but not necessarily to wider energy distribution. In fact, one can expect that the surface would exhibit maximum heterogeneity at some intermediate density r 0. Details of the distribution function x(e) would depend on the pair correlation function h 00 . Figure 3 compares the adsorbate–adsorbate and adsorbate– center correlation functions, g ij (r) 5 h ij (r) 1 1, for (i, j) 5
187
(1, 1) and (1, 0) evaluated from the solution of the ROZ equation and from simulation. These functions illustrate the probabilities of finding of an adsorbate particle at the distance r from another adsorbate particle (the functions g 11 (r)) and from an adsorbing center (the functions g 10 (r)). Solid lines and diamonds are for adsorbate–adsorbate correlations, (i, j) 5 (1, 1), whereas dashed lines and crosses denote adsorbate– adsorbing center correlations, (i, j) 5 (1, 0). The latter functions have been plotted for r . s 0 only, because the assumed intracore repulsive barrier of the adsorbate–adsorbing center potential [9] is high enough and the correlation functions g 10 (r) are close to zero for r , s 0 , especially at low adsorbate densities. Figure 3, part a, shows the results at a low two-dimensional adsorbate density (i.e., at a low value of adsorption), r 1 5 0.205. The adsorbate–adsorbing center correlation function exhibits a single maximum, indicating that the adsorbed particles are located in the nearest vicinity of the adsorbing centers. Because of the attractive adsorbate– center potential, the first peak of the adsorbate– center correlation function, g 10 (r), is significantly higher that the first peak of the adsorbate–adsorbate correlation function, g 11 (r). The inset in Fig. 3a shows the center– center total correlation function, g 00 (r) 5 h 00 (r) 1 1. Obviously, because of low concentration of the centers, g 00 (r) is close to unity everywhere. The curves shown in Fig. 3, part b, have been evaluated at r 1 5 0.469. The contact value of g 10 (r) is now lower than that in Fig. 2a; however, the second peak in g 10 (r) is well developed. The presence of this peak indicates the formation of a secondary “coordination shell” around adsorbing centers. This
FIG. 2. The energy distribution function x(e), evaluated as described in the text. The solid line denotes the function averaged over replicas, whereas the dashed line is the function for one replica, exhibiting the largest deviations from the average.
RZ˙YSKO ET AL.
188
FIG. 3. Fluid–fluid (solid lines and diamonds) and fluid– center (dashed lines and crosses) distribution functions obtained from theory (lines) and from Monte Carlo simulation (points). The amount of adsorbed fluid is r 1 5 0.205 (part a), 0.469 (part b), and 0.720 (part c). The inset in part a shows the center– center distribution function.
secondary shell is located at a distance about 1.1s 1 from the first coordination shell around the adsorbing center, i.e., at a distance slightly larger than the diameter of the adsorbate particles. In comparison with the situation presented in Fig. 3a, the adsorbate–adsorbate correlations are also more pronounced. Further increase of the adsorption up to r 1 5 0.720, cf. Fig. 3c, leads to increase of the adsorbate– center and adsorbate–adsorbate correlations—the distribution functions plotted in Fig. 3c exhibit three very visible maxima. In all cases the agreement between simulated and theoretical results is good. Small deviations are observed only at the highest twodimensional densities, i.e., at the highest values of adsorption. Figure 4 compares the adsorption isotherm obtained from simulation and the theory. We have also plotted here the submonolayer isotherm on a model homogeneous surface, i.e., on a surface without a deposited matrix of centers. The latter isotherm has been evaluated solving the two-dimensional OZ equation hom hom hom h hom 11 2 c 11 5 r 1 c 11 # h 11 .
[10]
To make these calculations consistent with those for a heterogeneous surface, the two-dimensional OZ equation has been supplemented by the HNC closure for the (1, 1) correlation functions, Eq. [4]. In comparison with the isotherm for the considered heterogeneous surface, the isotherm for a homogeneous surface is shifted toward higher values of the chemical potential. Obviously, energetic heterogeneity, i.e., the existence of highly energetic centers on the surface, leads to strong increase of adsorption, especially at values of the chemical potential. Again, we stress a quite satisfactory agreement between the isotherms predicted theoretically and obtained from simulations. Note that the isotherm would have been identical also in the case in which the initial surface (without frozen centers) had extorted an additional (constant) potential, e 0. The imposed constant value of the adsorption energy e 0 would only change the value of the Henry constant. The model under consideration is very simple. However, it permits several interesting extensions. One of them may be thought of in terms of the sophistication of adsorbate–adsorbate interaction. If the attractive forces, like in the Lennard–
189
REPLICA ORNSTEIN–ZERNLIKE EQUATIONS
simulation studies. We hope to address some of these issues in our future work. ACKNOWLEDGMENTS Z.S. thanks KBN for financial support under grant 3 P06B 029 14. The work of S.S. has been supported by KBN under grant 3T09 A 8216.
REFERENCES
FIG. 4. A comparison of theoretical (line) and simulated (points) adsorption isotherms, evaluated as a function of the chemical potential (shifted by the logarithm of the Henry constant). Dotted line shows the adsorption isotherm on a homogeneous surface, evaluated solving the two-dimensional OZ equation [10] with the HNC closure.
Jones potential, were included, a study of phase transitions in sub- and monolayer adsorbed phase on energetically heterogeneous surfaces could be attempted. One can also consider more a complicated structure of matrices obtained, for example, by freezing a polydisperse mixture. In such case the theory proposed in (17) would be useful. There is also the possibility of going beyond monolayer systems by using an inhomogeneous ROZ equation (18, 19), or approaches based on the Born– Green–Yvon equation that uses a concept of the averaged density (20). Moreover, one can consider different mechanisms of formation of the energetically heterogeneous surface, with a possibility of arrangement of the adsorbing sites in fractal structures (21), obtained, for example, via diffusion limited aggregation process. Evidently, the adsorption of fluid mixtures on energetically heterogeneous surfaces would be of interest for the study of a competition between different species. Theoretical progress in this area, either in the framework of the ROZ equations or of more difficult inhomogeneous ROZ equations and related approaches, however, requires extensive
1. Jaroniec, M., and Madey, J., “Physical Adsorption on Heterogeneous Surfaces.” Elsevier, Amsterdam, 1988. 2. Rudzinski, W., and Everett, D., “Adsorption of Gases on Heterogeneous Surfaces.” Academic Press, New York, 1992. 3. Mamleev, M. S., Zolotareev, P. P., and Glagyshev, P. P., “Heterogeneity of Adsorbents. Phenomenological Models.” Nauka, Alma-Ata, 1989. [In Russian] 4. Madden, W. G., and Glandt, E. D., J. Stat. Phys. 51, 537 (1988). 5. Madden, W. G., J. Chem. Phys. 96, 5422 (1988). 6. Given, J. A., and Stell, G., Physica A 209, 495 (1994). 7. Given, J. A., and Stell, G., J. Chem. Phys. 97, 4573 (1992). 8. Lomba, E., Given, J. A., Stell, G., Weis, J. J., and Levesque, D., Phys. Rev. E 48, 233 (1993). 9. Meroni, A., Levesque, D., and Weis, J. J., J. Chem. Phys. 105, 1101 (1996). 10. Vega, C., Kaminsky, R. D., and Monson, P. A., J. Chem. Phys. 99, 3003 (1993). 11. Page, K. S., and Monson, P. A., Phys. Rev. E 54, R29 (1996). 12. Nicholson, D., and Parsonage, N. G., “Computer Simulation and the Statistical Thermodynamics of Adsorption.” Pergamon, London, 1982. 13. Ford, D. M., and Glandt, E. D., J. Chem. Phys. 100, 2391 (1994). 14. Oscik, J., “Adsorption.” PWS/Ellis Horwood, Chichester, 1982. 15. Ibarra-Bracamontes, L., Pizio, O., and Sokołowski, S., Rev. Mexicana Fisica 48, 1858 (1999). 16. Ibarra-Bracamontes, L., Pizio, O., Sokołowski, S., and Trokhymchuk, A., Molec. Phys. 96, 1341 (1999). 17. Ilnytskij, J., Pizio, O., Patrykiejew, A., and Sokołowski, S., J. Phys. Chem. B 103, 868 (1999). 18. Pizio, O., and Sokołowski, S., Phys. Rev. E 56, R63 (1997). 19. Kovalenko, A., Henderson, D., Pizio, O., and Sokołowski, S., Phys. Rev. E 57, 1824 (1998). 20. Perez, L., Sokołowski, S., and Pizio, O., J. Chem. Phys. 109, 1147 (1998). 21. Jensen, P., Barabasi, A-L., Larralde, H., Havlin, S., and Stanley, H. E., Fractals 4, 231 (1996).