- Email: [email protected]
Substrate response to atomic scattering from solid surfaces
J. Phys. Chem. Solids Vol. 48, No. 7, pp. 621-627, 1987 Printed in Great Britain.
0022-3697/87 $3.00+0.00 Pergamon Journals Ltd.
SUBSTRATE RESPONSE TO ATOMIC SCATTERING FROM SOLID SURFACES MADHU MENON Department of Physics, Texas A&M University, College Station, Texas 77843, U.S.A. and CHARLES W. MYLES Department of Physics and Engineering Physics, Texas Tech University, Lubbock, Texas 79409, U.S.A.
(Received 28 October 1986; accepted 16 December 1986) Abstraet--A technique for calculating the response of the substrate to atomic scattering from solid surfaces is presented. It utilizes a time-dependent phonon Green's function approach to calculate the propagation into the substrate of a disturbance induced by an atom-surface collision and it can, in principle, be incorporated into a recently introduced subspace Hamiltonian technique for calculating atomic forces from electronic energies. It therefore promises to be useful in conjunction with computer simulations of atomic motion at solid surfaces. The method is illustrated in the case of a simple model where the propagation of the disturbance from the collision can be calculated in terms of simple functions.
Keywords: Lattice dynamics, computer simulations, phonons. 1. INTRODUCTION The interaction of molecules and atoms with solid surfaces occurs in many important physical and chemical processes. Examples include catalysis, the growth of exotic materials by molecular beam epitaxy (MBE) techniques, and the formation of metal-semiconductor contacts [1-5]. The development of a technique for the realistic simulation of atomic motion at surfaces would thus be useful from both the fundamental physics and the technological viewpoints. Recently, Menon and Allen [6] have introduced a technique which is, in principle, capable of accurately treating such problems. It is based upon the Hellmann-Feynman theorem [7] and the subspace Hamiltonian technique [8]. The primary advantage of this new method is that it enables one to accurately calculate intermolecular forces in covalently bonded materials directly from the electronic structure of the material, rather than from an assumed potential. When used in conjunction with computer simulations of atomic motion at semiconductor surfaces, this technique thus promises to give a more realistic picture of such motion than the widely used simulations based upon Lennard-Jones potentials [9-11], which can provide only a qualitative description. In order to illustrate their method, Menon and Allen [6] calculated trajectories for AI and As atoms scattering off the relaxed (110) surface of GaAs. Although their method is general, these calculations were not realistic because the response of the substrate to the incoming atom was not included. In particular, the substrate atoms were assumed to be in
fixed positions and the energy lost from the incoming atom to the substrate was simulated by artificially removing a fraction of that atom's kinetic energy at each time-step. This model is thus clearly not sufficient to accurately describe the scattering of atoms from semiconductor surfaces. This situation can be partially remedied within the context of the Menon-Allen [6] method by allowing the atoms in the surface to move under the mutual influence of forces from both the incoming atom and the surrounding surface and subsurface atoms. Calculations using this approach are currently under way and will be reported elsewhere [12]. Even this generahzed approach has its practical limitations, however. In particular, the many-body nature of the problem requires one to include the forces due to and the motion of a large number of atoms in order to obtain an accurate picture of the atomic scattering from the surface. Clearly, what is needed is a method which is capable of treating the response of the substrate without the necessity of including the individual motions of huge numbers of atoms. It is the purpose of the present paper to present a formalism which will form the foundation for such a method and to illustrate its use in a simple case. Treatment of the motion of atoms at surfaces has received considerable attention in the literature. Most prominent among these studies have been calculations which have made use of the generalized Langevin equation approach [13], developed originally by Zwanzig [14], Mori [15] and Kubo [16]. These approaches incorporate Lennard-Jones type interactions to describe the dynamics of molecular motion. In this paper, we present a formalism which 621
622
MADHU MENON and CHARLESW. MYLES
will, with suitable generalization, enable one to properly include the response of the substrate atoms to atomic scattering from surfaces by allowing propagation of energy from the surface into the substrate interior. This formalism employs a time-dependent phonon Green's function method which will enable it to be incorporated into the general approach introduced by Menon and Allen [6] in a straightforward manner. Our method assumes at the outset that the motion of the substrate atoms within the immediate vicinity (say several lattice spacings) of the incoming atom position can be adequately treated within an appropriate generalization of the Menon-Allen method [12]. In such a generalized treatment, these atoms would be allowed to move under the influence of the forces caused by the disturbance of the incident atom without any imposed restriction. Our Green's function approach is most useful at the (arbitrary) "boundary" between this group of atoms and the rest of the substrate. At that point, we assume that the disturbance caused by the scattering atom travels into the substrate in a manner that can be described by the usual techniques of lattice dynamics [17]. In what follows, we shall refer to the part of the lattice where this is a valid approximation as "the substrate" and, for the purposes of this paper, ignore the part of the disturbance for which this approximation breaks down (implicitly assuming that motion in that region has been accounted for by the Menon-Allen technique) [6]. To illustrate our method, in the present paper we consider only a simple model of the substrate, where the atoms are coupled by harmonic, nearest-neighbor forces so that standard, harmonic lattice Green's functions [17] may be used to describe the propagation of the energy received by the substrate from the scattering process as it propagates inward from the surface from one lattice plane to the next. Extensions of our method to more realistic models of the solid where, for example, higherneighbor coupling and anharmonic forces are accounted for is straightforward but tedious. The advantage of first considering a simple model as we do here is that, in this case, the response of the substrate to the disturbance caused by the incident atom and the propagation of that disturbance into the substrate can be expressed in terms of elementary functions. In later papers, we plan to make some of the generalizations just mentioned and to incorporate the generalized formalism into the Menon-Allen [6] method for computer simulations. The remainder of this paper is organized as follows. In section 2, the formalism for incorporating the response of the substrate to an atomic collision is outlined and the simple model we use in our illustration of the method is discussed. In section 3, we present some results obtained on the basis of this model for the propagation of a disturbance inward from the surface. Finally, in section 4, we make a brief summary and give a brief discussion of the
generalizations in the model which are necessary for its application to real systems.
2. THEORY
2.1. Model J The model we consider here for the substrate is that of a semi-infinite, simple cubic, monatomic lattice with masses M coupled harmonically by nearest-neighbor, central forces characterized by force constants K. More realistic models, which, for example, allow for more complicated lattice structures, allow for more than one atom per unit cell, and include second- and higher-neighbor coupling, can be treated using straightforward generalizations of the model we discuss here. In what follows, we adopt the convention that the positive z axis is in the inward normal direction from the substrate surface. Furthermore, while realizing that a real disturbance which is caused by the scattering of an atom with a solid surface will emanate in all directions into the substrate from the center of scattering, in the present paper we consider in detail only such disturbances which propagate primarily in the + z direction, that is, in the direction which is inward and perpendicular to the surface. This simplified model contains much of the essential physics of the response of the substrate to the disturbance, yet it contains none of the computational and formal complexities of more realistic models. In particular, this model enables one to calculate the response of the substrate to a disturbance in terms of integrals of elementary functions. We also note that, while our model for the propagation of the disturbance is essentially onedimensional, since we consider propagation normal to the surface only, it employs the full, threedimensional phonon Green's function for the semiinfinite crystal [17]. Clearly, a more realistic model for both the underlying lattice and the propagation of the disturbance will have to be employed in conjunction with computer simulations of real atom-surface scattering. 2.2. Equation of motion and formal solution Under the assumptions just discussed, and if one considers the solid as made up of layers of atoms in planes perpendicular to the z direction, the z component of the equation of motion for an atom in layer number n, driven by a time dependent driving force Fn(t), which represents the disturbance, is d2U~
M--~-=K(Un+~+Un_~-2Un)+Fn(t),
(1)
where Un is the z-directed displacement of that atom from equilibrium [17] and t is the time, Here, we have made use of our assumption that we are considering only propagation in the + z direction by labeling an atom in a particular layer with the index of that layer
Atomic scattering from solid surfaces and suppressing the other two indices, corresponding to the x and y coordinates of the atom, needed to specify the atomic position of that atom. Our numbering convention is that the surface layer is labeled by n = 0 and that n increases inward from the surface. Thus, eqn (1) is formally valid only for positive, nonzero n. However, in what follows we shall treat it as if it were valid for all n and build in the surface boundary conditions using a Green's function technique. In the following derivation, we find it convenient to use both the notation and some of the results from the calculus of discrete variables [18], which, to our knowledge, has not often been utilized in physics, despite the discrete nature of the physical systems which one usually has to treat. The use of this notation, while not necessary to achieve our final results, emphasizes the discreteness of the problem of interest (we desire to calculate, layer by layer, the time dependence of a disturbance propagating inward from the surface) and serves as a shorthand notation for some otherwise rather cumbersome expressions. A more detailed discussion of the calculus of discrete variables and its applications to solid state physics may be found in Ref. [18]. To begin with, we define the discrete analog of the spatial derivative operator by its operation on the displacement Un as [17, 18]: DU~= U ~ + , - U~,
m2~ f ~ dt'[G(n,n',t-t')D~Un,(t ') -- Un,(t')Dl~G (n, n', t - t')] d2
+~f~_ dt'IG(n,n',t--t')-d-~Un,(t')
-Un.(t')~i-aG(n, n', t - r) 2 j_ dt'G(n,n',t-t')f,.(t')- G(t). (5) In eqn (5), the sum on n' goes from the surface layer, n" = 0, to a layer n' = N, where we assume that layer N is far removed from the surface. In our final results we shall effectively take N --* oo. On using the discrete analog of Green's theorem [18] in the first term on the left-hand side in eqn (5) to evaluate the sum on n', integrating by parts in the second term on the lefthand side to evaluate the integral over t', and rearranging the resulting expression, one obtains for the displacement at site n and time t:
U~(t)=~ f? dt'G(n,n',t-t')A,(t') +co
(2b)
This notation allows us to rewrite eqn (1) in the more compact form:
(c020Dl~ + ff-~2t2) Un(t) = f~(t ),
one obtains:
(2a)
and its adjoint operator as DU~ = - ( U ~ - U~_,).
623
(3)
we have defined fn(t)=F,(t)/M and co~= K/M. Equation (3) is obviously analogous to
dt'[G(n,N + l,t-t')DUN(t')
- G (n, l, t -- t')DUo(t') -- UN+,(t')DG (n, N, t - t') + UI(t')DG (n, O, t - t')]
- ~ G(n,n',t-t )-d-~go,(t')
where
the wave equation in a continuous system, with DD playing the role of the negative of the Laplacian [18]. The Green's function G ( n , n ' , t - t ' ) for the differential equation [eqn (3)] satisfies
(
d2)
co~Dl~ + - ~
G(n,n',t-t')=fn.,,6(t-t').
(4)
Upon changing n to n' and t to t' in eqn (3), multiplying that equation by G (n, n', t - t'), multiplying eqn (4) by U,,(t'), summing both equations over all n', integrating both equations over all t', subtracting the resulting equations, and noting that d2G d2G dt 2 = d t , 2,
- U , , ( t ' ) d - ~ a ( n , n ' , t - t ' ) ] ~ ~.
(6)
The terms in the second integral in eqn (6) play the role of the surface terms in application of Green's theorem to continuous systems. At this point, we assume that the disturbance decreases with time sufficiently that we can set the last term on the right-hand side of eqn (6) to zero. The resulting equation is rather general within the context of our model. In particular, we have not yet specified the form of the Green's function G (n, n', t - t'). We note in passing that, for example, layer N could correspond to a second surface layer in the problem, so that the present formalism could be applied to the problem of a "slab". However, in our model, layer N is far removed from the surface (i,e. it is deep in the substrate). We thus assume that as N becomes very
624
MADI-IU MENON and CHARLES W. MYLES
large (strictly N --* oo), the Green's functions G (n, N + 1, t - t ' ) and G (n, N, t - t' ) become very small (strictly --*0) for small n (i.e. for layer n near the surface layer). Finally, we note that because of the definition of the operator D [eqn (2a)] we can write:
UI(t')DG (n, O, t - t') - G (n, 1, t - t')DUo(t') = Uo(t')G(n, 1, t - t ' ) -
Ul(t')G(n,O,t -t').
(7)
With these simplifications, eqn (6) becomes
U.(t) = X ~' j _
d t ' G ( n , n ' , t - t')f.,(t')
+e~g
dt'[Uo(t')G(n, 1, t - t ' )
-- U, ( t ' ) G (n, O, t -
t')].
(8)
Equation (8) is the general expression for the displacement (in the z direction) of an atom in the nth layer, under the approximations of our model. Clearly, given the Green's function G and the disturbance f,(t), the displacement U,(t) of any layer below the surface can be calculated as a function of time. It should also be noted that the displacements for the first two layers (n = 0 and n = 1) are coupled, but that once solutions for them are obtained, the displacement for any other layer is, in principle, obtained in a straightforward manner. In this paper, we consider solutions to eqn (8) only in the case of a delta function impulse at the surface layer. While such an approximation is not necessary to implement the formalism just described, it enables us to obtain solutions which are of a particularly simple form, but which still elucidate much of the basic physics of the propagation of a disturbance into the substrate. In obtaining the results discussed in the next section, we have therefore assumed that:
f~(t) =f3n,o 3 (t),
(9)
were f is a constant. With this assumption, eqn (8) becomes
U,(t) = f G ( n , O, t) +oa~ - gl(t')G
dt'[Uo(t')G(n, l , t - t ' ) (n, O, t - - t ' ) l .
(lO)
2.3. Green's function f o r the semi-infinite crystal In order to solve eqn (10), one obviously needs the Green's function for the semi-infinite crystal. The dynamic, frequency dependent Green's function for this case (or its analog for electrons in the tightbinding approximation) has been obtained by several authors [17, 19, 20]. Most derivations of this function
take a similar approach. In particular, the surface is treated as a two-dimensional defect in the infinite crystal, where the "defect potential" is obtained by cutting bonds on the opposite sides of a hypothetical plane which bisects the crystal but which contains no atoms itself. One then solves the Dyson equation to obtain the Green's function for the semi-infinite crystal from that for the infinite crystal. For the purposes of the present paper, it is particularly convenient to use the results obtained for this function by Kalkstein and Soven [19], suitably modified for the present nearest-neighbor phonon problem (the authors of Ref. [19] treated the electronic problem in the semi-infinite crystal in the nearest-neighbor tightbinding approximation). In the present notation, the expression for this function in the mixed Bioch-Wannier representation is:
G (k, n, n', 6o2) =(i/o0~){l - [r (k, o ' ) 1 : } - ' ~ X E(r ~ , 022) + i {1 i [r ~ , (.02)]2}1/2 >ln-n'l - ( r Ot, co2) + i{1 -- [r(k, r~2)]2 }'/2)" +"' + 2],
(li) where k is a wavevector in the plane parallel to the surface (the integral over the z directed wavevector has already been evaluated), ~o2 is the square of a frequency, n and n' (as before) have the meaning of layer indices in the z direction and we have defined o2 2 r (k, co2) = co--~- 3 + cos(kxa) + cos(kya). (12) In this expression, the last two terms are related to the phonon dispersion relation parallel to the surface, kx and ky are the cartesian components of the planar wavevector k, and a is the lattice spacing. In equation (11), the square roots should be interpreted in their complex sense. That is, they are positive and real if (r (k, o2))2< 1 and they are purely imaginary if (r (k, ¢o2))2 > 1, with the overall sign of the square root in the latter case being the same as that of r 01, Co2). The time-dependent Green's function needed in this paper is obtained by Fourier transforming eqn (11) as
'M
G(k,n,n',t)=~-~n
dcoG(k,n,n',co2)e i~'. (13)
This integral cannot be evaluated in closed form for abitrary k. However, since in the present paper we are considering propagation only in the z direction, we take k = 0 in what follows. Thus, eqn (12) simply becomes r(k=O, to2)=o92/o) g - 1. We remark in passing that, in this case, the Green's function of eqn (11) reduces to that for the semi-infinite monatomic
Atomic scattering from solid surfaces chain, as obtained by Maradudin [17] and others. In what follows we suppress any further reference to k. Under these conditions, it is convenient to change variables of integration in eqn (13) by setting
co = ,~2coo cos ~b,
(14)
and integrating over all ~. Then, it is not difficult to show that, if one assumes n > n', the integral can be transformed to
G(n, n', t) = ~/20 (t) 2nt9 o
x
d~ ei'/2~'~°~e2J(~+J)~sin [(2n' 4- 2)~ ] sin ~b sin (2~b)
fo
(15)
where the symmetry between the indices n and n' has been broken by the assumption n > n'. If instead we had assumed n ' > n, the above formula would hold with n and n' interchanged. Here, O (t) is the unit step function. From eqn (10), it is clear that the Green's functions of primary interest in this paper are those that represent propagation from the surface layer (n' = 0) and the first subsurface layer ( n ' = 1) to layer n. Using eqn (15) and the integral representation of the Bessel function [21],? it is easily shown that these two functions can be written as
o (t) G(n, 0, t ) = ~ _ [J~+3(T )-J~+t(T)] (16a) VZWo and
625 3. RESULTS
3.1. Approximate solution for the time-dependent displacement Clearly, given the general form of the Green's function [eqn (1 l)] the equation of motion, eqn (10), can be solved exactly using Fourier transform techniques. Briefly stated, to carry out such a solution, one first Fourier transforms eqn (10) in time and obtains two coupled algebraic equations for the Fourier components U0(co) and Ul(co) of the displacements of an atom in the surface (n = 0) and the first subsurface (n = l) layers by setting n = 0 and n = 1 in the resulting general expression for the Fourier component of the n th layer displacement U,(m). One then solves these algebraic equations for U0(co) and Ul(CO) and substitutes those solutions back into the expression for Un(co). The displacement U,(t) of any layer n can then be obtained by taking the inverse Fourier transform of the resulting expression. This procedure is briefly outlined in the Appendix. For the purposes of the present paper, however, it is more instructive to obtain a solution using an interation procedure in conjunction with the time-dependent Green's functions in eqns (16a) and (16b). Then, one can obtain an approximate solution which involves only the integrals of the Bessel functions J~,(T). In the simplest approximation, this reduces to a solution which is a series of such functions. Substituting eqns (16a) and (16b) in eqn (10), an exact expression for the time-dependent displacement of an atom in the nth layer is:
U " ( t ) = ~d [ J zo~ +T 3 ( T ) -' J ~ + l v( T ) ] + ~ f :
o (t) G (n, 1, t) = ~ [Jz,+s(r) -- J~+ 3(T)
T"
+J2~+I(T)- J2n_I(T)], (16b) where we have defined the characteristic dimensionless time
r = ~/2~o0t,
(16~)
and J.,(t) is the Bessel function of order m.
t In order to obtain eqns (16a) and (16b) from eqn 05) we have used the integral representation of the Bessel function in the form
J,.(x)
(-1)mf~dOe~°,~°°s° It
:~In order to obtain eqns (18a), (18b) and (18c) from the first iterative solution to eqn (17), we have used the following relation:
:f
dxJv(x)J~(z - x) = 2 ~. (-- 1)'J7+~+2,+ ,(z).
P.C.S. 48/7--C
I~0
x [J:,+3(r- r') - J~+l(r-
) r')] t,
(17)
where T ' = x/2co0t'. Although it is not difficult to solve eqn (17) numerically (we do so in the next subsection), it is also useful to obtain approximate solutions by an iteractive procedure. In this procedure, the first term in eqn (17) is taken as the zeroth-order approximation to U,(t). That result (for n = 0 and 1) is then substituted into the integrand of the second term of that equation to generate the next approximation, etc. If one makes the approximation that T ~ 1, which physically corresponds to the large mass, short-time limit, this procedure can be truncated after the first interation. Then, using the properties of Bessel functions [21],:[: one can obtain an approximate form for Un(t) as a series of such
626
MADHU MENON and CHARLES W. MYLES
functions. This results in the following approximate expressions for the displacements of an ato~n in the first three layers: - f
Uo(t) = ~
[JI(T) +
J3(T) + 4Js(T ) -
4s(T)], (18a)
U, (t) = ~
[2& ( r ) + 4s (T)]
(18b)
clear from Figs l(a)-(d) that the initial peak of the disturbance becomes more and more damped out as the disturbance reaches each successive layer inward from the surface. Thus, as one intuitively expects, there is damping of the disturbance spatially as it propagates inward from the surface, as well as damping of the displacements with time within a layer. Thus, while the model we have used to obtain the results shown in Figs l(a)-(d) is quite simple, we (a)
and
U2(t)= ~2~oS(T),
(18c)
Z
where we have defined the quantity
s(T)
= ~,, ( - 1)tJ2I+5(T).
2
(18d)
4
$
8
10
T
I=l
(b) Because of the nature of the Bessel functions as the order increases [21], in practice, it is necessary to keep only the first few terms in the series defined by eqn (18d). As is discussed below, for short times, eqns (18a)-(18d) are rather accurate representations of the exact solution to eqn (17). 3.2. Numerical results We have calculated the time dependence of the displacements U,(t) for the surface layer and the first few subsurface layers both by solving the exact expression [eqn (17)] numerically and also by utilizing eqn (18). Although we find that there are some differences between the two solutions, the approximate solutions obtained from eqn (18) are qualitatively very similar to the exact solutions obtained from eqn (17). Not surprisingly, the approximate solutions are very accurate for small values of the dimensionless time parameter (T) and become less accurate as time progresses. In Figs l(a)-(d), we plot the results of our numerical solutions to eqn (17) for Un(t) for n = 0, 1, 2 and 3, respectively. These solutions for the successive layers behave qualitatively as one would expect. At t = 0, the impulse is given to the surface layer (n = 0) [Fig. l(a)] and the resulting disturbance in that layer displays a damped, oscillatory behavior with time. Figure l(b) shows that the initial peak of the disturbance reaches the first subsurface layer (n = 1) a characteristic time to= l/x/2co0 after it peaks in the surface layer. This must, of course, correspond to the time it takes the disturbance to propagate from the surface layer to the subsurface layer or to tc = a/f, where a is the lattice spacing and f is the characteristic velocity of the impulse, as indicated in eqn (9). Likewise, Figs l(c) and l(d) indicate that the initial peak of the disturbance reaches each successive layer a delay time tc after it has reached the previous layer. Furthermore, it is
2
4
6
8
10
T (c)
2
4
6
8
0
T (d)
Z
2
4
6
8
10
T Fig. 1. Displacement of surface atom (l = 0), first subsurface atom (l = 1), second subsurface atom (l = 2), and third subsurface atom (I = 3) in semi-infinitelinear chain of atoms coupled by springs as a function of the scaled time T = x/2~o0t. Here, t is the time and to0 = K x / ~ , where K is the force constant of each spring and m is the atomic mass. An impulsive forcef 6 (t - to) has been exerted on the surface atom at t = t0. Notice that the atom vibrates with decreasing amplitude as energy propagates into the interior of the chain.
Atomic scattering from solid surfaces feel that these solutions display most of the qualitative physics one would expect if such a calculation were carried out using more realistic models of the substrate and the atom-surface collision. 4. SUMMARY
We have described a technique for incorporating the response of the substrate to atomic scattering from solid surfaces. It utilizes a time-dependent phon o n Green's function approach to calculate the propagation of a collision-induced disturbance into the substrate and it can, in principle, be incorporated into the M e n o n - A l l e n [6, 8] technique for computer simulation of such processes. While we have illustrated the method only in the case of a simple model of both the substrate and of the atom-surface collision, the solutions for the substrate response in this case display most of the qualitative physics one expects for a more realistic model. Clearly, in order to use the present technique in conjunction with computer simulations of real atom-semiconductor surface scattering, more realistic models of both the substrate and the scattering process need to be used. In particular, the substrate needs to be modeled using a realistic (zincblende) lattice structure for such a material and more than nearest-neighbor force constants must be used to characterize the lattice vibrations. Furthermore, the model of the atom-surface collision as a delta function impulse which causes a disturbance travelling only perpendicular to the surface is clearly not realistic and more elaborate models which allow for propagation of energy in all directions away from the scattering center are needed to treat real systems. Despite the simplicity of the model we have used to illustrate the method presented here, we feel that the general approach that we have outlined promises to enable one to accurately treat the substrate response to atom-surface scattering and that it thus has numerous potential applications in the computer simulation of such processes.
627
3. Monch W., Surf. Sci. 132, 92 (1983) and Refs therein. 4. Williams R. H., Surf. Sci. 132, 122 (1983) and Refs therein. 5. Wieder H. H., Inst. Phys. Conf. Ser. 50, 234 (1980); Appl. Phys. Lett. 38, 170 (1980). 6. Menon M. and Allen R. E., Phys. Rev. B 33, 7099 (1986). 7. Feynman R. P., Phys. Rev. 56, 340 (1939). 8. Allen R. E. and Menon M., Phys. Rev. B 33, 5611 (1986). 9. Allen R. E., deWette F. W. and Rahman A., Phys. Rev. 179, 887 (1969). 10. Tully J. C., Surf. Sci. 125, 282 (1983). I I. Abraham F. F., J. Vac. ScL TechnoL B 2, 534 (1984) and Refs therein. 12. Menon M., Myles C. W. and Allen R. E., unpublished. 13. Tully J. C., J, chem. Phys. 73, 1975 (1980). 14. Zwanzig R., 3". Chem. Phys. 32, i173 (1960). 15. Mori H., Prog. Theor. Phys. 33, 423 (1965). 16. Kubo R., Rep. Prog. theor. Phys. 29, 255 (1966). 17. Maradudin A. A., Montroll E. W., Weiss G. H. and Ipatova I. P., "Theory of lattice dynamics in the harmonic approximation" in Solid St. Phys,, Suppl. 3 (Edited by H. Edrenreich, F. Seitz and D. Turnbull), 2nd edn. Academic Press, New York (1971). 18. Menon M. and Allen R. E., unpublished. 19. Kalkstein D. and Soven P., Surf Sci. 26, 85 (1971). 20. Allen R. E., Phys. Rev. B 19, 917 (1979). 21. Abramowitz M. and Stegun I., Handbook of Mathematical Functions. Dover, New York (1965). APPENDIX
Fourier transform solution to eqn (10) As is mentioned in the text, it is possible to solve the equation of motion [eqn (I0)] exactly using Fourier transform techniques. Here, we briefly outline this method of solution. One begins by Fourier transforming eqn (10) in time. This results in the expression U,(co) = f G (n, O, co2) + 2nco]iUo(co) G (n, 1, 0 5) -U,(co)a(n,O, co2)J, (A.1) where G(n,n',co 2) (for n ' - 0 and 1) is given by eqn (11) with k ~ 0. In order to obtain the solution to eqn (A.1), one must first write that equation for U0(co) and U,(co) and solve the resulting coupled algebraic equations for those quantities. The solutions which are obtained upon carrying out this procedure are easily shown to be Uo(co) = f G (0, O, co2)[1 + (2nco02)2D(co)]
Acknowledgements--We thank the Office of Naval Research (N00014-82-K-0447) for their support, which made this work possible. C. W. M. is grateful to the Department of Physics at Texas A&M University for their hospitality while the majority of this work was carried out. It is a pleasure to acknowledge stimulating discussions with Professor R. E. Allen.
(A.2)
and U, (co) = f [ G (I, 0, co2) +2ncogD (co)l[1 + (2nto~)2D(co)l,
(A.3)
where we have defined D (co) = [G (0,0, co2)G(l, l, co2)
REFERENCES
1. Spicer W. E., Chye P. W., Skeath P. R., Su C. Y. and Lindau I., Vac. Sci. Technol. 16, 1422 (1979). 2. Spicer W. E., Lindau I., Skeath P. R. and Su C. Y., J. Vac. Sci. Technol. 17, 1019 (1980).
-G(0, I, co2)G (I, 0, coz)].
(A.4)
Upon combining eqns (A.Iy--(A.4) and inverse Fourier transforming, the solution to the time-dependent displacement U,,(t)can be obtained for arbitrary n.
0022-3697/87 $3.00+0.00 Pergamon Journals Ltd.
SUBSTRATE RESPONSE TO ATOMIC SCATTERING FROM SOLID SURFACES MADHU MENON Department of Physics, Texas A&M University, College Station, Texas 77843, U.S.A. and CHARLES W. MYLES Department of Physics and Engineering Physics, Texas Tech University, Lubbock, Texas 79409, U.S.A.
(Received 28 October 1986; accepted 16 December 1986) Abstraet--A technique for calculating the response of the substrate to atomic scattering from solid surfaces is presented. It utilizes a time-dependent phonon Green's function approach to calculate the propagation into the substrate of a disturbance induced by an atom-surface collision and it can, in principle, be incorporated into a recently introduced subspace Hamiltonian technique for calculating atomic forces from electronic energies. It therefore promises to be useful in conjunction with computer simulations of atomic motion at solid surfaces. The method is illustrated in the case of a simple model where the propagation of the disturbance from the collision can be calculated in terms of simple functions.
Keywords: Lattice dynamics, computer simulations, phonons. 1. INTRODUCTION The interaction of molecules and atoms with solid surfaces occurs in many important physical and chemical processes. Examples include catalysis, the growth of exotic materials by molecular beam epitaxy (MBE) techniques, and the formation of metal-semiconductor contacts [1-5]. The development of a technique for the realistic simulation of atomic motion at surfaces would thus be useful from both the fundamental physics and the technological viewpoints. Recently, Menon and Allen [6] have introduced a technique which is, in principle, capable of accurately treating such problems. It is based upon the Hellmann-Feynman theorem [7] and the subspace Hamiltonian technique [8]. The primary advantage of this new method is that it enables one to accurately calculate intermolecular forces in covalently bonded materials directly from the electronic structure of the material, rather than from an assumed potential. When used in conjunction with computer simulations of atomic motion at semiconductor surfaces, this technique thus promises to give a more realistic picture of such motion than the widely used simulations based upon Lennard-Jones potentials [9-11], which can provide only a qualitative description. In order to illustrate their method, Menon and Allen [6] calculated trajectories for AI and As atoms scattering off the relaxed (110) surface of GaAs. Although their method is general, these calculations were not realistic because the response of the substrate to the incoming atom was not included. In particular, the substrate atoms were assumed to be in
fixed positions and the energy lost from the incoming atom to the substrate was simulated by artificially removing a fraction of that atom's kinetic energy at each time-step. This model is thus clearly not sufficient to accurately describe the scattering of atoms from semiconductor surfaces. This situation can be partially remedied within the context of the Menon-Allen [6] method by allowing the atoms in the surface to move under the mutual influence of forces from both the incoming atom and the surrounding surface and subsurface atoms. Calculations using this approach are currently under way and will be reported elsewhere [12]. Even this generahzed approach has its practical limitations, however. In particular, the many-body nature of the problem requires one to include the forces due to and the motion of a large number of atoms in order to obtain an accurate picture of the atomic scattering from the surface. Clearly, what is needed is a method which is capable of treating the response of the substrate without the necessity of including the individual motions of huge numbers of atoms. It is the purpose of the present paper to present a formalism which will form the foundation for such a method and to illustrate its use in a simple case. Treatment of the motion of atoms at surfaces has received considerable attention in the literature. Most prominent among these studies have been calculations which have made use of the generalized Langevin equation approach [13], developed originally by Zwanzig [14], Mori [15] and Kubo [16]. These approaches incorporate Lennard-Jones type interactions to describe the dynamics of molecular motion. In this paper, we present a formalism which 621
622
MADHU MENON and CHARLESW. MYLES
will, with suitable generalization, enable one to properly include the response of the substrate atoms to atomic scattering from surfaces by allowing propagation of energy from the surface into the substrate interior. This formalism employs a time-dependent phonon Green's function method which will enable it to be incorporated into the general approach introduced by Menon and Allen [6] in a straightforward manner. Our method assumes at the outset that the motion of the substrate atoms within the immediate vicinity (say several lattice spacings) of the incoming atom position can be adequately treated within an appropriate generalization of the Menon-Allen method [12]. In such a generalized treatment, these atoms would be allowed to move under the influence of the forces caused by the disturbance of the incident atom without any imposed restriction. Our Green's function approach is most useful at the (arbitrary) "boundary" between this group of atoms and the rest of the substrate. At that point, we assume that the disturbance caused by the scattering atom travels into the substrate in a manner that can be described by the usual techniques of lattice dynamics [17]. In what follows, we shall refer to the part of the lattice where this is a valid approximation as "the substrate" and, for the purposes of this paper, ignore the part of the disturbance for which this approximation breaks down (implicitly assuming that motion in that region has been accounted for by the Menon-Allen technique) [6]. To illustrate our method, in the present paper we consider only a simple model of the substrate, where the atoms are coupled by harmonic, nearest-neighbor forces so that standard, harmonic lattice Green's functions [17] may be used to describe the propagation of the energy received by the substrate from the scattering process as it propagates inward from the surface from one lattice plane to the next. Extensions of our method to more realistic models of the solid where, for example, higherneighbor coupling and anharmonic forces are accounted for is straightforward but tedious. The advantage of first considering a simple model as we do here is that, in this case, the response of the substrate to the disturbance caused by the incident atom and the propagation of that disturbance into the substrate can be expressed in terms of elementary functions. In later papers, we plan to make some of the generalizations just mentioned and to incorporate the generalized formalism into the Menon-Allen [6] method for computer simulations. The remainder of this paper is organized as follows. In section 2, the formalism for incorporating the response of the substrate to an atomic collision is outlined and the simple model we use in our illustration of the method is discussed. In section 3, we present some results obtained on the basis of this model for the propagation of a disturbance inward from the surface. Finally, in section 4, we make a brief summary and give a brief discussion of the
generalizations in the model which are necessary for its application to real systems.
2. THEORY
2.1. Model J The model we consider here for the substrate is that of a semi-infinite, simple cubic, monatomic lattice with masses M coupled harmonically by nearest-neighbor, central forces characterized by force constants K. More realistic models, which, for example, allow for more complicated lattice structures, allow for more than one atom per unit cell, and include second- and higher-neighbor coupling, can be treated using straightforward generalizations of the model we discuss here. In what follows, we adopt the convention that the positive z axis is in the inward normal direction from the substrate surface. Furthermore, while realizing that a real disturbance which is caused by the scattering of an atom with a solid surface will emanate in all directions into the substrate from the center of scattering, in the present paper we consider in detail only such disturbances which propagate primarily in the + z direction, that is, in the direction which is inward and perpendicular to the surface. This simplified model contains much of the essential physics of the response of the substrate to the disturbance, yet it contains none of the computational and formal complexities of more realistic models. In particular, this model enables one to calculate the response of the substrate to a disturbance in terms of integrals of elementary functions. We also note that, while our model for the propagation of the disturbance is essentially onedimensional, since we consider propagation normal to the surface only, it employs the full, threedimensional phonon Green's function for the semiinfinite crystal [17]. Clearly, a more realistic model for both the underlying lattice and the propagation of the disturbance will have to be employed in conjunction with computer simulations of real atom-surface scattering. 2.2. Equation of motion and formal solution Under the assumptions just discussed, and if one considers the solid as made up of layers of atoms in planes perpendicular to the z direction, the z component of the equation of motion for an atom in layer number n, driven by a time dependent driving force Fn(t), which represents the disturbance, is d2U~
M--~-=K(Un+~+Un_~-2Un)+Fn(t),
(1)
where Un is the z-directed displacement of that atom from equilibrium [17] and t is the time, Here, we have made use of our assumption that we are considering only propagation in the + z direction by labeling an atom in a particular layer with the index of that layer
Atomic scattering from solid surfaces and suppressing the other two indices, corresponding to the x and y coordinates of the atom, needed to specify the atomic position of that atom. Our numbering convention is that the surface layer is labeled by n = 0 and that n increases inward from the surface. Thus, eqn (1) is formally valid only for positive, nonzero n. However, in what follows we shall treat it as if it were valid for all n and build in the surface boundary conditions using a Green's function technique. In the following derivation, we find it convenient to use both the notation and some of the results from the calculus of discrete variables [18], which, to our knowledge, has not often been utilized in physics, despite the discrete nature of the physical systems which one usually has to treat. The use of this notation, while not necessary to achieve our final results, emphasizes the discreteness of the problem of interest (we desire to calculate, layer by layer, the time dependence of a disturbance propagating inward from the surface) and serves as a shorthand notation for some otherwise rather cumbersome expressions. A more detailed discussion of the calculus of discrete variables and its applications to solid state physics may be found in Ref. [18]. To begin with, we define the discrete analog of the spatial derivative operator by its operation on the displacement Un as [17, 18]: DU~= U ~ + , - U~,
m2~ f ~ dt'[G(n,n',t-t')D~Un,(t ') -- Un,(t')Dl~G (n, n', t - t')] d2
+~f~_ dt'IG(n,n',t--t')-d-~Un,(t')
-Un.(t')~i-aG(n, n', t - r) 2 j_ dt'G(n,n',t-t')f,.(t')- G(t). (5) In eqn (5), the sum on n' goes from the surface layer, n" = 0, to a layer n' = N, where we assume that layer N is far removed from the surface. In our final results we shall effectively take N --* oo. On using the discrete analog of Green's theorem [18] in the first term on the left-hand side in eqn (5) to evaluate the sum on n', integrating by parts in the second term on the lefthand side to evaluate the integral over t', and rearranging the resulting expression, one obtains for the displacement at site n and time t:
U~(t)=~ f? dt'G(n,n',t-t')A,(t') +co
(2b)
This notation allows us to rewrite eqn (1) in the more compact form:
(c020Dl~ + ff-~2t2) Un(t) = f~(t ),
one obtains:
(2a)
and its adjoint operator as DU~ = - ( U ~ - U~_,).
623
(3)
we have defined fn(t)=F,(t)/M and co~= K/M. Equation (3) is obviously analogous to
dt'[G(n,N + l,t-t')DUN(t')
- G (n, l, t -- t')DUo(t') -- UN+,(t')DG (n, N, t - t') + UI(t')DG (n, O, t - t')]
- ~ G(n,n',t-t )-d-~go,(t')
where
the wave equation in a continuous system, with DD playing the role of the negative of the Laplacian [18]. The Green's function G ( n , n ' , t - t ' ) for the differential equation [eqn (3)] satisfies
(
d2)
co~Dl~ + - ~
G(n,n',t-t')=fn.,,6(t-t').
(4)
Upon changing n to n' and t to t' in eqn (3), multiplying that equation by G (n, n', t - t'), multiplying eqn (4) by U,,(t'), summing both equations over all n', integrating both equations over all t', subtracting the resulting equations, and noting that d2G d2G dt 2 = d t , 2,
- U , , ( t ' ) d - ~ a ( n , n ' , t - t ' ) ] ~ ~.
(6)
The terms in the second integral in eqn (6) play the role of the surface terms in application of Green's theorem to continuous systems. At this point, we assume that the disturbance decreases with time sufficiently that we can set the last term on the right-hand side of eqn (6) to zero. The resulting equation is rather general within the context of our model. In particular, we have not yet specified the form of the Green's function G (n, n', t - t'). We note in passing that, for example, layer N could correspond to a second surface layer in the problem, so that the present formalism could be applied to the problem of a "slab". However, in our model, layer N is far removed from the surface (i,e. it is deep in the substrate). We thus assume that as N becomes very
624
MADI-IU MENON and CHARLES W. MYLES
large (strictly N --* oo), the Green's functions G (n, N + 1, t - t ' ) and G (n, N, t - t' ) become very small (strictly --*0) for small n (i.e. for layer n near the surface layer). Finally, we note that because of the definition of the operator D [eqn (2a)] we can write:
UI(t')DG (n, O, t - t') - G (n, 1, t - t')DUo(t') = Uo(t')G(n, 1, t - t ' ) -
Ul(t')G(n,O,t -t').
(7)
With these simplifications, eqn (6) becomes
U.(t) = X ~' j _
d t ' G ( n , n ' , t - t')f.,(t')
+e~g
dt'[Uo(t')G(n, 1, t - t ' )
-- U, ( t ' ) G (n, O, t -
t')].
(8)
Equation (8) is the general expression for the displacement (in the z direction) of an atom in the nth layer, under the approximations of our model. Clearly, given the Green's function G and the disturbance f,(t), the displacement U,(t) of any layer below the surface can be calculated as a function of time. It should also be noted that the displacements for the first two layers (n = 0 and n = 1) are coupled, but that once solutions for them are obtained, the displacement for any other layer is, in principle, obtained in a straightforward manner. In this paper, we consider solutions to eqn (8) only in the case of a delta function impulse at the surface layer. While such an approximation is not necessary to implement the formalism just described, it enables us to obtain solutions which are of a particularly simple form, but which still elucidate much of the basic physics of the propagation of a disturbance into the substrate. In obtaining the results discussed in the next section, we have therefore assumed that:
f~(t) =f3n,o 3 (t),
(9)
were f is a constant. With this assumption, eqn (8) becomes
U,(t) = f G ( n , O, t) +oa~ - gl(t')G
dt'[Uo(t')G(n, l , t - t ' ) (n, O, t - - t ' ) l .
(lO)
2.3. Green's function f o r the semi-infinite crystal In order to solve eqn (10), one obviously needs the Green's function for the semi-infinite crystal. The dynamic, frequency dependent Green's function for this case (or its analog for electrons in the tightbinding approximation) has been obtained by several authors [17, 19, 20]. Most derivations of this function
take a similar approach. In particular, the surface is treated as a two-dimensional defect in the infinite crystal, where the "defect potential" is obtained by cutting bonds on the opposite sides of a hypothetical plane which bisects the crystal but which contains no atoms itself. One then solves the Dyson equation to obtain the Green's function for the semi-infinite crystal from that for the infinite crystal. For the purposes of the present paper, it is particularly convenient to use the results obtained for this function by Kalkstein and Soven [19], suitably modified for the present nearest-neighbor phonon problem (the authors of Ref. [19] treated the electronic problem in the semi-infinite crystal in the nearest-neighbor tightbinding approximation). In the present notation, the expression for this function in the mixed Bioch-Wannier representation is:
G (k, n, n', 6o2) =(i/o0~){l - [r (k, o ' ) 1 : } - ' ~ X E(r ~ , 022) + i {1 i [r ~ , (.02)]2}1/2 >ln-n'l - ( r Ot, co2) + i{1 -- [r(k, r~2)]2 }'/2)" +"' + 2],
(li) where k is a wavevector in the plane parallel to the surface (the integral over the z directed wavevector has already been evaluated), ~o2 is the square of a frequency, n and n' (as before) have the meaning of layer indices in the z direction and we have defined o2 2 r (k, co2) = co--~- 3 + cos(kxa) + cos(kya). (12) In this expression, the last two terms are related to the phonon dispersion relation parallel to the surface, kx and ky are the cartesian components of the planar wavevector k, and a is the lattice spacing. In equation (11), the square roots should be interpreted in their complex sense. That is, they are positive and real if (r (k, o2))2< 1 and they are purely imaginary if (r (k, ¢o2))2 > 1, with the overall sign of the square root in the latter case being the same as that of r 01, Co2). The time-dependent Green's function needed in this paper is obtained by Fourier transforming eqn (11) as
'M
G(k,n,n',t)=~-~n
dcoG(k,n,n',co2)e i~'. (13)
This integral cannot be evaluated in closed form for abitrary k. However, since in the present paper we are considering propagation only in the z direction, we take k = 0 in what follows. Thus, eqn (12) simply becomes r(k=O, to2)=o92/o) g - 1. We remark in passing that, in this case, the Green's function of eqn (11) reduces to that for the semi-infinite monatomic
Atomic scattering from solid surfaces chain, as obtained by Maradudin [17] and others. In what follows we suppress any further reference to k. Under these conditions, it is convenient to change variables of integration in eqn (13) by setting
co = ,~2coo cos ~b,
(14)
and integrating over all ~. Then, it is not difficult to show that, if one assumes n > n', the integral can be transformed to
G(n, n', t) = ~/20 (t) 2nt9 o
x
d~ ei'/2~'~°~e2J(~+J)~sin [(2n' 4- 2)~ ] sin ~b sin (2~b)
fo
(15)
where the symmetry between the indices n and n' has been broken by the assumption n > n'. If instead we had assumed n ' > n, the above formula would hold with n and n' interchanged. Here, O (t) is the unit step function. From eqn (10), it is clear that the Green's functions of primary interest in this paper are those that represent propagation from the surface layer (n' = 0) and the first subsurface layer ( n ' = 1) to layer n. Using eqn (15) and the integral representation of the Bessel function [21],? it is easily shown that these two functions can be written as
o (t) G(n, 0, t ) = ~ _ [J~+3(T )-J~+t(T)] (16a) VZWo and
625 3. RESULTS
3.1. Approximate solution for the time-dependent displacement Clearly, given the general form of the Green's function [eqn (1 l)] the equation of motion, eqn (10), can be solved exactly using Fourier transform techniques. Briefly stated, to carry out such a solution, one first Fourier transforms eqn (10) in time and obtains two coupled algebraic equations for the Fourier components U0(co) and Ul(co) of the displacements of an atom in the surface (n = 0) and the first subsurface (n = l) layers by setting n = 0 and n = 1 in the resulting general expression for the Fourier component of the n th layer displacement U,(m). One then solves these algebraic equations for U0(co) and Ul(CO) and substitutes those solutions back into the expression for Un(co). The displacement U,(t) of any layer n can then be obtained by taking the inverse Fourier transform of the resulting expression. This procedure is briefly outlined in the Appendix. For the purposes of the present paper, however, it is more instructive to obtain a solution using an interation procedure in conjunction with the time-dependent Green's functions in eqns (16a) and (16b). Then, one can obtain an approximate solution which involves only the integrals of the Bessel functions J~,(T). In the simplest approximation, this reduces to a solution which is a series of such functions. Substituting eqns (16a) and (16b) in eqn (10), an exact expression for the time-dependent displacement of an atom in the nth layer is:
U " ( t ) = ~d [ J zo~ +T 3 ( T ) -' J ~ + l v( T ) ] + ~ f :
o (t) G (n, 1, t) = ~ [Jz,+s(r) -- J~+ 3(T)
T"
+J2~+I(T)- J2n_I(T)], (16b) where we have defined the characteristic dimensionless time
r = ~/2~o0t,
(16~)
and J.,(t) is the Bessel function of order m.
t In order to obtain eqns (16a) and (16b) from eqn 05) we have used the integral representation of the Bessel function in the form
J,.(x)
(-1)mf~dOe~°,~°°s° It
:~In order to obtain eqns (18a), (18b) and (18c) from the first iterative solution to eqn (17), we have used the following relation:
:f
dxJv(x)J~(z - x) = 2 ~. (-- 1)'J7+~+2,+ ,(z).
P.C.S. 48/7--C
I~0
x [J:,+3(r- r') - J~+l(r-
) r')] t,
(17)
where T ' = x/2co0t'. Although it is not difficult to solve eqn (17) numerically (we do so in the next subsection), it is also useful to obtain approximate solutions by an iteractive procedure. In this procedure, the first term in eqn (17) is taken as the zeroth-order approximation to U,(t). That result (for n = 0 and 1) is then substituted into the integrand of the second term of that equation to generate the next approximation, etc. If one makes the approximation that T ~ 1, which physically corresponds to the large mass, short-time limit, this procedure can be truncated after the first interation. Then, using the properties of Bessel functions [21],:[: one can obtain an approximate form for Un(t) as a series of such
626
MADHU MENON and CHARLES W. MYLES
functions. This results in the following approximate expressions for the displacements of an ato~n in the first three layers: - f
Uo(t) = ~
[JI(T) +
J3(T) + 4Js(T ) -
4s(T)], (18a)
U, (t) = ~
[2& ( r ) + 4s (T)]
(18b)
clear from Figs l(a)-(d) that the initial peak of the disturbance becomes more and more damped out as the disturbance reaches each successive layer inward from the surface. Thus, as one intuitively expects, there is damping of the disturbance spatially as it propagates inward from the surface, as well as damping of the displacements with time within a layer. Thus, while the model we have used to obtain the results shown in Figs l(a)-(d) is quite simple, we (a)
and
U2(t)= ~2~oS(T),
(18c)
Z
where we have defined the quantity
s(T)
= ~,, ( - 1)tJ2I+5(T).
2
(18d)
4
$
8
10
T
I=l
(b) Because of the nature of the Bessel functions as the order increases [21], in practice, it is necessary to keep only the first few terms in the series defined by eqn (18d). As is discussed below, for short times, eqns (18a)-(18d) are rather accurate representations of the exact solution to eqn (17). 3.2. Numerical results We have calculated the time dependence of the displacements U,(t) for the surface layer and the first few subsurface layers both by solving the exact expression [eqn (17)] numerically and also by utilizing eqn (18). Although we find that there are some differences between the two solutions, the approximate solutions obtained from eqn (18) are qualitatively very similar to the exact solutions obtained from eqn (17). Not surprisingly, the approximate solutions are very accurate for small values of the dimensionless time parameter (T) and become less accurate as time progresses. In Figs l(a)-(d), we plot the results of our numerical solutions to eqn (17) for Un(t) for n = 0, 1, 2 and 3, respectively. These solutions for the successive layers behave qualitatively as one would expect. At t = 0, the impulse is given to the surface layer (n = 0) [Fig. l(a)] and the resulting disturbance in that layer displays a damped, oscillatory behavior with time. Figure l(b) shows that the initial peak of the disturbance reaches the first subsurface layer (n = 1) a characteristic time to= l/x/2co0 after it peaks in the surface layer. This must, of course, correspond to the time it takes the disturbance to propagate from the surface layer to the subsurface layer or to tc = a/f, where a is the lattice spacing and f is the characteristic velocity of the impulse, as indicated in eqn (9). Likewise, Figs l(c) and l(d) indicate that the initial peak of the disturbance reaches each successive layer a delay time tc after it has reached the previous layer. Furthermore, it is
2
4
6
8
10
T (c)
2
4
6
8
0
T (d)
Z
2
4
6
8
10
T Fig. 1. Displacement of surface atom (l = 0), first subsurface atom (l = 1), second subsurface atom (l = 2), and third subsurface atom (I = 3) in semi-infinitelinear chain of atoms coupled by springs as a function of the scaled time T = x/2~o0t. Here, t is the time and to0 = K x / ~ , where K is the force constant of each spring and m is the atomic mass. An impulsive forcef 6 (t - to) has been exerted on the surface atom at t = t0. Notice that the atom vibrates with decreasing amplitude as energy propagates into the interior of the chain.
Atomic scattering from solid surfaces feel that these solutions display most of the qualitative physics one would expect if such a calculation were carried out using more realistic models of the substrate and the atom-surface collision. 4. SUMMARY
We have described a technique for incorporating the response of the substrate to atomic scattering from solid surfaces. It utilizes a time-dependent phon o n Green's function approach to calculate the propagation of a collision-induced disturbance into the substrate and it can, in principle, be incorporated into the M e n o n - A l l e n [6, 8] technique for computer simulation of such processes. While we have illustrated the method only in the case of a simple model of both the substrate and of the atom-surface collision, the solutions for the substrate response in this case display most of the qualitative physics one expects for a more realistic model. Clearly, in order to use the present technique in conjunction with computer simulations of real atom-semiconductor surface scattering, more realistic models of both the substrate and the scattering process need to be used. In particular, the substrate needs to be modeled using a realistic (zincblende) lattice structure for such a material and more than nearest-neighbor force constants must be used to characterize the lattice vibrations. Furthermore, the model of the atom-surface collision as a delta function impulse which causes a disturbance travelling only perpendicular to the surface is clearly not realistic and more elaborate models which allow for propagation of energy in all directions away from the scattering center are needed to treat real systems. Despite the simplicity of the model we have used to illustrate the method presented here, we feel that the general approach that we have outlined promises to enable one to accurately treat the substrate response to atom-surface scattering and that it thus has numerous potential applications in the computer simulation of such processes.
627
3. Monch W., Surf. Sci. 132, 92 (1983) and Refs therein. 4. Williams R. H., Surf. Sci. 132, 122 (1983) and Refs therein. 5. Wieder H. H., Inst. Phys. Conf. Ser. 50, 234 (1980); Appl. Phys. Lett. 38, 170 (1980). 6. Menon M. and Allen R. E., Phys. Rev. B 33, 7099 (1986). 7. Feynman R. P., Phys. Rev. 56, 340 (1939). 8. Allen R. E. and Menon M., Phys. Rev. B 33, 5611 (1986). 9. Allen R. E., deWette F. W. and Rahman A., Phys. Rev. 179, 887 (1969). 10. Tully J. C., Surf. Sci. 125, 282 (1983). I I. Abraham F. F., J. Vac. ScL TechnoL B 2, 534 (1984) and Refs therein. 12. Menon M., Myles C. W. and Allen R. E., unpublished. 13. Tully J. C., J, chem. Phys. 73, 1975 (1980). 14. Zwanzig R., 3". Chem. Phys. 32, i173 (1960). 15. Mori H., Prog. Theor. Phys. 33, 423 (1965). 16. Kubo R., Rep. Prog. theor. Phys. 29, 255 (1966). 17. Maradudin A. A., Montroll E. W., Weiss G. H. and Ipatova I. P., "Theory of lattice dynamics in the harmonic approximation" in Solid St. Phys,, Suppl. 3 (Edited by H. Edrenreich, F. Seitz and D. Turnbull), 2nd edn. Academic Press, New York (1971). 18. Menon M. and Allen R. E., unpublished. 19. Kalkstein D. and Soven P., Surf Sci. 26, 85 (1971). 20. Allen R. E., Phys. Rev. B 19, 917 (1979). 21. Abramowitz M. and Stegun I., Handbook of Mathematical Functions. Dover, New York (1965). APPENDIX
Fourier transform solution to eqn (10) As is mentioned in the text, it is possible to solve the equation of motion [eqn (I0)] exactly using Fourier transform techniques. Here, we briefly outline this method of solution. One begins by Fourier transforming eqn (10) in time. This results in the expression U,(co) = f G (n, O, co2) + 2nco]iUo(co) G (n, 1, 0 5) -U,(co)a(n,O, co2)J, (A.1) where G(n,n',co 2) (for n ' - 0 and 1) is given by eqn (11) with k ~ 0. In order to obtain the solution to eqn (A.1), one must first write that equation for U0(co) and U,(co) and solve the resulting coupled algebraic equations for those quantities. The solutions which are obtained upon carrying out this procedure are easily shown to be Uo(co) = f G (0, O, co2)[1 + (2nco02)2D(co)]
Acknowledgements--We thank the Office of Naval Research (N00014-82-K-0447) for their support, which made this work possible. C. W. M. is grateful to the Department of Physics at Texas A&M University for their hospitality while the majority of this work was carried out. It is a pleasure to acknowledge stimulating discussions with Professor R. E. Allen.
(A.2)
and U, (co) = f [ G (I, 0, co2) +2ncogD (co)l[1 + (2nto~)2D(co)l,
(A.3)
where we have defined D (co) = [G (0,0, co2)G(l, l, co2)
REFERENCES
1. Spicer W. E., Chye P. W., Skeath P. R., Su C. Y. and Lindau I., Vac. Sci. Technol. 16, 1422 (1979). 2. Spicer W. E., Lindau I., Skeath P. R. and Su C. Y., J. Vac. Sci. Technol. 17, 1019 (1980).
-G(0, I, co2)G (I, 0, coz)].
(A.4)
Upon combining eqns (A.Iy--(A.4) and inverse Fourier transforming, the solution to the time-dependent displacement U,,(t)can be obtained for arbitrary n.