Improved Method for Analysis of Energetic Heterogeneity of Surfaces

Langmuir 1996, 12, 441-449

441

Improved Method for Analysis of Energetic Heterogeneity of Surfaces from Adsorption Isotherms Vadim Sh. Mamleev* and Esen A. Bekturov A. B. Bekturov Institute of Chemical Sciences, Kazakh Academy of Sciences, Valikhanov str. 106, Kazakhstan, Almaty 480100, CIS Received July 8, 1994. In Final Form: September 14, 1995X A simple and fast algorithm (IRA: improved regularization algorithm) is proposed for evaluation of the distribution of adsorption sites (centers) as a function of their adsorption energies. The experimental adsorption isotherms are used in the algorithm as initial information. In the IRA, the regularization provides reproducibility of numerical results, but also nonnegativity of the distribution function is taken into account for physical correctness of the solution. After testing on a simulated distribution, the algorithm has been applied to the treatment of the adsorption isotherm of nitrogen on the surface of maximally hydroxylated silica. The Langmuir isotherm and two isotherms, which include corrections taking into account the laterial interactions for the patchwise and random heterogeneity models, were used for the description of local equilibrium in this system. For all models the calculated distributions show two pronounced peaks.

1. Introduction. It is well-known that the surfaces of all solid adsorbents are more or less heterogeneous. It is also obvious that heterogeneity has an influence on physicochemical properties of adsorption systems. Nevertheless, quantitative estimations of this influence are difficult even for the simplest systems. The complexity of solving the heterogeneity problem lies in the fact that only the overall properties of the adsorbents can be measured as a rule in the usual thermodynamic experiments, whereas for a characterization of heterogeneity, it is necessary to know the local properties which are changeable in the different microscopic regions of the surface. In the theory of localized monolayer adsorption on heterogeneous surfaces from the gas phase, the connection between overall (integral) and local properties is usually characterized by using the distribution function of adsorption site energies. In particular, the overall (experimental) adsorption isotherm is represented in the form

∫ θ(P,) F() d ) n(P)

na

b

a

(1)

where P is the partial pressure of adsorbate in the bulk phase; , a, and b are the adsorption energy and its minimal and maximal values; θ(P,) is the fraction of occupation of sites with adsorption energy at pressure P (local adsorption isotherm); F() d is the fraction of sites with energy between and + d on the surface; n(P) is the amount of adsorbate on the surface per unit area; and na is the adsorption capacity (maximal amount of adsorbate in the monolayer per unit surface area). If any restrictions for the function F() are absent, then F() may be found as a solution of the linear Fredholm equation of the first kind (eq 1). But in correspondence with the physical sense of the problem the F() must conform to the following constraints:

∫ F() d ) 1 b

a

F() g 0

(a e e b)

(2) (3)

If na can be measured in the experiment, constraint 2 * To whom correspondence should be addressed. Phone: (327) 61-18-38. Fax: 7-3272-61-57-65. E-mail: [email protected]. X Abstract published in Advance ACS Abstracts, November 15, 1995.

0743-7463/96/2412-0441$12.00/0

must be taken into account in the calculation of F() from eq 1. Otherwise the function naF() should be found instead of F(); further, na is calculated by integrating the function naF(). Constraint 2 is especially important if the measurements are carried out far from the plateau n(P)/na ≈ 1. If the maximum experimental pressure PM is not large enough for the monolayer to attain appreciable saturation (n(PM)/na , 1), it is not possible to determine na from eq 1. However, na may be estimated, e.g., by using the theory of Brunauer, Emmett, Teller (BET). Constraint 2, accordingly, should be provided for in calculations. In the last 20 years considerable progress1-4 in the analysis of surfaces heterogeneity has been achieved due to the appearance of numerical methods for the inversion of eq 1. A number of the most famous and successful algorithms have been named by their authors as HILDA,5 ALINDA,6 CAEDMON,7,8 CAESAR,9 and EDCAIS.10 Although while composing each new algorithm the authors usually have taken into consideration the shortcomings of the previous ones, an ultimate procedure for the solution of eq 1 has not yet been elaborated. So, new efforts are needed to improve the methods of heterogeneity analysis. We have suggested11 the use of the Tikhonov regularization method for the solution of eq 1, improved by using constraints 2 and 3. The IRA (improved regularization algorithm) method was employed, mainly11 in the ionexchange theory, although it can be applied with the same success in solving the problems of physical adsorption. Recently Szombathely et al.12 and Koopal and Vos13 described the INTEG12 program packet and the REMEDI13 (1) House, W. A. Colloid Sci. 1983, 4, 1. (2) Jaroniec, M.; Braüer, P. Surf. Sci. Rep. 1986, 6, 65. (3) Sircar, S.; Myers, A. L. Surf. Sci. 1988, 205, 353. (4) LamWan, J. A.; White, L. R. J. Chem. Soc., Faraday Trans. 1991, 87, 3051. (5) House, W. A.; Jaycock, M. Colloid Polym. Sci. 1978, 256, 52. (6) Koopal, L. K.; Vos, K. Colloids Surf. 1985, 14, 87. (7) Ross, S.; Morrison, I. D. Surf. Sci. 1975, 52, 103. (8) Sacher, R. S.; Morrison, I. D. J. Colloid Interface Sci. 1979, 70, 153. (9) Vos, C. H. W.; Koopal, L. K. J. Colloid Interface Sci. 1985, 105, 183. (10) Braüer, P.; Fassler, M.; Jaroniec, M. Thin Solid Films 1985, 123, 245. (11) Mamleev, V. Sh.; Zolotarev, P. P.; Gladyshev, P. P. Heterogeneity of sorbents: Phenomenological models. Nauka: Alma-Ata, 1989 (in Russian). (12) Szombathely, M. v.; Braüer, P.; Jaroniec, M. J. Comput. Chem. 1992, 13, 17. (13) Koopal, L. K.; Vos, C. H. W. Langmuir 1993, 9, 2593.

© 1996 American Chemical Society

442

Langmuir, Vol. 12, No. 2, 1996

Mamleev and Bekturov

method where regularization is used in combination with the NNLS (nonnegative least squares14) method that allows constraint 3 to be taken into account. This variant of regularization was known to us from the monograph of Tikhonov et al.,15 but we preferred to use the method of barrier functions16 for accounting for constraint 3. First, this method, being probably faster, enables one to attain approximately the same accuracy. Second, it can be applied in solving the problems of nonquadratic minimization. In particular, we used it for computation of discrete distributions.17 We shall try to substantiate the ideas of this numerical procedure in the second section of this paper. In the third section several new remarks concerning the alternative methods are made to facilitate the choice of appropriate numerical procedure for future investigations. In the fourth section is considered the special example in which the resolution of IRA can be tested by using the analytical solution. The results of employing IRA for treatment of an experimental isotherm of nitrogen on the surface of silica are shown in the fifth section. When the F() for this system was calculated, three models of local equilibrium were used:1-4

field theory (Bragg-Williams approximation) for a twodimensional lattice gas. Note that the parameter na must be known for the calculation of Θ(P) from the experimental function n(P), but often it is sufficient to have only a rough estimate for n a. For IRA, local isotherms from many other models could be used for the description of local equilibrium1-4 besides those enumerated above. These are not considered here because the purpose of the work is not the development of a thermodynamic theory but only the verification of IRA for the solution of eq 1. 2. Numerical Method for Evaluation of the Energetic Distribution. It is well-known that eq 1 belongs to the class of so-called ill-posed20 (or “incorrect” 15,21) problems. The classical solutions of such equations are very sensitive to the little errors which are inavoidable in physical experiments. Besides that, errors may arise when the integral is approximated by the finite sum, although errors of the latter kind may be decreased as much as one wants by diminishing the integration step. The criterion of solution accuracy for eq 1 is the functional

∆[F] ) (PM - P1)-1

Langmuir model θ(P,) ) P/[P + A0 exp(-/RT)]

(4)

(5)

model of random heterogeneity θ(P,) ) P/[P + A0 exp(-(ZwΘ+)/RT)]

M

1

∫

[na/n(P)]

model of patchwise heterogeneity (Fowler-Guggenheim isotherm18) θ(P,) ) P/[P + A0 exp(-(Zwθ + )/RT)]

∫PP {1 -

(6)

where R is the gas constant, T is temperature, Z is the coordination number of the surface lattice (at a close packing, e.g., Z ) 6), w is the energy of the pair interaction of nearest-neighbor admolecules (which may be estimated from the heat of condensation of the gaseous adsorbate into liquid), Θ ) n/na is the overall degree of occupation of the adsorption sites (the integral isotherm), and A0 is a constant which may be calculated from statistical physics.1 It is convenient to carry out the practical calculations in the scale of reduced energy η ) - RT ln A0. When the energy η is introduced, one can avoid manipulations of the constant A0. If this constant is known, then the function F(η) is transformed easily into F() because the scales η and are merely shifted from each other by RT ln A0. We note also that it is convenient to calculate the energies , η, and Zw in units of RT. Equations 5 and 6 approximately take into account the lateral interactions on the surface. They are very much alike, but their right-hand parts include different isotherms: the local isotherm in eq 5 and the integral one in eq 6. Both equations are derived19 on the basis of mean (14) Lawson, C. L.; Hanson, R. J. Solving Least Squares Problems; Prentice-Hall: Englewood Cliffs, NJ, 1974. (15) Tikhonov, A. N.; Goncharskii, A. V.; Stepanov, V. V.; Yagola, A. G. Regularization algorithms and a priori information; Nauka: Moscow, 1983 (in Russian). (16) Gill, P. E.; Murray, W.; Wright, M. H. Practical Optimization; Academic Press: London, 1981. (17) Mamleev, V. Sh.; Bekturov, E. A. Langmuir, in press. (18) Fowler, R. H.; Guggenheim, E. A. Statistical Thermodynamics; University Press: Cambridge, 1949. (19) Hill, T. L. J. Chem. Phys. 1949, 17, 762.

b

a

θ(P,) F() d}2 dP (7)

where P1 and PM are the limits restricting the experimental interval of pressure. Other functionals can be used depending on the features of the distribution of the experimental data on the scale of pressure. In particular, if the measurements extend over a very wide range of P, an integration on a logarithmic scale seems preferable, i.e.,

∫PP {1 [na/n(P)]∫ θ(P,) F() d}2P-1 dP

∆[F] ) [ln(PM/P1)]-1

M

1

b

a

(8)

If the number M of experimental measurements n(Pj) (1 e j e M) of the isotherm is not large enough to estimate accurately the mean integral error, then a summation instead of an integration should be used, M

∆[F] ) M-1

∑ {1 - [na/n(Pm)]∫ m)1

b a

θ(Pm,) F() d}2 (9)

where n(Pm) are the experimental values of n(P) known for M points. In the case of a representative sample of n(Pm) (1 e m e M), the computations with eqs 7, 8, and 9 should lead to quantities of the same order of magnitude. Any standard procedure of minimization allows one to choose the function F() corresponding to a zero minimum of ∆[F]. This function is the exact (classical) solution of eq 1; but owing to the errors in n(P) it will most probably fluctuate so sharply that it will lose any physical sense. Let us assume that n j (P) is the exact right-hand part of eq 1. Though the function n j (P) is unknown, we shall try to approach the classical solution Fh() corresponding to n j (P). Since the experimental function n(P) contains errors and differs from n j (P), the substitution Fh() in eq 7 (or in eq 8 or 9) will lead to a nonzero value of the functional (20) Merz, P. J. Comput. Phys. 1980, 38, 64. (21) Tikhonov, A. N.; Arsenin, V. Ya. The methods of solution of incorrect problems; Nauka: Moscow, 1979 (in Russian).

Analysis of Energetic Heterogeneity of Surfaces

Langmuir, Vol. 12, No. 2, 1996 443

∫PP [1 - nj (P)/n(P)]2 dP ≈ ∆h M

∆[Fh] ) (PM - P1)-1

1

where ∆ h is the mean square of a relative error in n(P). If peculiarities of the physical experiment are known, the magnitude of ∆ h usually may be estimated. If some function F() provides the inequality ∆[F] < ∆ h and conforms to constraints 2 and 3, it completely exhausts the initial physical information. However, an infinite multitude of such functions exists. Since these functions may differ considerably from each other, the problem arises of choosing the “best” solution which can be correctly physically interpreted. First of all the “best” solution should be reproducible. In other words, the experimental error in n(P) must not exert an essential influence on the calculated F(). Reproducibility of the solution, in its turn, may be attained only in those cases when supplementary properties are ascribed to F() a priori. If some function F() corresponds to these properties (and, consequently, is reproducible), it may be considered as the “best” solution. The requirement of a smooth solution is taken to be a priori information in the method of regularization15,21 which we are using here. This is realized by finding the minimum of the new functional ∆R[F] instead of ∆[F]. The functional ∆R[F], e.g., may be

∆R[F] ) ∆[F] + RΩ[F]

(10)

where R is a small nonnegative parameter of regularization and Ω[F] is Tikhonov’s stabilizer, e.g., having the form

Ω[F] )

∫ [F2 + q(dF/d)2] d b

a

where q is some weighting coefficient; if there are no special considerations to choose its value, one can take q ) 1. The function F() ) FR() giving the minimum of functional 10 depends on R. This parameter must be chosen so that, on one hand, a smoothed and reproducible solution would be obtained from eq 1 and, on the other hand, an insertion of this solution into eq 1 would not lead to an error which essentially exceeds the known experimental error in n(P). For instance, it is reasonable to find h . This the parameter R providing the equality ∆[FR] ) ∆ variant of the regularization procedure is a typical problem of conditional minimization, viz., the finding of FR() requires

Ω[F] w min

(11)

∆[F] ) ∆ h

(12)

subject to

In this case one can see that R-1 is a Lagrange multiplier that allows one to seek an absolute minimum of the functional,

R-1∆R[F] ) Ω[F] + R-1∆[F] w min

The reglarized solution FR() is stable to small errors in n(P). Consequently, in correspondence with physical sense, if the function Fh() is smooth enough, then F() ) FR() = Fh(). Condition 2 may be taken into account by using the following functional:

∆R,R1[F] ) ∆[F] + RΩ[F] + R1(1 -

b

a

where R1 is the nonnegative parameter of the “penalty”. When R1 f ∞, condition 2 holds exactly. However, a more correct choice of R1 is that value when condition 2 has a relative error which is comparable with the experimental error in the measurement of na. Finally, condition 3 may be taken into account by the method of barrier functions.16 It is realized when one more item is added to functional 10,

∆R,R1,Rs[F] ) ∆R,R1[F] + Rs

∫ φ()[ln F() + D]2 d b

a

(14)

where Rs is a small positive parameter, D is a positive number (D = 5-15), and φ() is a positive weighting function, by assuming that the searched for solution is the formally known function, it is convenient to choose φ() to be F() so that φ() can be excluded from the final equations. As F f 0 we have ln F f -∞; hence the minimum ∆R,R1,Rs[F] is reached for F() > 0. When Rs f ∞, the solution providing the minimum of functional 14 obviously is F() ) exp(-D). In principle, for Rs . 0 the minimization of functional 14 may be a stable procedure without the stabilizer Ω[F], i.e., for R ) 0. In other words, at finite values of Rs the second term in eq 14 possesses a regularizing effect. However, we shall consider the other case, when Rs f 0. Negative logarithms with large absolute values correspond to the negligibly small values of F(). Therefore, as Rs f 0, the second term in the right-hand part of 14 is mainly determined by the result of the integration in intervals where F() , exp(-D). At the same time the magnitude of ∆R,R1[F] is not changed appreciably if we actually have F() , exp(-D) in intervals of integration where F() should be 0. In other words, the barrier term in eq 14 perturbs the solution only in close proximity to the boundary F() ) 0 and, consequently, allows one to find the conditional minimum of the functional ∆R,R1[F] subject to constraint 3 with a very high accuracy. The numerical minimization of ∆R,R1,Rs[F] is implemented after replacement of the functional by its discrete analog, i.e., after calculation of the integrals by using the method of finite sums, M

∆R,R1,Rs[F] ≈ ∆R,R1,Rsd[{Fj}] )

(13)

N

R1(1 instead of the conditional minimum 11. For the optimal R to be determined, minimization 13 should be carried out repeatedly. Then the parameter R may be found as a numerical solution of the algebraic equation 12, i.e., h. ∆(R) ) ∆[FR] ) ∆ Practical calculations, however, show that reproducibility of FR() is reached when R lies within a very wide range. Thus, there is no necessity to secure equality ∆ ) ∆ h ; to attain the desirable accuracy it is usually enough to orient oneself on the approximate correspondence ∆ = ∆ h.

∫ F() d)2 dP

∑

N

(1 -

m)1

N

cjFjUm,jh)2lm + ∑ j)1

N-1

cjFjh)2 + R[∑cjFj2h + q ∑ (Fj+1 - Fj)2/h] + ∑ j)1 j)1 j)1 N

Rs

cjφ(j)(ln Fj + D)2h ∑ j)1

where Um,j ) θ(Pm,j)na/n(Pm); Fj are the desired values of the function F() at N equidistant points ({j} (1 e j e N)); h ) j+1 - j and cj are the step and weighting coefficients for the integration (e.g., for the trapezoidal rule: c1 ) 1/2, cN ) 1/2, cj ) 1 (2 e j e N - 1)); and lm is the step of

444



integration. Thus, for eq 7,

l1 ) (P2 - P1)/(2d) lM ) (PM - PM-1)/(2d) lm ) (Pm+1 - Pm-1)/(2d) 2emeM-1

d ) PM - P1

and for eq 8,

l1 ) ln(P2/P1)/(2d) lM ) ln(PM/PM-1)/(2d) lm ) ln(Pm+1/Pm-1)/(2d) 2emeM-1

d ) ln(PM/P1)

and for eq 9,

lm ) 1/M

(1 e m e M)

For reducing the problem in question to an algebraic system it is necessar to find the partial derivatives of ∆R,R1,Rsd[{Fj}] for every varialbe Fk, to prescribe φ(j) ) Fj, and to substitute Fj ) exp(yj). After these operations we can obtain the following system of N nonlinear equations for yk: N

M

exp(yj)cj(R1 + ∑ Um,jUm,klm) + R{exp(yk) ∑ j)1 m)1

h

(q/ck)Λ[exp(yk)]} + Rsyk ) R1 - RsD + M

∑ lmUm,k m)1

(1 e k e N) (15)

where Λ is the finite difference analog of the second derivative of F():

{

(Fk+1 - 2Fk + Fk-1)/h2 Λ[Fk] ) (F2 - F1)/h2 (FN-1 - FN)/h2

(2 e k < N) (k ) 1) (k ) N)

The substitution Fj ) exp(yj) is made for two reasons. First, it allows one to avoid undetermined expressions when Fj e 0 and, hence, it is impossible to calculate the logarithm; and second, exponentials are linearized better in comparison to logarithms. The last circumstance makes it possible to use the fast iterative method of Newton for solving system 15. The method is realized by representing the exponents using the formula (c+1) exp(y(c+1) ) ) exp(y(c) - y(c) j j )(1 + yj j )

where c is the counter of iteration. If initial values of F() are given in accordance with the condensation approximation,1,2

F() ) (RTna)-1(dn/d ln P)|P)A0exp(-/RT) then almost complete convergence of the iteration process is reached after three to five cycles. There are many parameters in the problem under consideration, but their specific values only weakly affect the final result. They can be chosen, for example, as R1 ) 103, Rs ) 10-6, D ) 10, q ) 1. Only the choice of R

requires a numerical experiment. This parameter was varied according to R ) 10-m (m ) 1, 2, 3, ...). The minimal R which causes the chaotic oscillations of solution to be eliminated was taken as an optimal R. If the ∆ h is known and small enough and all parameters (besides R) are set, then, in principle, repeatedly solving the problem at different R will provide the equality ∆[FR] )∆ h . Essential complications, however, arise when ∆ h is unknown or when the minimum of ∆ exceeds ∆ h on account of constraints 2 and 3. Although we did not encounter such situations, the consideration of both of them is deemed important. It was shown by Merz20 that, for a choice of R with an unknown ∆ h , the method of generalized cross-validation (GCV) can be successfully used. The idea of the method is based on the assertion that the regularized solution FR() may serve for a prediction of the exact dependence of n j (P); the farther an experimental point n(Pm) is situated j (Pm) may be predicted. from n(Pm), the worse the value of n The optimal R is found as result of global minimization of the function M

∑ m)1

Hg(R) )

[n(Pm) - na

N

Fjm(R)cjhθ(Pm,j)]2lm ∑ j)1

where Fjm(R) (1 e j e N) denotes the regularized solution obtained when the mth data point is deleted from the experimental set. The minization of Hg(R) is a timeconsuming method because for each trial value of R the computation of the solution {Fjm(R)} should be repeated M times. Merz suggested20 an ingenious technique that allows one to reduce the computational work. However, the method cannot be employed with constraint 3. The method of finding a quasioptimal R15,21 does not require the value of ∆ h as well, and besides this, it is faster than the minimization of Hg(R). The quasioptimal R is a value for which the solution FR() manifests the least sensitivity on R. The measure of such sensitivity may be expressed as the following function: N

Hq(R) )

∫ [dFR()/d ln R]2 d ≈ R2∑cjh[dFj(R)/dR]2 b

a

j)1

(16)

Having computed the solution Fj(R) (1 e j e N) at different R one can find the derivatives in eq 16 by numerical differentiation. Constraints 2 and 3, as well as any other a priori information, should increase the stability of the inversion of eq 1. However, each constraint increases the minimum of ∆[F]. If the exact minimum of ∆[F] under constraints 2 and 3 is equal to ∆* and ∆* > 0, eq 1 has no classical solution Fc() corresponding to the equality ∆[Fc] ) 0. As to choice of R, it can be made from the system 11 and 12 only for the case in which ∆[Fc] ) 0 < ∆* < ∆[FR] ) ∆(R) ηb/RT. The EDCAIS10 algorithm is based on the formula derived by Temkin and Levich29 for the Langmuir model:

Fd(E)|E)ξ ) RTF()|P)A0exp(-/RT) ≈ -[dΘ/dξ - (π2/3!)d3Θ/dξ3 + (π4/5!)d5Θ/dξ5] The algorithm computes the derivatives of the function Θ(-ln P) up to and including the fifth order. However, even an accurate third derivative of the experimental function Θ(-ln P) is doubtful. Nevertheless, EDCAIS yields good results for artificial distributions simulated for testing the algorithm. At the same time, if a local isotherm is chosen unsuccessfully, the result of EDCAIS actually cannot be interpreted.3 The method of LamWan and White4 uses a Fourier transform of convolution 17. The method was tested for simulated distributions and showed good results. Unfortunately, too many points were taken for calculating the integral isotherm: M ) 1000. If perturbations were introduced for each point by the randomizing generator, the illustrations shown4 allow one to study only the influence of high-frequency noise on the stability of the inversion of eq 17, whereas the available experimental data usually are characterized by low-frequency noise. The special form of the local isotherm restricts considerably the sphere of application of this method, but it appears to suffer also from other shortcomings. First, the Θ(ξ) and 1/[1 + exp(ξ)] functions have no classical Fourier transforms. In order to overcome this difficulty the authors recommend calculating the Fourier integral for (29) Temkin, M.; Levich, V. Zh. Fiz. Khim. 1946, 20, 1441.

the difference

τ(ξ) ) n(ξ)/na - 1/[1 + exp(ξ)]

(18)

However, for this operation to be used in practice one needs to know the na value exactly, since otherwise the calculated function Fd(E) may lose physical sense. For example, it may occur that the calculated function Fd(E) (or F()) does not conform to eq 2 in spite of participation of na in the calculations according to eq 18. In the method described in the previous section this possibility is excluded owing to the self-consistency of eq 2 in the minimization, viz., the calculated Fd(E) (or F()) does conform to eq 2. Second, for the Fourier convolution method to be numerically performed, the limits for the pressure and the energy should be set in strong correspondence with each other, viz., a/RT - ln A0 ) Ea ) -ln PM, b/RT - ln A0 ) Eb ) -ln P1. In other words, the interval of integration [RT ln(A0/PM), RT ln(A0/P1)] is actually used in eq 1 instead of the interval [a,b]. Unfortunately, most experiments are carried out within too narrow a range of pressure to provide acceptable accuracy for integrating eq 1 within this interval. Authors have advised the use of an extrapolation of the integral isotherm to extend artificially the range of pressure, but from our viewpoint, the extrapolated isotherm cannot include additional physical information in comparison to the original isotherm. Moreover, the extrapolation may distort the form of the calculated distribution and may lead to spurious structure. The question of the choice of integration limits is discussed in the fifth section here. Experimental isotherms are as a rule smooth functions. Therefore 20-50 points measured in a logarithmic scale of pressure completely reflect the physical features of an adsorption system. A subsequent increase of the number of points on the isotherm does not lead to additional information. On the other hand, 40-200 nodes on a scale of energy are sufficient to approximate the majority of experimental isotherms by using the method of regularization. For the present we do not see any reasons to use other methods for the heterogeneity analysis. Results and Discussion 4. Testing the Algorithm. The testing of the algorithm usually includes two stages. First a “synthetic” overall isotherm, which corresponds to a known distribution, is analyzed by computing the integral in eq 1. The evaluation of the integral may be called a direct problem. Thereupon the inversion of the equation is carried out by the algorithm under test. This procedure with regard to the first one is a reverse problem. For verification of the reliability of IRA it is reasonable to dwell upon the special example in which integral 1 is calculated analytically. The particular case with the local Langmuir isotherm seems to be the most simple and suitable for such calculations. In this case equation 1 may be written in the form

Θ(P) )

∫0∞[P/(P + K)]β(K) dK

(19)

where K ) A0 exp(-/RT) ) exp(-η/RT) is the local adsorption coefficient and β(K) is the distribution density of adsorption sites with respect to K (the energy distribution is connected with β(K) by the relationship F() ) Kβ(K)/(RT)). We assumed that in eq 19 the energy may vary between its “mathematical” limits: from -∞ to ∞.


Langmuir, Vol. 12, No. 2, 1996 447

Figure 2. Isotherm calculated from 23 (solid line) and the superposition of two Langmuir isotherms (dashed line). Pressure is expressed as a dimensionless variable. Figure 1. Distribution corresponding to eq 22. The points are the result of solving the reverse problem.

By using properties of a Laplace transform30 we are able to write the solution of integral eq 19 in the form

2

Θ(P) )

a +i∞ +i∞ exp(Ks){(2πi)-1∫a -i∞ ∫aa-i∞ 1

β(K) ) (2πi)-1

1

(πP/Kjσj)1/2 exp[(P + Kj)2/PKjσj] × ∑ j)1 erfc[(P + Kj)/(PKjσj)1/2] (23)

2

2

exp(Ps)[Θ(P)/P] dP} ds ) ∞ a +i∞ +i∞ {∫0 {(2πi)-1 ∫a -i∞ exp(Ps)(K + ∫aa-i∞ a +i∞ ∞ P)-1 dP}β(K) dK} exp(Ks) ds ) (2πi)-1 ∫a -i∞ {∫0 1

(2πi)-1

After substituting eq 22 in eq 21 and using the tables of Laplace transforms,30 we have

where

2

1

2

1

∫0xexp(-x2) dx

erfc(x) ) 1 - 2/π1/2

1

exp(-Ks) β(K) dK} exp(Ks) ds (20) where i is the imaginary unity, a1 and a2 are arbitrary positive constants, and s is an intermediate parameter. It follows from eqs 20 that double integration may be used for finding the isotherm from the known function β(K) (or F():

∫0∞{∫0∞exp[-s(P + K)] ds}β(K) dK ) ∞ ∞ P∫0 exp(-sP){∫0 exp(-sK) β(K) dK} ds

Θ(P) ) P

(21)

The convenience of formulas 20 and 21 is that they allow one to solve eq 1 (or 19) by using tables of Laplace transforms.30 In particular, analytical calculations are possible for the distribution 2

∑ j)1

RTF() ) RT

2

Fj() )

Kβj(K) ) ∑ j)1

2

gj(K/πKjσj)1/2 exp[-(K - Kj)2/σjKKj] ∑ j)1

(22)

where gj, Kj, and σj are some positive parameters. On the scale of ln(K) the dependence of eq 22 is close to the superposition of two Gaussian distributions (see Figure 1). One can readily see that for distribution 22 the necessary normalization 2 is valid:

∫-∞∞Fj() d ) gj

∫0∞βj(K) dK ) gj

Besides this, when P f 0, distribution 22 conforms to Henry’s law: 2

∫0 K-1βj(K) dK ) P(g1/K1 + g2/K2) ∑ j)1

Θ(P) ) P

∞

(30) Korn, G. A.; Korn, T. M. Mathematical Handbook; McGrawHill: New York, 1968.

In a number of works,1,2,4,9 test isotherms have been generated by the numerical calculation of integral 1. An analytical expression for the isotherm allows one to test the algorithm more objectively due to the lack of integration errors at the stage of solving the direct problem. The following parameters have been assigned for calculation: g1 ) 0.2, g2 ) 0.8, σ1 ) σ2 ) 1, K1 ) 200, K2 ) 0.5, and -3.5 e ln K e 8.5. With these parameters the overall isotherm is close to the superposition of two Langmuir isotherms (see Figure 2) corresponding to the case σ1 ) σ2 ) 0. Before solving the reverse problem, isotherm 23 was multiplied by (1 + (5 × 10-3)Y), where the Y are quasirandom numbers31 distributed almost evenly on the interval [-1, 1]. Calculations have been carried out with R ) 10-4, M ) 41, N ) 41, P1 ) 5 × 10-1, and P2 ) 104. The points on the integral isotherm have been distributed equidistantly in the scale of ln P. The results of the reproduction of RTF() are represented in Figure 1. It was taken that R1 ) 0; nevertheless, the normalization 2 for the approximate solution is fulfilled with high accuracy (0.998 instead of 1). 5. Investigation of a Nitrogen Adsorption on Maximally Hydroxylated Silica. In Figure 3 the adsorption isotherm of N2 on maximally hydroxylated silica25 at 78 K is represented using the scale of ln P. To solve the reverse problem local isotherms 4-6 were used. The magnitude Zw/RT ) 3.8 has been chosen25 to take account of the lateral interaction effects for isotherms 5 and 6. Since the value na was considered as an unknown, eq 1 was solved without constraint 2 by setting R1 equal to 0. First the function naF was calculated; further it was integrated, i.e., na was computed; and finally naF was divided by na. Thus the solution turned into a normalized function. The reverse problem was solved on the scale of reduced energy η ) - RT ln A0 ) -RT ln K. For plotting F on the scale of η it is not necessary to know A0. The limits of energy for all models were set equal: -7 e η/RT e 8. (31) Forsithe, G. F.; Malcolm, M. A.; Moler, C. B. Computer Methods for Mathematical Computations; Prentice-Hall: Englewood Cliffs, NJ, 1977.

448


Figure 3. Calculated isotherm (solid line) of nitrogen adsorption on maximally hydroxylated silica at 78 K and the data of House25 (points, M ) 35), P in Torr. The dashed line shows the adsorption capacity calculated by using the Langmuir isotherm.

Figure 4. Results of solving the reverse problem (N ) 61) by IRA: (1) distribution obtained for the Langmuir local isotherm (eq 4); (2) result for a patchwise heterogeneity model with the local isotherm of Fowler-Guggenheim (eq 5); (3) result for the model of random heterogeneity with the local isotherm (eq 6). The distributions are plotted on a scale of the reduced energy η ) -RT ln K.

The results of a F(η) calculation are represented in Figure 4. As seen for all models of local equilibrium, the distributions have two well-resolved peaks. For employing the regularization method there is no necessity that the limits of pressure and energy be brought into accord with each other. In particular, the width of the energy interval may be chosen with some reserve to attain the desirable accuracy in calculating integral 1 for any pressure P within the interval P1 e P e PM. However, the physical interpretation of the solution should be made with the understanding that the function F(η) loses accuracy in the regions η < -RT ln PM and η > -RT ln P1. If the range P1 e P e PM is so wide that the overall isotherm is sensitive to the form of the distribution within the whole region ηa e η e ηb, then the function F(η) must decay at the ends of the energy interval. Solving the problem with regard to naF(η) and integrating the function, one can find a realistic value of na. It is seen in Figure 4 that distributions obtained with isotherms 5 and 6 do not vanish at the points ηa and ηb; consequently, with these distributions the parameter na cannot be found. Only for the Langmuir model can a reasonable value of na be evaluated. The parameter na ) 15.5 µmol/m2 (10.7 Å2 per molecule) results from the normalization of F(η). This value was employed for calculation of Θ(P) in the model of random heterogeneity. The na value calculated by the HILDA25,32 algorithm is a little smaller: 14.4 µmol/m2 (11.5 Å2/molecule). House adduced the mean relative deviations between theoretical (32) House, W. A. J. Chem. Soc., Faraday Trans. 1 1978, 74, 1045.


and experimental overall isotherms. The calculations with HILDA lead to the deviation 0.34% for the local isotherm of Langmuir. In our calculations with this local isotherm, being estimated according to formula 9, the relative error is 0.28%. Hence, we carried out the minimization of the functional ∆[F] approximately up to the value that House obtained by using HILDA. To evaluate the energetic distributions for models 5 and 6, House employed the technique of regularization. The mean deviations in this case were 1.07% and 0.14% for isotherms 5 and 6, respectively. We obtained by using IRA the deviations 0.25% and 0.12%. Our results for local isotherms 5 and 6 are close to those obtained by house25 because the corresponding distributions within the interval ηa e η e ηb are positive functions and, consequently, constraint 3 may be ignored. However, some differences between our results and those of House25 are observed. Our curves appear to be more smooth. Probably, stronger smoothing is reached due to the larger amount of nodes on the energy scale (in our calculations N ) 61 while House used N ) 20). Any additional disagreements may be explained by different choices of the integration interval. This question needs to be discussed. If the pressure region is so narrow that, at minimal pressure P1, the adsorption sites with energy -RT ln P1 , η2 < η e ηb are saturated, while, at maximal pressure PM, for sites with energy ηa e η < η1 , -RT ln PM Henry’s law is observed, then, without any loss of accuracy, we can rewrite eq 1 in the form

∫ηη F(η)[P exp(η/RT)] dη + na∫ηη θ(P,η) F(η) dη + η na∫η F(η) dη ) n(P) (24)

na

1

2

a

1

b

2

In the region ηa e η < η1 we have P exp(η/RT) , 1; therefore, the first term in eq 24 is negligibly small and may be omitted. In other words, we have no information about F(η) in the region ηa e η < η1 and we cannot estimate, consequently, the full adsorption capacity na. It is preferable to solve the problem with regard to naF(η); this function may be reproducible for the energy within the interval η1 e η e η2. It can be normalized, just as in our calculations for convenience of comparison of the distributions (see Figure 4), but the evaluated na will depend on the choice of the limit ηa. With decreasing ηa, the na value will increase. Let us assume that within the region -∞ < η e η1 the function F(η) has a uniform distribution F(η) ) Fu/RT, where Fu is dimensionless constant Fu , 1. If the maximal pressure PM is so small that PM exp(η1/RT) , 1, then the integral

∫-∞η F(η) exp(η/RT) dη ) naFuP exp(η1/RT) 1

naP

may be less than the experimental error in determining the overall isotherm, whereas the physical incorrectness of the function F(η) is obvious because

(η1 - ηa)naFu/RT f ∞ ∫ηη F(η) dη ) ηlim f-∞

lim na

ηaf-∞

1

a

a

This simple consideration allows one to conclude that in the absence of any information about ηa we cannot find the parameter na as an a posteriori result and, consequently, the extrapolation of the integral isotherm toward high pressure is impossible. Again, when solving eq 1 at fixed na without accounting for constraint 2, we cannot find the function F() that satisfies the normalization. Thus, if PM exp(ηa/RT) , 1, then it is preferable to set the


Langmuir, Vol. 12, No. 2, 1996 449

lower limit of the integration ηa (or a) on the basis of the estimations of statistical physics. Only in the case in which ηa is set correctly and the function F() vanishes near the point ηa can one find a realistic value of na. Another way out for correct normalization of F() is the use of constraint 2 when solving the problem. When -RT ln P1 , η2 < ηb, we have no information about behavior of the function F(η) in the region η2 e η e ηb. One can assert, however, that in the region [η2,∞) the function F(η) has the finite integral

∫η∞F(η) dη ) na[1 - F(η2)]

na

2

(25)

Any form of the function F(η) within the region η2 e η < ∞ for which integral 25 has fixed magnitude will conform with the experimental isotherm. If the upper limit ηb is set arbitrarily as some large constant, then in the region η2 < η e ηb the regularization algorithm prefers to choose the uniform distribution (see curve 2 in Figure 4) for which dF/dη ) 0. We note that for local isotherms 4 and 6 in comparison with the case of the Fowler-Guggenheim local isotherm 5, the overall isotherm in the range of low pressures carries more information about the part of the distribution with high adsorption energy, since the saturation of sites on patchwise surfaces occurs at lower pressures due to the lateral interactions between admolecules (compare curves 1 and 3 with curve 2 in Figure 4), the distribution for the Fowler-Guggenheim local isotherm relative to that calculated for local Langmuir isotherm being shifted by about Zw toward lower energy of adsorption.1-3 Because of the constancy of the value of integral 25 with an increase in the upper limit of integration, the distribution density will decrease at the ηb end:

lim F(η2) ≈ lim F(ηb) ) lim [1 - F(η2)]/(ηb - η2) ) 0

ηbf∞

ηbf∞

ηbf∞

but for the F(ηb) to vanish, an unrealistically large limit ηb needs to be set. With the overall isotherm having been measured within a narrow range of pressure, the upper limit of integration ηb, as well as the lower limit ηa, should be chosen on the basis of additional physical estimations. All distributions calculated by IRA are smoothed; therefore, at fixed ηa and ηb they are reproducible in spite of error in n(P). It is most interesting that for each of the three models of local equilibrium the distributions obtained are seen as two peaks. The maxima on the distribution curves must be attributed to definite adsorption centers of two different kinds. To interpret this fact on a molecular level one should know the structure of the atomic environment surrounding a N2 molecule which is localized at a point of minimum potential energy for its interaction with the surface. The existence on a silica surface of two kinds of adsorption centers was explained32,33 by two kinds of surface hydroxy groups: “free” and “hydrogen bonded”. Indeed, the existence of such vicinal hydroxyls was confirmed by an infrared spectroscopic study of silica surfaces.34 It was supposed32,34 that, owing to the high (33) Waksmundzki, A.; Toth, J.; Jaroniec, M.; Rudzinski, W. Rocz. Chem. 1975, 49, 1003. (34) Davydov, V. Ya.; Kiselev, A. V.; Zhuravlev, L. T. Trans. Faraday Soc. 1964, 60, 2254.

quadrupole moment, N2 moleculs exhibit specific interactions with hydroxyls. Besides this, the hypothesis was proposed32 that interaction of N2 with “free” groups is much stronger than with “bonded” ones. Adsorption of N2 on silica with different concentrations of surface hydroxyls was thoroughly investigated by Aristov and Kiselev.35-37 After degassing at 400 °C silica loses almost all “bonded” hydroxyl groups; their content is about 50% of the overall amount of hydroxy groups on the surface. However, the adsorption isotherm of N2 is not changed appreciably after such treatment of silica. After degassing of silica at 600 °C, when “free” hydroxyls are also leaving the surface, adsorption of N2 is essentially decreased. Since in the presence or absence of “bonded” hydroxyls the overall isotherm and, consequently, the energetic distribution are not changed, the peak at low energy on the distribution curve may be attributed to sites of the surface either containing “bonded” hydroxyls or containing no hydroxyls at all, whereas the peak at high energy, apparently, corresponds to sites containing “free” hydroxyls which are able to form weak hydrogen bonds with N2.38 It is very difficult to estimate the influence of hydroxyl groups on adsorption of N2 quantitatively. Unfortunately, the specific interactions of N2 with the surface hydroxyls are manifested against a background of its nonspecific (dispersion) interactions with the surface. Additional investigations are needed to ascertain whether a correlation between the energies of specific and nonspecific interaction exists or not. 6. Conclusions. We have considered employment of IRA for treatment of the isotherm of adsorption of nitrogen on a maximally hydroxylated silica surface. For all of three models used for localized adsorption, the calculated energetic distributions show two pronounced peaks. Apparently, the high-energy peak may be explained by the existence of “free” hydroxy groups on the surface. The algorithm described here yields reproducible results; this is the main advantage of IRA in comparison with the many other known algorithms. Reproducibility is achieved because IRA selects the most smooth solution from a great number of those which are consistent with the experimental data. Resolution of IRA depends on experimental error. If a solution is a sharply oscillating function (contains many peaks), then one should be careful in interpretation of the results. An increase of the regularization parameter leads to strong smoothness and to disappearance of some peaks in the distribution. If, in spite of this, theoretical and experimental isotherms differ within experimental error, then the peaks must be considered as computational artifacts. Only in the case in which the possibility of smoothness is restricted by a small value of experimental error can one proceed to the physical interpretation of an energetic distribution. So, with the development of the technique of adsorption experiments and with more accurate definitions of the models for local equilibrium, it will be possible to investigate smaller and smaller details of the energetic heterogeneity of the surface. LA940534B (35) Aristov, B. (36) Aristov, B. (37) Aristov, B. (38) McDonald,

G.; Kiselev, A. V. Zh. Fiz. Khim. 1963, 37, 2520. G.; Kiselev, A. V. Colloid J. 1965, 27, 246. G.; Kiselev, A. V. Zh. Fiz. Khim. 1964, 38, 1984. R. S. J. Am. Chem. Soc. 1957, 79, 850.

Improved Method for Analysis of Energetic Heterogeneity of Surfaces

Recommend Documents