Chapter 8. Computer simulation of surface diffusion in adsorbed phases

W. Rudzifiski, W.A. Steele and G. Zgrablich (Eds.) Equilibria and Dynamics of Gas Adsorption on Heterogeneous Solid Surfaces Studies in Surface Science and Catalysis, Vol. 104 9 1997 Elsevier Science B.V. All rights reserved.

451

C o m p u t e r S i m u l a t i o n of Surface Diffusion in A d s o r b e d P h a s e s William Steele Department of Chemistry, The Pennsylvania State University, 152 Davey Laboratory, University Park, PA 16802, U.S.A.

1

Introduction

Computer simulation is a now well-known technique for the evaluation of the thermodynamic properties of gases adsorbed on solid surfaces or in pores. The basic idea is to use the power of the computer to bypass the simplifying assumptions made in the past to pass from a molecular model of the adsorption system to the evaluation of its observable, macroscopic properties. In essence, one defines the solid adsorbent in terms of its interaction energies with gas molecules that come in contact with its surface. This model is completed by defining the interactions of the adsorbate molecules with each other. The simulation which can then be performed relies on one of two general algorithms: Monte Carlo or molecular dynamics [1, 2]. For both, tricks are needed to make such systems mimic the macroscopic adsorption of a gas on (or in) a solid: the potential energies of interaction must be modeled accurately, but not so accurately as to prevent the computer from evaluating them in a reasonable time; and periodic boundary and minimum image conditions must be applied to allow one to associate the properties of a few hundred molecules with those of a macroscopic fluid. These well-known ideas have been utilized in a very large number of problems explored using this technique. Features that are perhaps unique to adsorption are that the solid adsorbent is taken to be rigid so that one has a mechanical system of adsorbate particles with translational and rotational degrees of freedom moving in the external force field of the solid. Vibrational degrees of freedom for the adsorbed molecules are usually neglected, at least when the molecule is small- obviously, intramolecular motions are an important aspect of the simulations of long-chain adsorbed molecules. The molecular dynamics algorithm produces molecular trajectories for interacting molecules over periods of time ranging up to a nanosecond. This is accomplished by numerically integrating the classical equations of motion for each molecule of adsorbate to determine the time-dependent trajectories followed. Time-averages of quantities such as the average potential energy of the entire sample of adsorbate molecules can then, by the postulates of statistical mechanics, be equated to ensemble averages. Of more importance to the present discussion is the fact that the information about the time-dependent positions of the molecules in the system allow one to evaluate dynamical properties such as mean square displacements. It is well-known that these displacements are related to self-diffusion constants (or tracer diffusion constants, as they are sometimes called) in a straightforward but exact manner. An alternative calculation of the diffusion constant relies on the fact that molecular velocities are also known at each instant in time. This allows one to evaluate velocity autocorr~elation functions and, by integration of this function over time, to evaluate the self-diffusion coefficient in a calculation which is completely equivalent to the mettiod based on the long time limit of the slope of a plot of mean square displacement versus time. The formalism and the explicit equations for the self-diffusion constant

452

are written out in the chapter in this book by Zgrablich. In fact, other dynamical properties of a fluid such as the shear viscosity and/or the thermal conductivity can also be evaluated using suitable analogues to the mean square displacements or the velocity autocorrelation functions [2]. An important point is that all these quantities can be evaluated for fluids in geometries where the confining dimensions are of the order of nanometers. Furthermore, for inhomogeneous fluids, the dependence of the transport coefficients upon position within the fluid and upon the direction of the flux (of mass, momentum or energy, for diffusion, shear viscosity and thermal conduction, respectively) relative to the confining geometry can also be evaluated. It should also be emphasized here that these very detailed results can be generated by either molecular dynamics or by a suitably modified form of the Monte Carlo algorithm. Here, we will discuss some of the results that have been obtained over the past decade or so for diffusional motion in fluids on surfaces. In the application of the Monte Carlo algorithm to transport, dynamics are not followed per se, but the passage of the system through phase space is evaluated by making a lengthy series of small random displacements of the molecules. These trim displacements are accepted or rejected by a probability rule (hence, Monte Carlo) that has been shown rigorously to generate the statistical ensemble of interest. Many different types of ensemble can be produced (Grand Canonical, isobaric-isothermal, canonical, etc.) by appropriate modification of the probability rules that govern the random process. The tricks of periodic boundary conditions, rigid solid, etc. mentioned above for molecular dynamics are also utilized in Monte Carlo simulations. Since there is no time variable in this algorithm, one might think that evaluation of transport properties via Monte Carlo simulation is not possible. However, many workers have argued that this is unnecessarily restrictive. The "kinetic Monte Carlo" algorithm is based on a master equation with phenomenological rate constants [3]. In this way, molecules trace out trajectories in space and time that allow one to evaluate mean square displacements (in particular) and thus, a kind of diffusion coefficient. The problem is that the units of time remain undefined; nevertheless, many useful results have been obtained in this way [4]. In fact the kinetic Monte Carlo algorithm is indispensable for adsorbed phases that have been modeled as lattice gases. In this case, motion is due to instantaneous jumps between lattice sites and molecular dynamics is not applicable. This is a particularly appropriate model for the diffusive motion of strongly adsorbed (i.e., chemisorbed) species, especially on clean metallic surfaces - see for example, [5]. However, the lattice gas model is widely used for systems other than chemisorption - specifically, nearly all current models of heterogeneous surfaces involve the concept of adsorption on distributions of sites with varying adsorption energy. In this case, jump models plus kinetic Monte Carlo are employed to evaluate surface diffusion - the only added assumption needed is the specification of the energy barriers to jumps between neighboring sites. Simulation studies of such model systems are rather wide-spread [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17]. Over the past decade or two, molecular dynamics simulations of transport in adsorbed fluids have been reported for a large number of model systems. We give here a brief and not exhaustive review of this work and subsequently discuss in detail simulations of self- (or tracer-) diffusion in thin films (one or two atomic layers at most) of gases adsorbed on model surfaces. Many of the molecular dynamics simulations of transport in confined fluids have been concerned with gases sorbed in pores. Initially, model pores with a simple geometry were studied. Two of the most popular are slit pores with planar parallel walls and straight-walled cylinders.

453

For such models, one must then specify the roughness of the pore walls. Perfectly smooth walls can be studied, but in more realistic systems, one attempts to include the atomic structure of the wall in a way that is simple enough to allow one to evaluate fluid-wall interactions reasonably efficiently. Within the category of transport of pore fluids, one should distinguish two rather different cases: 9 Pores that are essentially filled with adsorbed fluid [18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32] 9 Pores that have submonolayer or at most, one or two strongly held layers of fluid on the pore walls (see below). In the first case, one is interested in the behavior of either diffusive or shear flow in the small-size limit where the nature of the boundary conditions has a significant effect on the hydrodynamics of flow, i.e., on the solutions to Fick's equation for diffusion or to the Navier-Stokes equation for viscous flow down a microscopic tube (Poiseuille flow) or between two planar walls (Couette flow). In the second case, the molecular motion in a thin layer on the wall tends to be rather similar to that for a fluid adsorbed on a free surface - that is, the presence of another wall at some distance from a very thin adsorbed film has a minimal influence on the dynamics of the molecules in the film except in the case of extreme constriction where the second wall is at a separation distance of a couple of atomic diameters. Well known examples of transport through constrictions are fluids in the pore systems of zeolites, [33, 34, 35, 36, 37, 38, 39] where the sizes of the openings to the cages inside the solid are of the order of a molecular diameter, thus producing the molecular sieve effect. As an example of the simulation of molecular motions of rare gas atoms, Figure 1 shows the trajectories followed by xenon atoms at 300 K in the pore system of a zeolite known as zeolite rho [33]. This crystal contains two interlaced non-connecting cubic lattices of cages with diameters _~ 12 ~t. The cages are connected by short tubes with diameters of "~ 3.6 )1. When xenon is brought into contact with this system, the atoms first are adsorbed quite strongly in the connecting tubes and then begin to fill the cages. Figure 1 shows three cages as the volumes within the hexagonal rings that actually define one end of the connecting tubes to the six neighbors of each cage. The small dots show the limited trajectories followed by atoms trapped in the tubes and vibrating on their sites at this point. All these sites are filled, but the periodic b o u n d a r y conditions cause some of them to be invisible (each atom of this kind is half on one face and half on the face on the opposite side, but is shown here as a single atom). If we focus our attention on the pair of cages vertically aligned in Fig. 1, one of which is a periodic image, we see t h a t the 6 atoms that are in each cage are translating across the inner wall of their cage in a random walk that is produced by collisions with neighboring xenon atoms. During the duration of this simulation, one atom leaves the tube that connects the central cage and its image above. It moves into the upper image cage, where it collides with the image atoms (not shown) and diffuses within the second cage. As it makes the change of cages, an atom from the cage below moves into the vacated tube. In this way, one sees that inter-cage diffusion controls the overall diffusive motion of the xenon sorbed in this pore system. Since only one such event has occurred in 37 pic0seconds, it is evident that a very long simulation will be needed to observe a sufficient number of inter-cage jumps to give a statistically satisfactory estimate of the jump rate in this system. A problem that frequently arises in this and other zeolite studies is that the rigid lattice model assumed in this work may give results for transport that are significantly

454

different from those expected if the atoms making up the solid lattice are allowed to undergo vibrational displacements. This problem is magnified in importance in porous systems that are "tight" - i.e., where the sorbate and the pore dimensions are sufficiently close to be at the edge of the molecular sieving that occurs when the gas atoms are too large to pass through the cage

F i g u r e 1 Computer generated trajectories for 16 xenon atoms in cages taken from the two independent networks of connected cages in zeolite rho. Also shown is a periodic image of one of the cages. Six of the atoms are tightly held in the openings to the cages (not all are shown) and the remaining 10 are evenly divided between the two cages. During the time of this simulation, one atom jumps from its position in the window into a neighboring image cage and is replaced in the window by a different atom from the cage. openings. If we now focus on the molecular dynamics simulations of self-diffusion of molecules confined to few-layer films on flat surfaces, one can again divide the problem into homogeneous and heterogeneous model adsorbents. In the homogeneous case, we begin with the fact that a very large effort has been devoted to experimental, theoretical and simulation studies of gases adsorbed on graphite. This work has included a number of simulation studies of diffusion on this surface together with a few experimental reports [40]. As is well-known, the surface of this solid can readily prepared so that it is nearly entirely made up of large-area perfect basal planes. Furthermore, the basal plane is believed to be remarkably smooth compared to ionic solids such as MgO or oxides such as Ti02. Thus, one expects very small barriers to translation across such a surface and diffusion that is very different from the standard descriptions of activated jumps over energy barriers with height equal to an activation energy Ea. It would appear that the best description of diffusive motion in monolayers on such surfaces might be based on

455 a two-dimensional gas approach, slightly modified to take account of the effects of the small variations in gas-solid energy as a molecule translates across the surface [41]. In fact, even for the rougher surfaces presented by adsorbents with atomic structure such as those presented by model amorphous or crystalline oxides, the barriers to translation of a molecule across the surface tend to be small enough that a description of diffusive motion in such systems in terms of jump motion should be carried out with caution. For the highly irregular gas-solid potential functions presented by heterogeneous surfaces, one has a wide range of "activation energy" in a complicated surface geometry that leads to trajectories that exhibit frequent changes in direction and which do not necessarily pass near the bottoms of the saddle points in the energy surfaces [42]. Many of these points will be illustrated by the simulation results that will be presented here. Another group of dynamical simulations involves the diffusive motion of isolated molecules over an adsorbing surface [43, 44]. Such systems, which can be studied experimentally by techniques such as field emission microscopy, are of course primarily probes of the gas-solid potential function, which for molecules with rotational degree of freedom and/or internal flexibility (such as short hydrocarbon chains) can be quite complex. It should be noted that molecular trajectories can be used to determine a variety of dynamical quantities, some of which are essentially

!

+'!

mo r 99

N I 98 97 9 I~

4

,o~ N, ,O4hn

,mt=~n~

.........

'%-'";

io4 r ,o31., ,oz.'-

,o,F

I00(~

98

,

el

,

, IAI'

5

97

,

'

I0

15

5

IOTr

,o~ ,o~

.

,o' '~,o ;n r

r I ,

moF 2

I, ,

I 5 I0 timesteps/lO00

io3r io21_ ,o, I_ 99~)

_.11.

,I

9%

5

I

I0

109F

,o8~ ,oT~

,oe~ I ~,11- , 5

8

' 5 i0 ~ timest eps/lO00

losr ,o41_ ,o31_ I01~) lit

I0

9

u_,,

" ' 5 I0 timeslel~/lO00

-

F i g u r e 2 Number of atoms in a dense N2 monolayer on the graphite basal plane at 73.6 K as a function of time. Note that 1000 timesteps = 2.2 picosecond. The coverages, relative to commensurate, are 1.04, 1.10, 1.25, 1.35, 1.50, 1.75, 2.00 and 2.50 for the indexes on the panels increasing from 2 to 9. (Panel 1 corresponds to the commensurate layer - no interlayer transfers were observed for this case during the time of observation.

456

inaccessible to experiment. The rate of transfer of molecules in and out of a monolayer film can be monitored- for example, see Figure 2, which shows the time dependence of the number of molecules in several monolayer and submonolayer films of nitrogen adsorbed on the graphite basal plane at 73.6 K [45]. The fluctuations in the number of molecules in the first layer can be due either to transfer from (or to) the gas phase or from (or to) the second layer. Figure 2 indicates that the jumps occur rather infrequently on a molecular time-scale with a frequency on the order of a few per picosecond, depending upon coverage. The jumps themselves are completed very quickly on this time scale. Information of this kind is of interest to the determination of rates of adsorption and desorption, but there appears to be little use of simulations in this connection to date. Another type of dynamics that can be simulated is the rotational motion of polyatomic molecules in an adsorbed film. In this case, angular velocity and reorientational time correlation functions can be evaluated (and compared with experiment- see [40]). In this case, one must

1.0---

_

_

-

~= 0.6

9

-

I Z

~

O.

3 4 I

(~o

1.0

I

I

2.0

I

I

3.0

4.0

TIME/pS

1.0 O.6

O.

'

0.0

i.0

2.0

12

5.0

4.0

TIME/ps

1.0 o.6

-O. L Vl

(3.0

!

!

I

I.O 2.0 TIME/ps

I

!

3.0

!

J

4.0

F i g u r e 3 Time-correlation functions for N2 on the graphite basal plane at low temperature. Part (a) shows the in-plane tcfs < cos(m$r > for m = 1 - 4 (top curve to bottom) for a commensurate monolayer at 18.8 K; part (b) shows the same functions for 41.4 K; and part (c) shows angular velocity time correlation functions for 41.4 K - the solid line is for out-of-plane tilting motion and the dashed line is for the in-plane reorientation.

457

realize that in-plane and out-of-plane molecular reorientation can be very different from each other due to the relatively large torque which will, for linear molecules, usually have the effect of constraining them to lie mostly parallel to the surface. In the case of in-plane rotation, the torques are due to the variations in gas-solid potential energy as a molecule moves parallel to the surface and to the torques arising from interactions with other adsorbate molecules in the same layer as the rotating species. Simulations of motion in the dense nitrogen monolayer [46, 47, 48, 49] gave the curves in Figure 3 which show in-plane libration at very low temperature (usually, a damped more-or-less harmonic oscillation as in the solid curve of Fig. 3c). At high temperatures, the in-plane molecular reorientation may show an approximation to rotational diffusion. In this case, theory gives exponential time-decays of the angular correlation functions gin(t)=< cos(m~r >, where ~ r the change in the in-plane angle of the molecular axis in time t and m is an integer. In fact, the solution of the two-dimensional rotational diffusion equation yields gin(l) = exp [-m2~t], where ~ is a rotational diffusion constant. In this limit, the angular velocity time-correlation function has been assumed to have a negligibly short life-time compared to the decay time of the orientational correlation function. Note that the in-plane rotational motion can be primarily diffusive even when the out-of-plane reorientation is a damped harmonic torsional oscillation. Part (a) of Figure 3 illustrates the in-plane orientational behavior for a fairly dense (15.7 .~l2 per molecule, or v/3 x v/3 commensurate) monolayer of N2 adsorbed on the graphite basal plane at a temperature below the in-plane orientational ordering transition which occurs for this system at ~_ 27 K and part (b) shows the behavior at a temperature above this transition. The almost constant time correlation functions (tcfs) of part (a) are indicative of small in-plane torsional oscillations, while the tcfs in part (b) correspond to in-plane motion which is more or less diffusional. Figures 4 and 5 show simulated orientational tcfs at 40.3 K for a N2 layer on graphite that is 5% denser than commensurate. In Fig. 4, rotational diffusion

t/ps 0.0 0.0

0.4

0.8 9

.

4..*

E -0.4

o C

-0.8

F i g u r e 4 In Cm(t)/m 2 for in-plane motion of N2 at 40.3 K on the graphite basal plane in a layer that is slightly compressed relative to commensurate.

458

theory is tested by plotting the in-plane orientational tcfs for for m = l to m = 4 on a scale that should give a single straight line for all m if the theory is obeyed [50]. Deviations from linearity at short time are due to the persistence of angular velocity correlations. It is probable that the

t/ps 0 0.00

,

,

1

2

I,

i

3 ,.,

i

z __

,,p+ .,...

:

-0.04

c t_.,,,.J

-0.08

I

-0.12

F i g u r e 5 In 7>~(t)/l(l + 1) for out-of-plane reorientational motion in the system of Fig. 4. l = l corresponds to the lowest curve with the other curves in ascending order for l=2 - 4 and ~t(t) denotes a tcf of t h e / ' t h spherical harmonic function of the out-of-plane angle.

1.0

0.8

~

~

T

~6 r

0

o.4 ;7*

0.2 O0

I

2 TIME (PICOSEC)

3

F i g u r e 6 Reorientational tcfs for monolayer N2 on the basal plane of graphite. Temperatures are shown on the Figure. The nearly constant curve is for out-of-plane (tilting) motion, and the decaying curves are for in-plane reorientation.

459

differences shown in Fig. 4 between the observed and the predicted slopes at long time are due to the fact that reorientation is no longer strictly in-plane at this temperature and surface density, but has become three-dimensional. This idea is supported by the out-of-plane reorientational tcfs that are shown in Figure 5. The decays of these functions, which should be linear and identical for all l if the out-of-plane reorientation is diffusional, indicate torsional oscillation at short times and non-diffusional but significant exponential decays at long time. Figure 6 shows the in-plane orientational tcfs for molecules in the commensurate N2 monolayer on graphite as a function of temperature for the case m = l . Also shown is the out-of-plane reorientational tcf. Clearly, out-of-plane motion is minimal for all the temperatures considered, which range up to the bulk N2 boiling point, even though the in-plane motion is becoming ever more rapid as temperature increases. It should be noted here that .,N2 on graphite exhibits two-dimensional melting at a temperature that depends strongly upon coverage, but is close to 48 K for an isolated commensurate patch on the surface [51]. The conclusions to be drawn from the simulation studies of rotational motion in dense N2 monolayers on the graphite basal plane are that orientational ordering is almost complete at sufficiently low temperature, and that the in-plane reorientation is considerably less hindered than the out-of-plane at higher temperatures. In fact, the out-of-plane reorientational motion is essentially torsional oscillation that is induced by the strong gas-solid forces on a N~ molecule on the graphite surface (or more accurately, strong torques). The torques are much weaker in the in-plane direction because of the flatness of the basal plane. Thus, hindrance in that degree of freedom is due primarily to interactions with neighboring N2 molecules in the adsorbed film. Although simulations of orientational dynamics for other adsorbates on other surfaces are quite limited, it is not unreasonable to assume that this behavior of the N2/graphite system might be characteristic of small non-spherical molecules on relatively fiat surfaces. Furthermore, persistent velocity correlations can have a significant effect upon the relevant tcfs in the translation problem as well the rotational, as we will see in the following.

2

D i f f u s i o n on n e a r l y flat s u r f a c e s

An obvious but important feature of diffusional motion in an adsorbed film is that it is extremely anisotropic, with motion perpendicular to the surface tending to be mostly harmonic oscillation, interrupted by occasional desorption events a n d / o r collisional changes in the oscillatory motion. Thus, the three dimensional molecular motion in real physisorption systems must be treated as two distinctly different types. One expects that the mean square displacements (msds) < ~z2(t) > and < 5p2(t) > = < ~x2(t) + ~y~(t) > as well as the corresponding velocity tcfs, will behave quite differently. Here, z is the coordinate perpendicular to the surface and x, y (and p) are parallel to it. The brackets here denote an average over all (or over a specific subset - see below) of particles and time-origins and 5z2(t) is (z(t) - z(0)) 2, the square of the displacement in the z direction in time t. The tcf and the displacement are related by the equation [52]: < t~z2(t) > - 2

( t - r) < v'._(O)vz(,) > dr

(1)

Analogous equation hold for the x and y msds and their velocity tcfs. If the molecules that give rise to the statistical quantities in eq. 1 are unbounded in their z-displacements, one has the

460

the well known long-time limiting form: ~ m < ~z2(t) > = 2t t----* OO

jo

< ~z(0)~(~)>

d,

(2)

and Dz, the diffusion constant for motion in the z-direction is thus Dz =

< v~(O)vz(,)

(3)

> d~

We define:

Cz(t) = < vz(O)vz(r) >

(4)

At short times, a power series expansion for the velocity tcf gives [52]:

Cz(t) = k T

< F~ > t 2

m

m2

2 + .....

(5)

where m is the molecular mass and Fz is the z-component of the force on a molecule. When this expression is substituted in eq. 2, one finds that the short-time msd is given by: < ~z2(t)

>=

kTt2

< F~ > t 4

~

m

m2

12 +

...

""

(6)

This argument shows that simulated msds must be quadratic functions of time initially, with a curvature that depends only upon the mean square velocity of a particle at temperature T with mass m, and not upon the density or forces acting in the adsorbed film. At the other limit of time, diffusive motion yields msds that are linear in time with slopes that give the diffusion 6

2

0

I

0

0.5

~~

1

I

I

I

I

I

1.5

2

2.5

3

3.5

4

time (picosecond) F i g u r e 7 The solid curve shows the mean square displacement for a system with an exponentially decaying velocity tcf with the decay time r = 1 picosecond. The straight dashed line is for strictly diffusional motion with diffusion constant ~ = 1 picosecond-:. The value of (kT/m) has been set equal to (1 A/picosecond) 2 in this case.

461

constants I)z, T~z and I)y. The approach to linearity in the msds is governed by the decay time of the velocity tcf, as can be seen explicitly by supposing for the moment that Cz(t) = (kT/m)exp[-t/Tx], where 7"z is the correlation time for the tcf Cx(t) and is equal to ~-1 where ~z is the diffusion constant. An integration of eq. 1 for this case gives:

kT < &~(t) >= --~-~[t + T~(~xp(-t/T~)- 1)]

(7)

?72

This msd has the correct initial (quadratic) time dependence but becomes linear in t only for times much greater than the decay time of Cx(t). The approach to linearity given by eq. (7) is shown explicitly in Figure 7. Simulated translational tcfs for dense physisorbed films are far from exponential function of time, as will be shown shortly for mono- and bi-layer films of oxygen on the graphite basal plane. Nevertheless, there is a class of systems that might be expected to conform to this simple expression. The physical basis goes back to the kinetic theory of dilute three-dimensional gases. A molecule in such a gas undergoes a series of binary collisions with, on average, a frequency vcou. It is reasonable to assume that a molecular velocity is constant until a collision occurs and that the velocity vector of the molecule after a collision is completely uncorrelated with its previous value. Since the fraction of molecules which have not yet collided in a time interval t is just exp[-~,coU t], the decay of the velocity tcf is given by this expression. If one asks whether there is an analogous situation in adsorbed films, the answer is a qualified "yes" if one considers a dilute film at high temperature on a very flat surface. Although the out-of-plane motion can be essentially vibrational, the in-plane motion could be well approximated by a freely translating two-dimensional gas. In this case, one might begin with the simple two-dimensional version of the collisional expression for the in-plane diffusion coefficient T~ll which is: 1 v,I = ~

< c>

(8)

where < c > is the mean speed = (8kT/rcm) 1/2 and )~ is the mean free path in two-dimensions and approximately equals 1/2v/2ap, with a equal to a hard disk diameter and p equal to the density in mole/cm 2. Simulations of a system that might fit this simple picture have been reported by two groups [53, 54] who have studied methane at room temperature in slit pores with parallel walls made up of the graphite basal plane. However, the point of view taken was rather different in the two studies. In [54], the total flux of molecules down the pore and the associated msds were evaluated after considering carefully the boundary conditions that controlled the collisions of the gas phase molecules with the pore walls. (This was a non-trivial problem because the gas-wall interaction had been assumed to be that for a perfectly fiat wall, thus yielding specular reflection and a physically unlikely flow mechanism in the absence of extra assumptions for the boundary conditions.) In [53], the msds were evaluated separately for the molecules in the monolayer on the wall and for those molecules in the interior of the pore. The molecules that jumped between the two regions and those in the central volume were excluded from the calculation. In fact, the local density within the pore shown in Figure 8 indicates that nearly all the gas in these systems was in fact in the monolayer. In [62], the gas-wall interaction was taken to have the periodic variation with position along the wall that has been used in the most accurate representations of the energy and that had the

462

advantage of spoiling the specular nature of a gas-wall collisions. In addition, two other systems were considered in which the wall had been roughened by adding "sulfur" atoms amounting to 1/13 (in a v/7 x ~/7 lattice) and 1/7 (in a 2x2 lattice) of the wall atoms. The variation in minimum methane-pore wall energy as an atom translates parallel to the wall is small for the non-sulfided surface, amounting to -'~ 400 joule/mole, but is "~ 2100 joule/mole for the dilute sulfided case and _ 2900 joule/mole for the more concentrated example. Apparently, these variations were insufficient to affect the surface diffusion process very much. In spite of the differences in outlook of the studies in [53] and [54], the self-diffusion coefficients obtained were also numerically similar over most of the range of methane loading. Since only the motions of the molecules in the monolayer films were considered in [53], the displacements in this case are "surface diffusion". The variation of the in-plane T~ with surface coverage ranging up to the approximate monolayer of 8 x 10 -1~ mole/cm 2 is shown in Figure 9 and is compared with the curve calculated from the crude two-dimensional coUisional expression given in eq. 8. The agreement is as good as one could expect from such a simple treatment. Note that the rapid variation of T~ with coverage is ascribed to the decrease in the mean free path of the adsorbate 60

9

,

1.65 1 0 - 1 0 m o l e / e r a 2 ...... . . . . . ......

3.30. I 0 - 1 0 m o l c / c m 2 4.95.10- l O m o l c / c m 2 6.6010- lOmole/cm2 . 5.

-

i o

[~ []

.

cr

.<

"~ ~o

o

-5

0

5

F i g u r e 8 Local density plots for various surface coverages of methane at 300 K in a slit pore with sulfided graphite walls. The separation of the two peaks at the opposing walls is 7.6/1. molecules in the two-dimensional gas on the surface as the density of this gas increases. The question of what happens in the limit of very low density is left unanswered, although one anticipates that the surface roughness will be the controlling factor in determining 79 in this limit. In addition to the simulation at 300 K, the methane/graphite system was studied at the

463

24-

.

~

lpr=Phit'm

~J a. OI

u 12

-,,,'....x,. X

_,,,

,,

.

.

.

.

.

0.0

I

.

.

--

3.5

7.0

Coverage x 10m" ( m o l e / c : r n ~

F i g u r e 9 Surface diffusivities for methane at 300 K in a slit pore with pure graphite walls, and with walls that contain model sulfur atoms in lattices with spacings that are v/7 • x/'7 or 2 • 2 relative to the underlying graphite lattice. Also shown is the self-diffusion constant evaluated for a crude version of the two-dimensional hard-disk gas.

10

.

.

.

.

.

.

,

.~ \ \

\\ -

9

ra

t

.

.

.

.

.

.

.

.

z P = 14.8,;(,/Tx 4T) z p = 14.s~(2 ~ 2)

....

z~, = u.1,~(,/T~

,/?)

\

9,, \

Ol

9. 9

5

"

.

. .....

% %%

t

"~

t

%

X

"~., "*%%

0

.

0

.

.

.

.

!

..

,

4.5

=

9

(]overage x 101~ m o | e / c m 2 F i g u r e 10 Same as Figure 9, but at 150 K. Here, the dependence of the surface T) upon wall separation ZP in the slit pore is shown together with the dependence upon surface methane density. (The separation of the peak methane densities in Figure 8 corresponds to a pore wall separation of 14.8 ~t).

464

considerably lower temperature of 150 K. Such adsorbed films can still be treated as twodimensional gases rather than surface liquids, and one can see from Figure 10 that a) the results are not very sensitive to the spacing of the two walls that make up the slit pore; b) the rapid decrease of :D with increasing coverage is quite similar to that found at 300 K; c) the change of D with temperature at fixed coverage is not large, and is of the correct order of magnitude to fit the simple collisional kinetic theory result of a v/T dependence. We now consider surface diffusion in very thin films at much lower temperatures and considerably higher coverages than for the methane/room temperature case. In particular, simulations of oxygen adsorbed on the graphite basal plane at temperatures from 50 - 60 K in monolayer and bilayer films will be discussed. [55]. The oxygen was modeled as diatomic molecules made up of Lennard-Jones interaction sites separated by the bond length so that orientational changes were accounted for in the calculations. Molecular dynamics was utilized to evaluate thermodynamic properties as well as the velocity tcfs for the molecules in each layer separately - only the molecules which remained in their designated layer for the duration of the calculation were included. In principle, in-plane diffusion constants could be evaluated from the time-integrals of the tcfs. The first point to be noted is that the oxygen at low temperature forms an ordered monolayer with all molecules standing nearly perpendicular to the surface with a density corresponding to 9.3 )12 per molecule. At slightly below 60 K., the monolayer melts to an area of 12.5 ~2 per molecule, with a considerable loss in orientational and in-plane translational order. The in-plane velocity tcfs for the monolayer molecules in the monolayer and bilayer films are

....

55K

.......... 6 0 K 0.5 tcf

-0.5~0

I

2 time

3

4

(picosec)

F i g u r e 11 In-plane velocity tcfs for an oxygen monolayer on the graphite basal plane at three temperatures. plotted in Figures 11 and 12 for several temperatures in this region. It is evident that they are far from the exponential decays expected for the binary collision mechanism. Their initial decays are quadratic, as required by theory, and they all show a negative dip at times less than one picosecond. This dip is associated with a reversal of the initial direction of motion due

465

to collisions or equivalently, to the strong forces that tend to enclose the adsorbed molecules in "cages" under these conditions. Of course, in a calculation of the diffusion coefficient by integration, the positive contribution of the short-time decay of the velocity tcf will be partially canceled by the negative dip, producing a small D which is rendered more uncertain by being the resultant of relatively large positive and negative contributions. In fact, the areas under

I-

I

0

....

I

60K

.

~176

--0.

0

I

2.

3

4

time (picosec)

F i g u r e 12 Same as Figure 11, but for the molecules in the first layer of a bilayer film. the curves in Figs. 11 and 12 are all zero to within the considerable uncertainty in the long-time parts of the tcfs. The alternative route of evaluating the msds gives much better estimates of D for the in-plane motion. Such curves are shown in Figures 13 and 14. Evaluation of these slopes yield D of ~_ 0.08 ~2/picosec. at 65 K. for both systems (l~i2/picosec. = 10 -4 cm2/sec.). At 60 K., D ~ 0.07 and 0.03 ~12/picosec. for the monolayer and the first layer in a bilayer film respectively. It appears that the presence of a second layer over the monolayer can impede the diffusive motion of the molecules below; this effect is associated with an increase in the density of the monolayer when it is covered by a second layer [55]. In all cases considered here, the diffusion constants are much smaller than those obtained for the two-dimensional methane gas discussed above, where D is of the order of 10 ~12 per picosec, for surface density equivalent to roughly 1/2 monolayer. The N2/graphite system has also been extensively investigated [42, 56]. In addition to the low-temperature reorientational results discussed above, translational motion in the N2 monolayer has been simulated over a range of higher temperatures. For example, msds for molecules in a complete N2 monolayer are plotted in Figure 15. The curves for T<70 K are as expected for a solid film in which the molecules are vibrating around lattice points, but the diffusive curve obtained for 77 K indicates that melting has occurred, probably at a temperature very near to 77 K, consistent with other simulation data. The in-plane diffusion constant obtained from the slope of the curve at 77 K is 0.21 ~t2 per picosecond.

466

~=~

=

.=

0.3

I

i

!

i

i

I

..........

T=65 K

J

-~

0.25

'-

0.2

...............ii iii ::: :

=

...:.::..:::".:;:.."

0.15

,=,=

.~ "G

..*~ssSS

0.1

9

** ,,.,

...::::'" .3"

f_

=

..-"""T=60 K

..- .......

q,)

0.05

T--5~

II0

9 E

0

9

0

~'

I

I

I

I

I

i

I

0.5

1

1.5

2

2.5

3

3.5

4

time(picosecond) Figure 13 Mean square in-plane displacements in reduced units for the system of Figure 11. Reduced distance = distance i n / i divided by 2.46 ~.

|

.,=..

0.3

0.25 ~176

~,

0.2

=

0.15

9-~

0.1

&

o.os

**.~ o~176 ~176

T=65 K 9 .......

,,.,.,....,.,... ,,.,o - ' ' " " " ~

.......:-::::i: ............. ..: ::':::::';""":1"---60K

T=55 K

m

~ E

..

**.~ o~ ~176

0 0

0.5

1

1.5 2 2.5 time(picosecond)

3

3.5

4

Figure 14 Mean square in-plane displacements in reduced units for the system of Figure 12.

467

Another system where diffusion on a flat surface has been simulated using molecular dynamics is that of benzene on graphite over the range 85 - 298 K [57]. We will discuss here some results obtained for the complete monolayers in this system. To begin, note that the simulations are in good agreement.with thermodynamic experiments but have yielded information concerning

0.3

T= 7 7 ~

~J

=.

0.2

J

~.T= 67o

~

O.l

~T-

57 ~

L T - 48" e~ u



f 0

~.

I

1

I

!

2

5

4

time(picosecond) F i g u r e 15 Mean square in-plane displacements in nitrogen monolayers at several temperatures. Also shown is a single curve for the displacements perpendicular to the surface (denoted by < z 2 >) - it is essentially independent of T in this range). The same reduced units for distance are used here as in Figs. 13 and 14. molecular orientation and time correlation functions that have not yet been measured in detail. It is known that the monolayer films form two-dimensional solids at T < 140 K. At 85 K, these films are highly ordered, forming a v/7 • v/7 lattice on the basal plane substrate with all molecules lying nearly flat on the surface. As temperature rises, the adsorbed benzenes begin to reorient into on-edge (perpendicular) configurations. The plots of two-dimensional density as a function of distance from the graphite plane in Figure 16 show how the molecular centers move away from the surface as molecules change from flat to on-edge orientation, giving double-peaked density distributions with an on-edge fraction that increases as temperature increases.. Since the area occupied by a reoriented molecule is significantly less than that for one lying flat, there is more space available for surface self-diffusion. This is clearly shown in the slopes of the msd curves given in Figure 17 and in the resulting self-diffusion constants plotted versus temperature in Figure 18. One sees small, nearly temperature independent values for the two-dimensional solid and a rapid increase in D at temperatures greater than 150 K. It seems most likely that the increase is primarily due to the increase in the available space for diffusion with increasing

468 temperature even though the number of molecules in the monolayer is remaining nearly constant kinetic energy of the adsorbed film. Orientational tcfs were also evaluated for this system and provide support for this physical picture. For molecules such as benzene, one should consider two types of reorientation: around the molecular axis, which should be relatively free because of the high in-plane symmetry of the molecule; and reorientation of the axis, which describes the tilting motion mentioned above. Figures 19 and 20 show the tcfs for < c o s ~ ( t ) >, where r is the angle for reorientation either around the molecular axis (Figure 19) or of the axis (Figure 20). It is evident that motion around the axis is indeed much easier. More significant for the

|

|

...

l J01 I.

D

I t t I

T--85K I

r

w

T--140K .e

~-

2.00

.m

J

I

Lm

~

I I ,

z(X) -

l v

i

1o.00

u

z(~.)

|

"i '--

T=200K

,~,

I

J

IJI0'

T--250K

.m

A I t_

4 j

L~!

" k ~

'-

UO

.m m n

"N

I,~

m

,

' ...... XO,W

T

1

~

zCk) F i g u r e 16 Local densities for benzene monolayers on the graphite basal plane plotted as a function of distance from the surface z in ~t for a range of temperature. The peak at small z corresponds to those molecules lying fiat on the surface, and the second part of the split peak corresponds to molecules standing on edge.

469

20.00 18.00 16.00

~-298

=

14.OO

E

~r-fs6 ~r=f25

12.00 , mm

I0.OO 8.OO

6.OO

E

4.00 ZOO 0.00

]

2.0o

0.OO

!

t

4.00

6.00

time(picosecond) F i g u r e 17 Mean square displacements for benzene molecules in the monolayer on graphite are shown for various temperatures.

0.75 _

_

I

1

i

O.5O

u 0 u

N ~

0.25

_i-.----!

!

!-

I00

200 T ((leg. K)

300

F i g u r e 18 In-plane self-diffusion constants for monolayer benzene obtained from the slopes of the lines in Figure 17 are plotted as a function of temperature.

470

present purposes, one can see that the tilting motion is nearly frozen at the lowest temperatures, but becomes quite extensive on a picosecond time-scale as T increases to 298 K. In comparing the :P values for monolayer benzene with those reported above for monolayer oxygen, one should first realize that the temperatures of the two systems are rather different, especially on a comparative scale that takes the two-dimensional critical temperatures as reference points. These are

0

.

0

I

~

" t 1

i: --

-1.2 -

~

--

C

-1.6

%.

"

'%

_

-

~

-

-

-2.0 -

'~

%

-2.4

-

%I

1 I

0

__

! 2 t (picosec)

I 3

I 4

F i g u r e 19 Orientational tcfs for benzene on graphite at the same series of temperatures as in Figures 17 and 18. These curves illustrate the rapid reorientation of a vector fixed in the plane of the benzene molecule. i

"

I

~

-0.4

I~176

"~ I

!--

i

1

/

e

~176

"'~.

,,

~

/

x

0BI

"'""

-,.or /

-.,

70

"". -1 I

2

3

4

I (picosec)

F i g u r e 20 Same as Figure 18, but here the vector is the molecular axis perpendicular to the molecular plane, so the motion characterized here is the tilting motion of the adsorbed benzene molecules.

471

approximately 205 K and 60 K for benzene and oxygen, respectively. Thus, much of the benzene data shown is supercritical and that for oxygen is subcritical. (Of course, the densities of the monolayers are roughly three times the critical density, so that neither system is actually near the critical region at any temperature.) In any case, over most of the temperature range, the benzene diffusion constants are much larger than those for the oxygen layers. This is probably due to the relative lack of space for diffusional motion in the oxygen case.

3

D i f f u s i o n o n s u r f a c e s w i t h s t e p s or g r o o v e s

In the previous section, we have been discussing diffusion on surfaces that are geometrically quite flat. Even in those cases where atomic structure plays a significant role in determining the gassolid potential energy surface, the roughness is only of the order of the atomic size and is widely distributed over the surface. We now take up an important class of adsorption systems that are characterized by strong barriers to translation in specific regions of an otherwise flat surface. Examples include model pores with sharply defined constrictions, a simple version of which has been studied by Demi and Nicholson [58, 59], who considered atoms in spherical pores connected by short, straight-walled narrow cylinders. We do not give a detailed discussion of this interesting study here on the grounds that the simulations of the flux of molecules down these tubes included both surface and volume flow and is consequently not clearly relevant to surface diffusion, which is the primary subject of the present review. Another type of heterogeneity is encountered when one has surfaces with atomic-scale steps or even with straight-walled grooves cut into the surface. There has been considerable interest in the diffusion of strongly bound gases on metallic stepped surface. Uebling and Gomer have studied diffusion by Monte Carlo simulation of a lattice gas on stepped surfaces with various assumptions concerning the interaction of the gas atoms at the steps [60, 61]. These and similar studies have been reviewed by Gomer [5]. Simulations of the surface diffusion of krypton on grooved and stepped surfaces at temperatures ranging from 90- 150 K (boiling point of bulk Kr = 120 K) [62, 63, 64, 65] will be discussed here. These simulations are not subject to the limitations of the Monte Carlo algorithm as applied to the lattice gas model, but the main features of this process are qualitatively the same as those found by Uebing and Gomer, when the temperatures and surface coverages are comparable. The nature of the diffusion on these surfaces is determined by the same features as for other adsorption systems: the gas-solid potential and its variations across the surface relative to kT plus, for finite coverage, the adsorbate-adsorbate interactions. Figure 21 shows a side view of a model stepped surface [62] with a nominal monolayer of Kr adsorbed on it at 110 K. The straight lines denote graphite planes separated by 3.4 ~i, and the terraces are 25.5 ~t wide and 44.3 ~1 deep (in the direction perpendicular to the plane of the paper). The trajectories of the adsorbed krypton atoms generated in an isokinetic molecular dynamics simulation lasting for 71 picoseconds are also shown. The thick black regions are overlapping trajectories for the atoms in the monolayer, and the irregular lines are trajectories for those few atoms that have left the monolayer at 110 K. This graphitic adsorbent shows strikingly large variations in the minimum gas-solid potential energy (i.e., the adsorption energy) as an adsorbed atom moves perpendicular to the steps, as is illustrated in Figure 22. Clearly, there is a very large barrier to free translation

472 present at the top of each step equal to 6.3 Kjoule/mole for the Kr/stepped graphite system. This is a consequence of the model used for the gas-solid interactions of this non-metallic substrate in which it is assumed that this energy is a pair-wise sum of the interactions of a Kr atom with sites in the solid that have been taken to be carbon atoms with interaction parameters adjusted to give good agreement between the calculated adsorption energy and the voluminous experimental data available for the Kr/graphite system. As a result, an atom sitting at the top edge of the step sees what amounts to only half a surface, so that its gas-solid interaction energy

_ ~:--

F i g u r e 21 Side view of the trajectories followed by 288 Kr atoms on a stepped surface over a period of 70 picoseconds at 110 K. The surface has a step height of 6.8 ~l, a terrace width of 25.5 ~1 and a dimension perpendicular to the page of 44.3 ~i. (Periodic boundary conditions are imposed in the two directions parallel to the surface.) The Kr density of 0.0636 a t o m / ~ 2 is that of the monolayer on the flat surface ([67], but a few atoms have been thermally promoted into the second layer. ,,6OO

i -,,t20(

S S !

,.-~SG 0

51

,,

102

Y (A) F i g u r e 22 The potential energy of an adsorbed Kr atom on the surface shown in Fig. 21 is shown as a function of position relative to the steps (Y-coordinate). The parameters and potentim model used were drawn from the analogous calculations for Kr on the infinite flat graphite surface [42].

473

is roughly half the value of-11.2 Kjoule/mole for an atom on the infinite flat surface. (Note that other models used for the interactions of gases with metallic surfaces can yield no barriers or even wells at the tops of steps rather than the barriers obtained here [66].) Another feature of this model of the surface is the set of deep but narrow potential wells for atoms at the bottoms of the steps adjacent to the step risers. This amounts to 2.6 Kjoule/mole for the model Kr/stepped graphite system, which is quite large compared to kT (= 0.9 Kjoule/mole at 110 K). The energy variations have a large effect on the diffusion and the local densities of krypton adsorbed on this surface. The local densities are shown in Figure 23 as a function of Y, the coordinate perpendicular to the steps, for two temperatures and two Kr surface coverages (1/2 and 1 monolayer, based on the known monolayer density of 0.0626 atom/A 2 for Kr on the fiat graphite surface [42]). It is evident that the strong variations in adsorption energy have induced one-dimensional order in these layers. Although a localized (in Y and in Z, the coordinate perpendicular to the surface) line of atoms forms in the deep wells, the ordering in the 50

.

50

'

A

A

.

,

~

Q,I ,.,.. O

B

O

a25

25

>(

X

>~

0

51

o

102

0

51

102

5O

5O

C D _.o

o

o

0

o

o 25

o 25

x

X

w

m al

0

0

51 Y (A)

102

0

51

102

Y (A)

F i g u r e 23 The local density of the first layer Kr atoms on stepped graphite is shown as a function of the Y-coordinate defined in Figure 22. Densities are shown for the monolayer at 110 K (part (A)- see also Fig. 21) and 150 K (B), and for 144 adsorbed atoms (equal to 1/2 layer) at 110 K (C) and 150 K (D).

474

Y-direction actually extends from the bottoms of the steps out into the terrace, gradually decaying as the atom moves away from the step. The potential energy variations and the local densities of Kr on grooved graphitic surfaces [64, 65, 66] are much the same as for the stepped surface case, except that the symmetry of a rectangular groove is reflected in both the energies and the densities. This is illustrated in Figure 24 which shows the local density for Kr at 110 K. on a surface with a series of grooves of depth ~ 7 .~t, separation = 34 A (referred to as step width) and groove width = 68 .;t. Strong one-dimensional ordering is seen for those atoms confined in the grooves, even though the temperature is above that for freezing of the monolayer

0.300.25~- 0.20 ;

i

o.os 0.00 ~,

9- 6 0 - 4 S - 3 0 - t S

0

IS

[

30

9

45

60

,r (1)

F i g u r e 24 Local density is plotted in the top part of the Figure as a function of Y (in ~) for monolayer Kr on a grooved graphitic surface at 110 K. The lower part of the Figure shows the trajectories of the adsorbed atoms (for 70 picoseconds), both as a side view and as a top view. The top view shows the localized fines of atoms adsorbed at the bottoms of the groove edges - this corresponds to the highest pair of peaks in the local density plot. The relatively vacant regions near the tops of the edges are visible in the top view of the trajectories and in the local density.

475

on the flat graphite surface. Weak ordering is also seen in the layer adsorbed on the steps between the grooves. Of course, the degree of ordering increases rapidly as the total surface coverage increases, as can be seen in the trajectory plots for the simulations of the systems with coverages of up to 2 1/2 nominal monolayers. As temperature decreases, the order also increases and indeed, for high coverage and low temperature, one observes two-dimensional solid-like order in each of the grooves [65]. We now turn to a discussion of the time-dependent properties of the atoms adsorbed on such surfaces. Although we will limit the actual results presented to those for a nominal monolayer, a number of general features should be mentioned at the outset. The first point is that the large barriers to motion in the Y-direction will cause surface diffusion to be quite anisotropic. Of course, in the simulation it is not difficult to evaluate mean square displacements in each of the three directions separately. At the temperatures and coverage of interest here, motion in the Z-direction is nearly pure vibration except for the occasional jumps at the steps or the groove walls. These jumps differ greatly from the small-step translation that occurs in the dense layers on the fiat regions of the surface. Motion on the flat regions is likely to be diffusive, but should be quite different in the X-direction, where there are no barriers associated with the surface structure, and in the Y-direction, where the atoms encounter both deep (trapping) wells and repulsive barriers. In fact, at long times, the mean square Y-displacement should approach an almost constant value that is determined by the width of the step or the groove- any small increase over this limit will be due to those atoms that jump to the adjacent step or groove. It is therefore necessary to evaluate X and Y displacements separately and in addition, the jump rates from one confined area to the next should be evaluated if one wishes to obtain the full description of motion on these surfaces. The rates of transfer of atoms by jumping from one region to another can be monitored rather easily by evaluating the number of atoms in a region (step or groove) as a function of time. Figure 25 shows this for a nominal monolayer at 110 K on the grooved surface of Fig. 24. The changes in the number of atoms on the steps given in this Figure determines the step--groove transfer rates. (In the equilibrium layer, the forward and the reverse rates determined by counting over a long time are, of course, equal.) The rate is equal to the total i

r

0 0

oJ

,

i

I

I

'

~

n

i

i

~TLj,

r

I

I

!

I.

I

Timesteps]l(~30

F i g u r e 25 Changes in the number of atoms on the step of a grooved surface at 110 K are shown as a function of time (1000 timesteps = 7.4 picoseconds). The system is monolayer Kr on a graphitic surface with a total step edge length of 88.6 ~l.

476

number of jumps divided by the observation time and by the total length of the groove edge, but can be only roughly estimated because of the rather small number of events observed in the time intervals of the simulations. Table I shows rates estimated in this way from data gathered for a 71 picosecond time interval for a total edge length of 88.6 A. The rates amount to 0.013, 0.018 and 0.018 for T = 90, 100,and 130 K., respectively (in units of atoms/~t/picosecond). Since the rates depend upon the density of atoms near the groove edges, the expected rate increase in going from 100 to 130 K appears to be hidden by a decrease in this density. Mean square displacements for atoms in the Kr monolayer were evaluated for those atoms on the steps, in the grooves but not at the step edges, and at the bottoms of the steps (denoted by "edge" atoms). The time dependences of the displacements in the X and Y directions are shown in Figs. 25 and 26 respectively. First, note that the X-displacements are not very sensitive

1.01

1

i

I.

130 K

...-

0.8

((~y)z)

0.6 0.4 0.2

0.5 1.0 time (picosec) 1

i

!.5

[

0.8

.=

90 K

<(Ay}2>

0.6

==

.=

0.4 0.2

0.5

1.0

1.5

time (picosec)

F i g u r e 26 The time dependence of the mean square displacements of first layer Kr atoms on the grooved graphite surface are shown here for two temperatures. The curves given are for atoms inside the groove (excluding the edge site) - solid lines, on the step - dashed fine, and on the edge sites - dotted fine. The displacements here are in the direction parallel to the grooves.

477

to the initial location of the diffusing atoms, at least for the step and groove atoms. This is reflected in the values of the self-diffusion constant T)x that were calculated from the slopes of the linear portions of the X-displacement plots - see Table I. The Y-displacements differ greatly, depending upon starting position. The deep potential wells occupied by the edge atoms means that this atomic motion is primarily vibrational, and corresponds to a constant mean square displacement which is reached after a rather short time interval. (At 130 K, it appears that a few atoms escape from their initial positions at the bottom of the step). The repulsive potential barriers encountered during motion in the Y-direction reduce the displacements of atoms that are on the steps, relative to those in the grooves. The apparent diffusive motion of these atoms in the Y-direction is misleading: the linear dependence of the root mean square Y-displacements upon time that are shown in Figs. 25 and 26 is for displacements that are still quite small compared to the sizes of the confining boxes (34 and 68 Jl for a step and a groove, respectively). Thus, the slopes of these lines give :Dr listed in Table I that are only apparent, valid for small rms displacements, with true long-time values that are nearly zero. In fact, a reasonable approximation to this situation would be one-dimensional diffusion in a box with impermeable walls. The many solutions of the diffusion equation [68, 69] (or of its mathematical relatives, the Schroedinger equation [70] or the heat conductivity equation [71]) include those for this boundary condition so that comparisons of the simulated and the theoretical msds for motion in the Y-direction for large time intervals should give a good test of the model. |

I

I /

l.Oi 130K

/ i

~

A::-"

F > o-6t"

/. :;-"

c'2' o.41 0.2

..-"" 0.5 '

0.6 ((~x

12)

|

i.0 '

|,

1.5 9

|,,

.~,~. -

90K

0.4

0.2

-

,

0.5 t i m e (pir

L

,

1.0

1.5 I

F i g u r e 27 Same as Figure 26, except that the direction of motion here is perpendicular to the groove direction.

478

Table 1: Self-diffusion constants for monolayer Kr on a grooved graphite surface (units of ~!2 per picosec) Groove atoms 9O K 130 K Step Atoms 9O K 130 K Edgeatoms 90 K 130 K

Z)x

1)y

0.16 0.34

0.17 0.34

0.15 0.31

0.13 0.29

0.15 0.23

_~ 0 "" 0

Of course, surface diffusion is also dependent on the transfer of atoms between layers or into the vapor phase. In simulations, it is possible to keep track of the locations of all atoms at all times. Thus, the mean square displacements shown in Figs. 26 and 27 are actually for atoms that start in the first layer on steps, in grooves or at edges, and that stay there for the duration of the calculation. Any atoms that change layers or change regions have been excluded. Figure 28 shows a limited comparison of mean square displacements when the layer-changing atoms have or have not been included. It is evident that diffusive motion is enhanced by atoms that have been promoted. They move much more quickly in the very dilute second layers present in these nominal monolayer systems, as one might expect. The importance of this factor will depend very much on coverage, temperature and the magnitude of the adsorption energy. The I

454-

I

I

I

!

I

454-

A

I

t

t

/

B

I

/

/....-':.....-" -. 3.63-

~ -

o(l: ..,,. A N 2.72<] v L81 -

,3.63 -

~....-"'

-

~.-"'"

-

.................... OC)

I

" 2

3

time (picosec)

4

Z" "3

~,

/.....-"

272-

S""

1.81--

"'"'"

,,.,..-'""""

5

o ~. o

i

2

1

3

time ( plcosec}

!

4

i

5

F i g u r e 28 Same as Figure 27, except that here, the displacements are compared for atoms that stay in layer one throughout (A) and for those atoms that start in layer one but can be promoted to higher layers during the time of the simulation (B). T = 140 K.

479

point to make here is that the simulation technique has the power to elucidate the various complications that arise in atomic diffusion on a realistic model surface. Results for diffusion on the stepped surfaces modeled are more limited, in part because the atomic motions on these surfaces are expected to be rather similar to those for the grooved surfaces. However, jump rates over the steps have been simulated [62] for several Kr coverages and temperatures. The data shown in Fig. 29 are for terraces of width 25.5 .~1. Coverages range from 1/2 to 3/2 layers and the jump rates are shown for three temperatures. At the lowest coverage, the rates are very small, indicative of jumps over the high barrier shown in Fig. 22 for atoms on the surface. At coverages of one and 1 1//2 layer, most of the jumps are due to atoms that are not on the surface but are in the second layer where the step-crossing barrier is much lower. Finally, the reason why the rate for 1 1//2 layers at 110 K differs from those at the same coverage but at 130 and 150 K is that the atoms must first break free of the solid-like structure that exists at 110 K but not at 130 or 150 K (for this coverage) before they can begin to translate to the adjacent step.

6oI 50 -

-'" J..." /'...."

40 0

ol

,.,.:.'.-"

L

I

..~ ....

2;"

30

20

10

0 !

,

!

144

288

~

[

432

N u m b e r of atoms

F i g u r e 29 Jump rates for Kr on a stepped graphitic surface with terrace widths = 25.5 .~1 and total step length = 178 .~. Rates are shown for three coverages (1/2, 1, and 1 1//2 layers), with the simulated values connected by straight lines. The solid line is for 110 K, the dotted line, for 130 K, and the dashed line, for 150 K.

480

4

Conclusions

In this paper, surface motion that is diffusive and, in some cases, proceeds by jumps, has been simulated for models of varying levels of complexity. All cases considered have been representations of actual adsorption systems. Because of the emphasis on surface motion in the first couple of layers on the adsorbing solid, there is relatively little discussion of diffusion in porous solids, where a number of issues such as the connectivity of the pore system or the effect of variable pore dimensions, especially when they give rise to constrictions with sizes comparable to molecular diameters, complicate the problem of relating simulation to experiment. Even so, the work on the dynamics of fluids in confined geometries is just beginning. There are almost no molecular dynamics studies of the important case of atoms on heterogeneous surfaces where the heterogeneity extends over the entire surface. Also, relatively little effort has been put into the study of the effect of non-spherical shape upon surface mobility. Still, the results obtained so far provide the beginnings of a picture of the way in which molecules move when they are on or near to the surface of a solid. The values obtained for diffusion constants are consistent with those when the three-dimensional case is modified into two dimensions. For example, the results for the high temperature two-dimensional "gas" of adsorbed methane can be accounted for in terms of a simple twodimensional kinetic theory argument. Furthermore, the dense, two-dimensional "liquids" such as N2, 02 and benzene show self-diffusion constants that are reasonable consistent with the bulk liquid results [72, 73, 74]. Here one has the problem that the diffusion constants for both the twoand the three-dimensional fluids are quite sensitive to density, and it is not clear how one should estimate corresponding densities in two and three dimensions. For example, the experimental 7) for bulk liquid benzene at 300 K is 0.2 ~2/picosec, but Fig. 18 indicates that this value is reached at 200 K for the monolayer fluid. However, one should realize that 200 K is nearly equal to the estimated 2D critical temperature of 205 K for benzene on graphite [57], whereas 300 K is well before the 3D critical temperature. Therefore, a careful scaling of density and temperature would be required before a quantitative comparison of the 3D and 2D diffusivities can be made. Similarly, the self-diffusion constants are known for the bulk liquid nitrogen, both from computer simulation [75, 76] and experiment [77]. A value of 0.33 )i2/picosec for the saturated bulk liquid at its normal boiling point is comparable with that of 0.21 Jl 2/picosec obtained from the slope of the 77 K curve in Figure 15 (and 0.08 )i2/picosec for the oxygen monolayer at 65 K - see Figure 13). Of course, these simulations are not without problems. A point that is often overlooked in work of this kind is that the assumption that the the substrate is rigid, providing only a potential field for the fluid adsorbate means there is no energy interchange between adsorbent and adsorbate. Within this model, different molecular dynamics algorithms allow one to evaluate trajectories with several different constraints. Two of those most often used are: constant total energy (of all the atoms); and constant total kinetic energy. In the first case, motion of an atom that results in a large change of its adsorption energy causes the entire sample to "cool"; i.e., any increase in potential energy is associated with a decrease in kinetic energy. Since this seems to be unrealistic (and possibly important, in a simulation of atomic time-dependence), the second or isokinetic algorithm is preferred. (Both give the same thermodynamic properties but well-defined differences in the fluctuations of the thermodynamics.) Nevertheless, it remains to be proved

481

that conservation of kinetic energy in adsorption systems is actually a good representation of the actual situation. It should also be noted that the assumption of a rigid, inert adsorbent means that thermal equilibrium in the adsorbed phase is achieved via energy-exchanging collisions between the adsorbate molecules only. This can produce questionable dynamical results at low coverages where such collisions are infrequent. Furthermore, the errors that might be produced depend upon the time-scale of the observation. (Imagine a particle in the process of passing over a feature in the gas-solid potential that corresponds to a barrier to diffusion. Is it in thermal equilibrium at all points along the trajectory? If not, what are the consequences relative to, say, an analysis via transition state theory?) To investigate these questions would necessitate a much larger simulation in which the adsorbent atoms are allowed to take part in the dynamics, if only to the extent of (anharmonic) vibrations that allow energy interchange between solid and gas. Such calculations are costly in terms of computer time, but are nevertheless highly desirable.

Acknowledgements Support for this review and for many of the results reported here was provided by grant DMR 9022681 of the N.S.F. Division of Material Research. Many of the results discussed were obtained by Alexei Vernov, Mary J. Bojan and V. Bhethanabotla while at Penn State. Not all of this work has been previously published and their contributions are greatly appreciated. Finally, permission to republish many of the figures was kindly granted by the authors and initial publishers. These include: American Chemical Society: Figures 1 [33], 2 [45], 8 [53], 21,22,23 [62] and 25 [65]. Taylor and Francis: Figures 3 [47] and 4,5 [48]. Engineering Foundation: Figures 6, 15 [46]. Materials Research Society: Figures 9, 10 [53]. VCR Verlagsgesellschaft mbH: Figures 25, 26, 27, 29.29 [66].

482

References [1] M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Oxford University Press, Oxford, 1987. [2] J. M. Haile, Molecular Dynamics Simulation, John Wiley and Sons, New York, 1992. [3] F. F. Abraham and G. White, J. Appl. Phys. 41 (1970) 1470. [4] Applications of the Monte Carlo Method in Statistical Physics, ed. K. Binder, Springer, Berlin, 1984; K. Binder and D. W. Heerman, Monte Carlo Methods in Statistical Physics - An Introduction, Springer, Berlin, 1988; Monte Carlo Methods in Statistical Physics, ed. K. Binder, Springer, Berlin, 1992. [5] R. Gomer, Rep. Prog. Phys. 53 (1990) 917. [6] V. P. Zhdanov, Elementary Processes on Solid Surfaces, Plenum, New York, 1991. [7] V. Limoge and J. L. Bouquet, Phys. Rev. Lett. 65 (1990) 1. [8] C. H. Mak, J. L. Brand, A. A. Deckert and S. M. George, J. Chem. Phys. 85 (1986) 1676; C. H. Mak, B. G. Koehler and S. M. George, Surf. Sci. 191 (1987) 108; J. L. Brand, A. A. Deckert and S. M. George, Surfl Sci. 194 (1988) 457; C. H. Mak, H. C. Andersen and S. M. George, J. Chem. Phys. 88 (1988) 4052. [9] J. L. Riccardo, M. Chade, V. Pereyra and G. Zgrablich, Langmuir 9 (1993) 2670; ibid, 8 (1992) 1518. [10] K. Sapag, F. Bulnes, J. L. Riccardo, V. Pereyra and G. Zgrablich, J. Phys. Cond. Matter A 5 (1993)223.

[11]

K. Sapag, V. Pereyra, J. L. Riccardo and G. Zgrablich, Surf. Sci. 295 (1993) 433.

[12] P. Argyrakis, A. Milchev, V. Pereyra and K. W. Kehr, Phys. Rev. E, in press. F. Bulnes, J. L. Riccardo, G. Zgrablich and V. Pereyra, Surf. Sci. 260 (1992) 304. [14] A. Pekalski and M. Ausloos, J. Chem. Phys. 100 (1994) 3175.

[15]

X.-P. Jiang and H. Metiu, J. Chem. Phys. 88 (1988) 1891.

[16]

I. Avramov and A. Milchev, Phys. Rev. E, 47 (1993) 2303.

[17] S. Havlin and D. Ben-Avraham, Adv. Phys., 36 (1987) 695. [18] S. SokoIowski, Phys. Rev. A 44 (1991) 3732.

[~9]

U. Heinbuch and J. Fischer, Phys. Rev. A 40 (1989) 1144.

[20]

S. A. Somers, K. Ganapathy Ayappa, A. V. McCormick and H. T. Davis, Adsorption 2 (1996) 33.

483 [21] I. Bitsanis, S. A. Somers, H. T. Davis and M. Tirrell, J. Chem. Phys. 93 (1990) 3427. [22] S. A. $omers and H. T. Davis, J. Chem. Phys. 96 (1992) 5389. [23] A. V. Klochko, E. M. Piotrovskaya and E. N. Brodskaya, Langmuir 12 (1996) 1578. [24] W.-J. Ma, L. K. Iyer, S. Vishveshwara, J. Koplik and J. R. Banavar, Phys. Rev. E 51

(~99~) 441. [25] J. J. Magda, M. Tirrell and H. T. Davis, J. Chem. Phys. 88 (1207) 1207. [26] J. Koplik, J. R. Banavar and J. F. Willemsen, Phys. Fluids A 1 (1989) 781. [27] G. Mo and F. Rosenberger, Phys. Rev. A 42 (1990) 4688. [28] M. Schoen, J. H. Cushman, D. J. Diestler and C. L. Rhykerd, Jr., J. Chem. Phys. 88 (1988) 1394. [29] D. J. Diestler, M. Schoen, A. W. Hertner and J. H. Cushman, J. Chem. Phys. 95 (1991) 5432. L. K. S. [30] M. Schoen, J. H. Cushman and D. J. Diestler, Mol. Phys. 81 (1994) 475. [31] W. Dong and H. Luo, Phys. aev. E 52 (1995) 801. [32] V. N. Burganos, J. Chem. Phys. 98 (1993) 2268. [33] A. A. Vernov, W. A. Steele and L. Abrams, J. Phys. Chem. 97 (1993) 7660 and references to earlier simulations given therein. [34] A. K. Jameson, C. J. Jameson and R. E. Gerald, II, J. Chem. Phys. 101 (1994) 1775. [35] S. Bandyopadhyay and S. Yashonath, J. Phys. Chem. 99 (1995)4286; S. Yashonath and P. Santikary, Mol. Phys. 78 (1993) 1. [36] Y. Nakazaki, N. Goto and T. Inui, J. Catalysis 136 (1992) 141. [37] R. L. June, A. T. Bell and D. N. Theodorou, J. Phys. Chem. 96 (1992) 1051. [38] G. Schrimpf, M. Schlenkrich, J. Brickman and P. Bopp, Phys. Chem. 96 (1992) 7404. [39] P. R. Van Tassel, S. A. Somers, H. T. Davis and A. A. McCormick, Chem. Eng. Sci. 49 (1994) 2979. [40] M. Bienfait, J. P. Coulomb and J. P. Palmari, Surf. Sci. 182 (1987) 557; M. Bienfait, in Dynamics of Molecular Crystals, ed. J. Lascombe, Elsevier, Amsterdam, 1987, p. 353. [41] D. L. Koch, J. Chem. Phys. 101 (1994) 4391. [42] W. A. Steele, Chem. Rev. 93 (1993) 2355. [43] K. Fichthorn, Adsorption 2 (1996) 77.

484 [44] D. Cohen and Y. Zeiri, Surf. Sci. 274 (1992) 173. [45] A. V. Vernov and W. A. Steele, Langmuir 2 (1986) 219. [46] A. V. Vernov and W. A. Steele, in Fundamentals of Adsorption, ed. A. I. Liapis, Engineering Foundation, New York, 1987, p. 611. [47] J. Talbot, D. J. Tildesley and W. A. Steele, Mol. Phys. 51 (1984) 1331. [48] R. M. Lynden-Bell, J. Talbot, D. J. Tildesley and W. A. Steele,, Mol. Phys. 54 (1985) 183. [49] A. V. Vernov and W. A. Steele, unpublished. [50] R. M. Lynden-Bell and W. A. Steele, J. Phys. Chem. 88 (1984) 6514; W. A. Steele, in Spectroscopy and Relaxation of Molecular Liquids, ed. D. Steele and J. Yarwood, Elsevier, Amsterdam, 1991. [51] Y. P. Joshi and D. J. Tildesley, Molec. Phys. 55 (1985) 999. [52] J. P. Hansen and I. R. McDonald, Theory of Simple Liquids, 2nd Ed., Academic Press, New York, 1986, sec. 7.2. [53] M. J. Bojan and W. A. Steele, Mat. Res. Soc. Symp. Proc. 290 (1993) 127; R. van Slooten, M. J. Bojan and W. A. Steele, Langmuir 10, 1994, 542. [54] R. F. CrackneU, D. Nicholson and K. E. Gubbins, J. Chem. Soc. Faraday Trans. 91 (1995) 1377. [55] V. R. Bhethanabotla and W. A. Steele, Phys. Rev. B 41 (1990) 9480. [56] W. A. Steele, Langmuir, (1996) 12 (1996) 145. [57] A. V. Vernov and W. A. Steele, Langmuir 7 (1991) 3110; ibid 7 (1991) 2817; in Proc. Fourth International Conference on Fundamentals of Adsorption, ed. M. Suzuki, Kodansha Publishers, Tokyo, 1993. [58] W. Demi, J. Chem. Phys. 95 (1991) 9242. [59] W. Demi and D. Nicholson, Langmuir 7 (1991) 2342. [60] C. Uebing and R. Gomer, J. Chem. Phys. 95 (1991) 7626, 7636, 7641, 7648. [61] C. Uebing and R. Gomer, Surf. Sci. 306 (1994) 419, 427. [62] M. J. Bojan and W. A. Steele, Langmuir 9 (1993) 2569. [63] M. J. Sojan and W. A. Steele, Surf. Sci. 199 (1988) L395. [64] M. J. Bojan and W. A. Steele, Langmuir ~ (1989) 625. [65] M. J. Bojan and W. Steele, Ber. Bunsenges Phys. Chem. 94 (1990) 300. [66] J. A. Serri, J. C. Tully and M. J. Cardillo, J. Chem. Phys. 79 (1983) 1530.

485 [67] N. D. Shrimpton, M. W. Cole and W. A. Steele, in Surface Properties of Layered Materials, ed. G. Benedek Kluwer, Dordrecht, Netherlands, 1991. [68] M. P. Allen and A. J. Masters, Mol. Phys. 79 (1993) 435. [69] J. Crank, The Mathematics of Diffusion, Clarendon Press, Oxford, 1979. [70] H, Margenau and G. M. Murphy, The Mathematics of Physics and Chemistry, Second Edition, D. van Nostrand, New York, 1956, Chap. 7. [71] H. S. Carslaw and J. C. Jaeger, Conduction of Heat in Solids, Clarendon Press, Oxford, 1959. [72] A. F. Collings and I. L. McLaughlin, J. Chem. Phys. 73, 1980, 3390. [73] C. Hoheisel, Theoretical Treatment of Liquids and Liquid Mixtures, Elsevier, .New York, 1993, sec. 8.2. [74] H. J. V. Tirrell and K. R. Harris, Diffusion in Liquids, Butterworth's, London, 1984, sec. 7.3. [75] P. S. Y. Cheung and J. G. Powles, Mol. Phys. 30 (1975) 921; ibid, Mol. Phys. 32 (1976) 1383. [76] J. Barojas, D. Levesque and B. Quentrec, Phys. Rev. A 7 (1973) 1092. [77] K. Krynicki, E. J. Rahkamaa and J. G. Powles, Mol. Phys. 28 (1974) 853.