- Email: [email protected]
Computer simulations of chemisorption
Progress in Surface Science, Vol. 25(1-4), pp. 317-331,1987 Prinled In the U.S.A.
0079-6816187 SO.OO + .50 Copyright ~ 1988 Pergamon Press pic
COMPUTER SIMULATIONS OF CHEMISORPTION MADHU MENON and ROLAND E. ALLEN Center for Theoretical Physics, Oepartment of Physics, Texas A&M UnIversity, College Station, TX 77843, U.S.A.
Abstract
Green'. function techniques and tight-binding Hamiltonians provide a very versatile approach to the electronic properties of surfaces, as demonstrated in the pioneering work oC ProCessor KouteckY. Here we show that the same general approach is useful for studies of atomic and molecular motion at surfaces; i.e., one can reliably perform molecular dynamics simula.tions of the various processes that occur when atoms, atomic: dusters, or molecules impinge on a solid surface.
Since Professor Kouteckj's early work on the electronic properties of surfaces [1-4J, there have been a very large number of studies oC surCace states, chemisorption, and related phenomena employing similar methods and ideas. It is remarleable how reliable the tight-binding model has proved to be, if used responsibly. For example, it yields the correct geometries for the (110) and (111) surfaces of III-V semiconductors [5-8J. and the correct surface atates Cor at least the (110) surfaces oC III-V and II-VI materials with the zincblende structure [9-14J. It is al.o remarleable how versatile the Green's Cunction method has proved to be, as a tool for investigating many different phenomena involving surfaces, interfaces, and chemisorption. In this paper we will demonstrate this versatility atill further by showing that the Green'. 317
H. Henon and R.E. Allen
318
function approach can be applied in a new direction - namely, it can be used for efficient and reliable studies of atomic motion at surfaces. We will first describe some new results involving
Green'. functiona that provide a mathematical foundation [15,16J, and will then discuss some of the interesting physics and chemistry that emerges from such studies [17-21J. At the heart of our technique for computer simulations is the "subspace Hamiltonian" [15J, an energy-dependent effective Hamiltonian. (For concreteness, we will consider the on... electron Hamiltonian for electrons in a solid. The formalism also applies when many-body effect. are included, with H
H
-+
potential. It also applies when H
+E -
1', where E i. the self-energy and I' the chemical
D, where D is the dynamical matrix for phonon•. ) IT
....II is the N x N Hamiltonian matrix for a "large" system (e.g., an infinite solid), the Green'.
function for this system is
-+
.... ....
....
G(E) == (E 1 - Ht 1
....
(1)
,
where 1 is the N x N identity matrix. Let some n-dimensional subspace be chosen, and the rows and columns permuted so that i,i = 1,2, ... ,n within this subspace. We partition
Gin
the well-known fashion: U
G" G 021 (p' where G", GU, G", and
(J'2
)
(2)
t
are, respectively, n x n, n
X
(N - n), (N - n)
X
n, and
(N - n) X (N - n). That is, G" is the Green's function in the chosen subspace. For simplicity of notation, we write G
= G".
(3)
We then let
E - H, •• (E) == [G(E)t 1. The subspace Hamiltonian H, •• is only n
X
....
n, whereas the original Hamiltonian H is N
(4) X
N.
The price that we pay for this reduction is that H, •• is energy dependent. In addition, as will
....
be .een below, it is frequently non-Hermitian. However, the replacement of H by H ... (E) i. very advantageous despite these complications. We will also see that the subspace Hamiltonian technique has advantages over more conventional Green's·function methods.
The eigenvalue. E, and n-dimensional eigenvectors.p, are also energy dependent: H, •• (E).p, (E)
= E,(E).p,(E)
(Sa) (5b)
Computer Simulations of Chemisorption
319
Suppose that +-+
+-+
+-+
(6)
H=Ho+V,
where
+-+
Ho is
+-+
an unperturbed Hamiltonian and V has nonzero elements only in some n
subspace:
vo
0)
0
•
X n
(7)
With the partitioning of (2), the one-electron N x N Dyson equation +-+
+-+
G = Go
+-+
+-+ +-+
+ Go . V· G,
+-+
+-+
Go(E) == (E 1-
+-+
Ho)-'
(8)
(9)
is essentially equivalent to the n x n Dyson equation G= Go +GoVG,
(10)
.... is quite well known. Since (10) is equivalent to
G-' =
Go' - V,
(11)
(4) i. equivalent to E - H ••• (E) =
Go' - V.
(12)
In Ref. 15, we proved three theorems that make the subspace Hamiltonian a particularly useful and efficient tool for calculating surface and defect states, and for performing simulations of atomic motion: Theorem 1. The eigenvalues E, (E) of the subspace Hamiltonian decrease monotonically with the energy E,
ilE,(E) < 0 ilE - ,
(13)
whenever the energy E lies within a band gap - i.e., outside the continua of unbound stat...
In surface-state calculations, or other problems where there is two-dimensional translational invariance, we define G(ICE), E, (ICE), etc., where li = (k" k,) is the planar wave vector. Then (13) is valid within the band gaps at fixed li, which cover a larger range of energies than the intersection of the band gaps for allIC. Furthermore, (13) is a good approximation in energy regions where the bulk density of .tat.. is low.
M. Menon and R. E. Allen
320
The significance of (13) is that it limits the behavior of the eignevalues E.(E) of the subspace Hamiltonian. They cannot vary wildly and erratically within the band gaps, because (13) guarantees they will decrease monotonically. One consequently expects, and finds in numerical calculations, that E; (E) is a smooth function of E within the band gaps. Theorem Z. The subspace Hamiltonian reduces to the Hamiltonian of the decoupled subspace at large energies:
H••• (E)
-+
H as
lEI -+ 00,
(14)
where
(15) and
(16) just as in (2) and (3) . This theorem provides an even more important limitation on the behavior of H"", its
eigenvalues Ei(E), and its eigenvectors ,piCE), since it implies that they are all eonstant with resp'" to the energv E at large E. This allows a very great simplication in applications of the
aubspace Hamiltonian technique, as will be seen below. Theorem 3. Bound states and resonances correspond to
E,(E) = E.
(17)
This theorem appear. to provide the most effident method available for calculations of surface, interface, and deCect states and resonances. Two other useful theorems have also been proved, respectively, by ourselves 121 J and by Wang [22J:
Theor"m 4. The singularities corresponding to Ei(E) = E, E complex, are not found on the physical Riemann sheet: 1m E;(E) § 0 for 1m E ~ O.
(18)
The theorem is very useful because it allows the deformation of contours into the complex energy plane required below. Theorem 5. If the subspace region is expanded, the larger set of eigenvalues E:
Computer SImulatIons of ChemIsorptIon
321
eigenvalues E; (E) corresponding to the original subsp""e Hamiltonian II••• (E). For example, if we include 4 planes near the surf""e rather than 2, we obtain a new Bet of eigenvalues in addition to the old set that was obtained for 2 planes; but the old set is still present and unmodified. This f""t tells us that a calculation can be converged with only a relatively amall subsp""e. We also mention that the subspace Hamiltonian provides a spectral representation for the Green's function,
i i:
G(E) = ' " ojJ; (E) (E) , ~ E-E;(E)
•
(19)
and that for resonably narrow resonances, a convenient approximation to (17) is Re E;(E) = E.
(20)
This criterion has yielded surface-state resonances that are in agreement with experiment for
the (110) surfaces of ZnSe, InP, GaSb, and GaAs [91. In SijNiSi, interface-state calculations
[231, it has also been found to agree with the alternative theoretical criterion that a resonance corresponds to a peak in the local density of states,
1
PI.,(kE) = - - 1m Tr G(liE). ".
(21)
Now let us turn to the problem of molecular dynamics computer simulations. The force F. associated with some atomic coordinate x is given by F
•
=_au ax
(22)
where U is the total energy of the system. If F. can be calculated, numerical solution of the equation of motion
tPx
m-;w =
F.
(23)
is straightforward.
Let the electrons be divided into valence electrons, whose states will be regarded as adiabatically changing in response to the motion of the nuclei, and core electrons, whose atates will be regarded as rigidly following the nuclei. There are three contributions to the total energy U:
(24)
M. Menon and R.E. Allen
322
Since the present technique employs Green's functions, one could include many-body effects in the treatment of U. In this paper, however, we will assume a one-electron picture
throughout. Then U" is the sum of the energies
E.
of the occupied one-eledron states: (25)
Aiso, U.. is the energy associated with the electron-electron Coulomb repulsion, which is doubly counted in (25), and U; ••• is the ion-ion coulomb repulsion, with an ion defined to be the atomic nucleus plus core electrons.
The interaction between two atoms can be regarded as consisting ol (1) an ..ttractive bonding interaction involving the one-electron energies
(.II
and (2) a repulsive interaction rep-
resenting a combination of effects, including the Pauli exclusion principle and the direct ion-ion
repulsion. We will therelore model (24) by
(26) where U". is given by a repulsive pair potential ,,6(r) ,
(27)
U". = L,,6(r;;), i>l
and r,; is the separation of atoms i and j. Here, we calculate the electronic energies
E.
using .. semiempirical tight-binding Hamil-
tonian [24J, with the intra-atomic matrix elements
Hii
regarded as atomic energies, and the
interatomic elements Hii representing overlap of the electronic wavefunctions. The Hi; are
taken to be constant, and the H;; are taken to decrease exponentially with the interatomic separation r. To be precise, let H;;(a,{3) be the Hamiltonian matrix element coupling orbital a on atom i to orbital {3 on atom j. This element is given by the direction cosines l = I:l:z;/r, etc., and by primitive tight-binding matrix elements V.... etc. (See p. 481 of Ref. 24). For example, the matrix element H,;(p. ,p.l coupling a P. orbital to a P. orbital is given by
V••• (r) and V••• (r). We t ..ke each primitive element VCr) to decrease exponentially with r: VCr) = V(ro) exp[-a(r - roll. For the parameters employed in the present work - V, •• , Y'PIr I
a
= 2/ro.
(28) VPPIlI' I
and V'Ptr - we choose
Computer Simulations of Chemisorption
323
We fit each V(r.) (V... , etc.) to its value in Ref. 24 at some fixed distance r. which represents a typical bond length. To be specific, we take r. to be the bond length of the solid (2.45
A for GaAs and 2.54 Afor InP) . Since the present work involves only 0
for which the Bcaling of Ref. 24 is VH(r) and [dVa /d.j.,
= (-2/r'}VH (ra).
= V( •• }( •• /.)-·,
Our choice a
= 2/••
and p electrons,
we have [dV/d.j.,
=
-aV(.,)
thus matches the present logarithmic
derivatives at • = " to those of Ref. 24. The repulsive potential is also assumed to decrease exponentially with.,
4>(r) = 4>. exp[-P(. - d)j,
(29)
but in this case we take d to equal the sum of the covalent radii of the two interacting ato"",, i.e., we take the range of the repulsive interaction to scale with the size of the ..tom. U the attractive (electronic) interaction is to dominate at large distances, and the repulsive interaction at small distances, we simply choose
P=
p should be larger than a. One could try to optimize p, but
2a.
This leaves the parameter 4>., which is determined by the requirement that the tot.l energy for the bulk solid be a minimum at • = '., where '. is the experimental bond length. Specifically, we fit 4>(r) to the repulsive potential of Harrison and co-workers [25j, 4>H (.)
C/r', at. = '" Notice that [d4>/drj., choice {J
= -M(',) and
=
[d4>H /drj •• = (-4/ •• )4>H (•• ), 80 our
= 4/•• matches the present logarithmic derivative of 4> at • = "
to that of Ref. 25.
Finally, since nearest-neighbor interactions are more chemically meaningful than distant interactions, we truncate both VIr) and
4>(.); specific..lly, we t ..ke VIr)
= 0 and 4>(.) = 0 for
• > R. = 1.25•• . One could 8mooth out the delta-function singularity in the radial force ..t R., but this would m..ke only a small difference in the present simulations. One could also, of course, choose a less conservative value for Re .
Notice that only three puameters distinguish a chemic ..1 species - the , and p atomic energies, f, and fp - which determine the atom's electronegativity - and the covalent radius '. _ which determines its size. We can represent
'0
and 'p approximately by an average or
"hybrid" energy
l
= (f' + 3'p )/4.
(30)
All electronic states below the Fermi energy ., of the solid are assumed to be occupied, and all states above ., to be empty. Our use of a cutoff R, makes this assumption physically reasonable: If an atom impinging on the surface is outside the range of interaction, there are
H. Henan and R.E. Allen
324
no forces and its occupancy is irrelevant. As Boon as it is within the interaction range, it is
interacting rather strongly, forming bonding and antibonding states with the solid, and its occupancy should be approximately correct. In Refs. 15 and 21, it is shown that the electronic force F:t, resulting from
Vorl
in (22)
and (26), can be simplified to F"•
F. '" ~Im {'" '" ~
= '" F· L..J'
(31)
;
aE;(E~) 10 az
[E! -
E;(E~)]}
(32)
gE!-E;(E~)
where the E~ are some set of !Sample energies in the complex-energy plane, along the contour C'J of Fig. 1 extending from (,.
+ ioo
down to the Fermi energy (,.. Also,
E!
and
E!
are
the lower and upper points of the interval along C, represented by E~. We have found that 5 such sample energies provide an accuracy of about 1% in the total energy and a few percent
in the electronic contribution to the force. The derivatives aE,jaz are calculated from the lIellman-Feynman theorem 1151
(33) Further details are given in Refs. 15 and 21.
ReE
Fig. 1. Complex energy plane. Integration over the rapidly varying structure between
f ..
(an
energy below the relevant electronic energy bands) a.nd t.,. (the Fermi energy of the solid) would require a very large number of sample energies for an accurate calculation. We
+ if mu and tr. Since the eigenvalues of the subspace Hamiltonian are 610wly varying functions of E for E complex, replace this interval by a. finite segment of C'J between t.,.
only a. few sample energies are required.
computer SImulatIons of ChemIsorptIon
325
Now let us turn to the results for various chemical species impinging on Ihe (UO) surface of GaAs. Figure 2 shows a top view of this surface, with the choice of coordinates indicated. Our unit of length is half the GaAs bond length, or 1.225
A. Then the As atoms at the lower
and upper right, respectively, have coordinates (1.63, -1.15, 0) and (1.63, 3.45, 0) before they are allowed to relax. Figure 2 also schematically shows a top view of 6 possible chemisorption sites on the surface, with bonding to a 'ingle As or Go. surface atom (lor 4), to a pair of As or Go. atoms (2 or 5), or in a Ga-As bridge .ite (3 or 6). These positions, which are observed as initial chemisorption sites in our simula.tions of time-dependent chemisorption, have also been observed as energy minima along the surface in calculations of energy surfaces associated with
static chemisorption [26,27J. . ,.) - - - - - (
GaAs(110)
Ga
0° ~
0
As,
0
CD
AS 2
Fig. 2. Top view or (110) surrace or GaAs, showing positions of Ga and As surface atoms and approxima.te positions of 6 initial chemisorption sites which have been observed in the present simulations.
In this paper, the unit oftime is the same as the time step,
t· =
t/~t.
~t
= 1.04 x 10-" sec; i.e.,
Also, in the simulations reported below, only two atoms are allowed to move-
the surface atom at the origin and the incoming atom. The other atoms in the semi-infinite solid affect the simulation, because they are electronically coupled to the moving atoms, and they are also represented by repulsive potentials. They are, however, motionless in these first simulations with the present technique. We have also performed some simulations with 4
atoms a.lIowed to move, and the differences are about what one expects. Namely, the results in most cases are fundamentally the same as those of the present paper, but the response of additional surface atoms does have some effect on the motion of the chemisorbing atom. One could also treat the propagation of vibrational energy into the solid with the time-
M. Menon and R. E. Allen
326
dependent phonon Green's function 117,28J. In the present work, however, we simply reduce each velocity by 0.5% at each time step. This estimate of the appropriate rate of energy dissipation is based on our study of vibrational relaxation in a semi-infinite one-dimensional chain of atoms, with a Budden impulse applied at the end 117,28J. In each case, the surface atoms are initially motionless at their relaxed positions. The initial position and velocity of
the incoming atom can be read off the graph.. Recall that the unit of length is r. /2, where r.
= 2.45A is the bulk GaA. bond length.
00 x -o . ~
-' 0
d--,,"oo"'-',"o'::o-::,o:::--:,oo o ::::--:,C: o"o '-',:o::o--::,oo :::--::!,oo
- 1.50
"
Fig. 3. Simulation of Cl chemisorbing on Go.Aa (110) . The coordinate system io defined in Fig. 2; also, t'
= t/ At, where At = 1.04 X 10-" sec.
Figure 3 shows a simulation for a Cl atom released with zero velocity one GaA. bond length above the surface, at the Ga-As bridge site 3 of Fig. 2. This atom moves away from the surface As, and bonds to a single Ga at site 4 of Fig. 2. Since CI is electronegative, this is the kind of behavior one expects, Notice that the z vibrations of the CI have large amplitude and a long period, because there is only a small, a.ngle-bending restoring force for motion in
that direction. We thus predict a low-frequency vibrational mode for this site. (Bonding to more distant neighbors will increase the frequency of this mode somewhat, but it should still be relatively low.) Notice that the present technique automatically gives angle-bending as well as bond-stretching forces, with both resulting from change. in the total electronic energy
Computer Slmulatlons of Chemlsorptlon
327
as the atomic positions are varied. Also notice that one can study the equilibrium vibra.tions
of surface and adsorbed atoms, even when such vibrations are highly anharmonic. Figure 4 shows & simula.tion for an 0 atom released with exactly the same initial conditions
as the CI above. Since 0 is light, it vibrates with high frequency in the II and
%
directions.
Again, however, it vibrates with a long period and large amplitude in the:z: direction, for which there is only a small angle-bending restoring force (in a nearest-neighbor approximation). 2.
0
2. N
,.
,.
•
0
-, - , ~o1--;:'OO;;;--;';;;OO;-;J;;;OO;-;"';;;;-"";;;OO~';;;O~O--;;'0-0-"""'" t'
....
Fig. 4. 0 on GaM (110). Figure 5 show a simulation for As, which bonds at the Ga-Ga bridge site 5. The initial conditions in the present simulations - zero velocity and the Ga-As bridge aite 3 - were chosen in order not to bias the motion of the chemisorbing atom; i.e .• it can choose to
bond either to the cation (Ga) or the anion (As). As one expects, the electronegative Group
V, VI, and VII atoms tend to bond to the cation. Furthermore, it will be seen below that the electropositive Group I, II, and III atoms tend to bond to the anion.
H. Menon and R.E . Allen
328 2.
....
1.
N
1. 1.
o.
,.,
o. o. -0.
-1. 1.
2.
"
1.
O.
a
l Oa
200
300
400
500
600
t'
700 8 00 GaAs
Fig. 5. A.!. For Group IV elements, one expects that Ga and As will be about equally attractive. This is confirmed by our simulations . In Fig. 6, Si is found to bond at the As-As bridge site 2. 51
2. N
1. 1.~=======:::=:::=====i
-o.s,
,.,
-1.
-1. -2.~============~ LOr o.
"
O. - 0. 0
Fig. 6. Si.
100
200
300
400
t·
500
600
700 Boe GaAa
Computer SImulatIons of ChemIsorptIon
329
In Figure 7, an AI atom is released 1.5 GaAs bond lengths above the surface, with an initial velocity downward. In this run alone, an As (rather than Ga) was placed at the origin. The AI bonds at site 4 with As at the origin, which is equivalent to site 1 with Ga at the origin. AI
·].i---:__
---J
o
100
200
300
400
500
600
t'
700
800
GoA. (Aol
Fig. 7. AI, with As at the origin.
2." . . . . - - - - - - - - - - - - ,
z.
o."~~========~
-o.5t;: -1 .
-0 ~o
Fig. 8. Zn.
100
200
100
400
t·
)00
600
100
800
CaA,
H. Henon and R.E. Allen
330
Figure 8 shows a simulation for Zn. As one expects, this electropositive metal atom is attracted by the anion, bonding at the As-As bridge .ite 2. Finally, Fig. 9 shows a simulation for Cu. Like Zn, this electropositive metal "tom bonds at the As-As bridge site 2.
2.
c.
2. N
,. 1. -0. -1.
>-
-, . -2.
o. x
o.
-O · 5!0--;-'07.0:-:: 20'-0-:'C:O-:: O -:.7. 00:-::'O'-O-:6C:O-::O~77. 00:-.::!OO t'
G.A.
Fig. 9. Cu. These results for various species chemisorbing on the (110) surface of GaA. are of conaiderable interest. Even more interesting, however, is the simple conclusion that the present method works, so that it is now possible to simulate atomic motion in real materials and at their surfa.ces.
Acknowledgment This work WM supported by the Office of N"val Research (NOOO14-82-K-0447) and by the Robert A. Welch Foundation.
1. 2. 3. 4.
J. J. J. J.
References Kouteay Phys. Rev. 108 (1957) 13. Koutecky, Trans. Faraday Soc. 54 (1958) 1038. Koutecky, Advan. Chern. Phys. !! (1965) 85. Koutecky, Progress in Surface "nd Membrane Science 11 (1976) 1.
Computer Simulations of Chemisorption
331
5. D. J. Chadi, Phy• . Rev. Lett. 41 (1978) 1062; Phy•. Rev. Lett. 52 (1984) 1911. 6. S. Y. Tong, A. R. Lubinsky, B. J. Mr.tik and M. A. van Hove, Phy•. Rev. Ill! (1978) 3303. 7. A. Kahn, E. So, P. Mark and C. B. Duke, J. Vac. Sci. Tech. 15 (1978) 580. 8. S. Y. Tong, G. Xu and W. N. Mei, Phys. Rev. Lett. 52 (1984) 1693. 9. R. P. Beres, R. E.Allen and J. D. Dow, Solid State Commun. 45 (1983) 13; Phys. Rev. ~ (1982) 769; Phys. Rev. B26 (1982) 5702. 10. G. P. Williams, R. J. Smith and G. J . Lapeyre, J. Vac. Sci. Tech. 15 (1978) 1249. 11. A. 11uijser, J. van Laar and T. L. van Rooy, Phys. Lett. 65A (1978) 337. 12. A. Ebina, T . Unno, Y. Suda, 11. Koinuma and T . Takahashi, J. Vac. Sci. Tech. 19 (1981) 301. 13. G. P. Srivastava, I. Singh, V. Montgomery and R . H. Williams, J. Phys. .Qlll (1983) 3627. 14. H. Carstensen, R. Manzke, I. Schafer and M. Skibowski, in Pree. Eighteenth Int. ConI. on the PhV.i.. 0/ Semicondudorl, Vol. 1edited by O. Engstrom (World Scientific, Singapore, 1987). 15. R. E. Allen and M. Menon, Phys. Rev. B33 (1986) 5611. 16. M. Menon and R. E. Allen, Phy•. Rev. B33 (1986) 7099. 17. M. Menon and R. E. Allen, Super lattices and Microstructures ~ (1987) 295. 18. M. Menon &nd R. E. Allen, Solid State Commun. M (1987) 353. 19. M. Menon and R. E. Allen, Proceedings 0/ the Second International Con,.,.nce on the Structure a' Sur'ace. (in press) . 20. M. Menon and R. E. Allen, Proceedings a' the 94th Annual Svmpo.ium the American
0'
21. 22. 23. 24. 25. 26. 27. 28.
Vacuum Societv, J. Vat. Sci. Tech. (in pres.) . M. Menon and R. E. Allen, to be published. R. Wang, J. Math. Phy•. 28 (1987) 2325,2329 and 2333. H. Lim and R. E. Allen, J. Vat. Sci. Tech. ~ (1985) 1221; 1M (1985) 2328. W. A. Harrison, Electroni. Strudu,. and the Properti.. 0/ Solids (W. H. Freeman, San Francisco, 1980). E. A. Kraut and W. A. Harrison, J. Vac. Sci. Tech. ID (1984) 409; R. Enderlein and W. A. Harrison, Phy •. Rev. B30 (1984) 1967. D. J. Chadi and R. Z. Bachrach, J . Vac. Sci. Tech. 16 (1979) 1159. J. lhm and J. D. Joannopoulos, Phy•. Rev. B26 (1982) 4429 &nd references therein. M. Menon &nd C. W. Myles, J. Phy•. Chem. Solids 48 (1987) 621.
0079-6816187 SO.OO + .50 Copyright ~ 1988 Pergamon Press pic
COMPUTER SIMULATIONS OF CHEMISORPTION MADHU MENON and ROLAND E. ALLEN Center for Theoretical Physics, Oepartment of Physics, Texas A&M UnIversity, College Station, TX 77843, U.S.A.
Abstract
Green'. function techniques and tight-binding Hamiltonians provide a very versatile approach to the electronic properties of surfaces, as demonstrated in the pioneering work oC ProCessor KouteckY. Here we show that the same general approach is useful for studies of atomic and molecular motion at surfaces; i.e., one can reliably perform molecular dynamics simula.tions of the various processes that occur when atoms, atomic: dusters, or molecules impinge on a solid surface.
Since Professor Kouteckj's early work on the electronic properties of surfaces [1-4J, there have been a very large number of studies oC surCace states, chemisorption, and related phenomena employing similar methods and ideas. It is remarleable how reliable the tight-binding model has proved to be, if used responsibly. For example, it yields the correct geometries for the (110) and (111) surfaces of III-V semiconductors [5-8J. and the correct surface atates Cor at least the (110) surfaces oC III-V and II-VI materials with the zincblende structure [9-14J. It is al.o remarleable how versatile the Green's Cunction method has proved to be, as a tool for investigating many different phenomena involving surfaces, interfaces, and chemisorption. In this paper we will demonstrate this versatility atill further by showing that the Green'. 317
H. Henon and R.E. Allen
318
function approach can be applied in a new direction - namely, it can be used for efficient and reliable studies of atomic motion at surfaces. We will first describe some new results involving
Green'. functiona that provide a mathematical foundation [15,16J, and will then discuss some of the interesting physics and chemistry that emerges from such studies [17-21J. At the heart of our technique for computer simulations is the "subspace Hamiltonian" [15J, an energy-dependent effective Hamiltonian. (For concreteness, we will consider the on... electron Hamiltonian for electrons in a solid. The formalism also applies when many-body effect. are included, with H
H
-+
potential. It also applies when H
+E -
1', where E i. the self-energy and I' the chemical
D, where D is the dynamical matrix for phonon•. ) IT
....II is the N x N Hamiltonian matrix for a "large" system (e.g., an infinite solid), the Green'.
function for this system is
-+
.... ....
....
G(E) == (E 1 - Ht 1
....
(1)
,
where 1 is the N x N identity matrix. Let some n-dimensional subspace be chosen, and the rows and columns permuted so that i,i = 1,2, ... ,n within this subspace. We partition
Gin
the well-known fashion: U
G" G 021 (p' where G", GU, G", and
(J'2
)
(2)
t
are, respectively, n x n, n
X
(N - n), (N - n)
X
n, and
(N - n) X (N - n). That is, G" is the Green's function in the chosen subspace. For simplicity of notation, we write G
= G".
(3)
We then let
E - H, •• (E) == [G(E)t 1. The subspace Hamiltonian H, •• is only n
X
....
n, whereas the original Hamiltonian H is N
(4) X
N.
The price that we pay for this reduction is that H, •• is energy dependent. In addition, as will
....
be .een below, it is frequently non-Hermitian. However, the replacement of H by H ... (E) i. very advantageous despite these complications. We will also see that the subspace Hamiltonian technique has advantages over more conventional Green's·function methods.
The eigenvalue. E, and n-dimensional eigenvectors.p, are also energy dependent: H, •• (E).p, (E)
= E,(E).p,(E)
(Sa) (5b)
Computer Simulations of Chemisorption
319
Suppose that +-+
+-+
+-+
(6)
H=Ho+V,
where
+-+
Ho is
+-+
an unperturbed Hamiltonian and V has nonzero elements only in some n
subspace:
vo
0)
0
•
X n
(7)
With the partitioning of (2), the one-electron N x N Dyson equation +-+
+-+
G = Go
+-+
+-+ +-+
+ Go . V· G,
+-+
+-+
Go(E) == (E 1-
+-+
Ho)-'
(8)
(9)
is essentially equivalent to the n x n Dyson equation G= Go +GoVG,
(10)
.... is quite well known. Since (10) is equivalent to
G-' =
Go' - V,
(11)
(4) i. equivalent to E - H ••• (E) =
Go' - V.
(12)
In Ref. 15, we proved three theorems that make the subspace Hamiltonian a particularly useful and efficient tool for calculating surface and defect states, and for performing simulations of atomic motion: Theorem 1. The eigenvalues E, (E) of the subspace Hamiltonian decrease monotonically with the energy E,
ilE,(E) < 0 ilE - ,
(13)
whenever the energy E lies within a band gap - i.e., outside the continua of unbound stat...
In surface-state calculations, or other problems where there is two-dimensional translational invariance, we define G(ICE), E, (ICE), etc., where li = (k" k,) is the planar wave vector. Then (13) is valid within the band gaps at fixed li, which cover a larger range of energies than the intersection of the band gaps for allIC. Furthermore, (13) is a good approximation in energy regions where the bulk density of .tat.. is low.
M. Menon and R. E. Allen
320
The significance of (13) is that it limits the behavior of the eignevalues E.(E) of the subspace Hamiltonian. They cannot vary wildly and erratically within the band gaps, because (13) guarantees they will decrease monotonically. One consequently expects, and finds in numerical calculations, that E; (E) is a smooth function of E within the band gaps. Theorem Z. The subspace Hamiltonian reduces to the Hamiltonian of the decoupled subspace at large energies:
H••• (E)
-+
H as
lEI -+ 00,
(14)
where
(15) and
(16) just as in (2) and (3) . This theorem provides an even more important limitation on the behavior of H"", its
eigenvalues Ei(E), and its eigenvectors ,piCE), since it implies that they are all eonstant with resp'" to the energv E at large E. This allows a very great simplication in applications of the
aubspace Hamiltonian technique, as will be seen below. Theorem 3. Bound states and resonances correspond to
E,(E) = E.
(17)
This theorem appear. to provide the most effident method available for calculations of surface, interface, and deCect states and resonances. Two other useful theorems have also been proved, respectively, by ourselves 121 J and by Wang [22J:
Theor"m 4. The singularities corresponding to Ei(E) = E, E complex, are not found on the physical Riemann sheet: 1m E;(E) § 0 for 1m E ~ O.
(18)
The theorem is very useful because it allows the deformation of contours into the complex energy plane required below. Theorem 5. If the subspace region is expanded, the larger set of eigenvalues E:
Computer SImulatIons of ChemIsorptIon
321
eigenvalues E; (E) corresponding to the original subsp""e Hamiltonian II••• (E). For example, if we include 4 planes near the surf""e rather than 2, we obtain a new Bet of eigenvalues in addition to the old set that was obtained for 2 planes; but the old set is still present and unmodified. This f""t tells us that a calculation can be converged with only a relatively amall subsp""e. We also mention that the subspace Hamiltonian provides a spectral representation for the Green's function,
i i:
G(E) = ' " ojJ; (E) (E) , ~ E-E;(E)
•
(19)
and that for resonably narrow resonances, a convenient approximation to (17) is Re E;(E) = E.
(20)
This criterion has yielded surface-state resonances that are in agreement with experiment for
the (110) surfaces of ZnSe, InP, GaSb, and GaAs [91. In SijNiSi, interface-state calculations
[231, it has also been found to agree with the alternative theoretical criterion that a resonance corresponds to a peak in the local density of states,
1
PI.,(kE) = - - 1m Tr G(liE). ".
(21)
Now let us turn to the problem of molecular dynamics computer simulations. The force F. associated with some atomic coordinate x is given by F
•
=_au ax
(22)
where U is the total energy of the system. If F. can be calculated, numerical solution of the equation of motion
tPx
m-;w =
F.
(23)
is straightforward.
Let the electrons be divided into valence electrons, whose states will be regarded as adiabatically changing in response to the motion of the nuclei, and core electrons, whose atates will be regarded as rigidly following the nuclei. There are three contributions to the total energy U:
(24)
M. Menon and R.E. Allen
322
Since the present technique employs Green's functions, one could include many-body effects in the treatment of U. In this paper, however, we will assume a one-electron picture
throughout. Then U" is the sum of the energies
E.
of the occupied one-eledron states: (25)
Aiso, U.. is the energy associated with the electron-electron Coulomb repulsion, which is doubly counted in (25), and U; ••• is the ion-ion coulomb repulsion, with an ion defined to be the atomic nucleus plus core electrons.
The interaction between two atoms can be regarded as consisting ol (1) an ..ttractive bonding interaction involving the one-electron energies
(.II
and (2) a repulsive interaction rep-
resenting a combination of effects, including the Pauli exclusion principle and the direct ion-ion
repulsion. We will therelore model (24) by
(26) where U". is given by a repulsive pair potential ,,6(r) ,
(27)
U". = L,,6(r;;), i>l
and r,; is the separation of atoms i and j. Here, we calculate the electronic energies
E.
using .. semiempirical tight-binding Hamil-
tonian [24J, with the intra-atomic matrix elements
Hii
regarded as atomic energies, and the
interatomic elements Hii representing overlap of the electronic wavefunctions. The Hi; are
taken to be constant, and the H;; are taken to decrease exponentially with the interatomic separation r. To be precise, let H;;(a,{3) be the Hamiltonian matrix element coupling orbital a on atom i to orbital {3 on atom j. This element is given by the direction cosines l = I:l:z;/r, etc., and by primitive tight-binding matrix elements V.... etc. (See p. 481 of Ref. 24). For example, the matrix element H,;(p. ,p.l coupling a P. orbital to a P. orbital is given by
V••• (r) and V••• (r). We t ..ke each primitive element VCr) to decrease exponentially with r: VCr) = V(ro) exp[-a(r - roll. For the parameters employed in the present work - V, •• , Y'PIr I
a
= 2/ro.
(28) VPPIlI' I
and V'Ptr - we choose
Computer Simulations of Chemisorption
323
We fit each V(r.) (V... , etc.) to its value in Ref. 24 at some fixed distance r. which represents a typical bond length. To be specific, we take r. to be the bond length of the solid (2.45
A for GaAs and 2.54 Afor InP) . Since the present work involves only 0
for which the Bcaling of Ref. 24 is VH(r) and [dVa /d.j.,
= (-2/r'}VH (ra).
= V( •• }( •• /.)-·,
Our choice a
= 2/••
and p electrons,
we have [dV/d.j.,
=
-aV(.,)
thus matches the present logarithmic
derivatives at • = " to those of Ref. 24. The repulsive potential is also assumed to decrease exponentially with.,
4>(r) = 4>. exp[-P(. - d)j,
(29)
but in this case we take d to equal the sum of the covalent radii of the two interacting ato"",, i.e., we take the range of the repulsive interaction to scale with the size of the ..tom. U the attractive (electronic) interaction is to dominate at large distances, and the repulsive interaction at small distances, we simply choose
P=
p should be larger than a. One could try to optimize p, but
2a.
This leaves the parameter 4>., which is determined by the requirement that the tot.l energy for the bulk solid be a minimum at • = '., where '. is the experimental bond length. Specifically, we fit 4>(r) to the repulsive potential of Harrison and co-workers [25j, 4>H (.)
C/r', at. = '" Notice that [d4>/drj., choice {J
= -M(',) and
=
[d4>H /drj •• = (-4/ •• )4>H (•• ), 80 our
= 4/•• matches the present logarithmic derivative of 4> at • = "
to that of Ref. 25.
Finally, since nearest-neighbor interactions are more chemically meaningful than distant interactions, we truncate both VIr) and
4>(.); specific..lly, we t ..ke VIr)
= 0 and 4>(.) = 0 for
• > R. = 1.25•• . One could 8mooth out the delta-function singularity in the radial force ..t R., but this would m..ke only a small difference in the present simulations. One could also, of course, choose a less conservative value for Re .
Notice that only three puameters distinguish a chemic ..1 species - the , and p atomic energies, f, and fp - which determine the atom's electronegativity - and the covalent radius '. _ which determines its size. We can represent
'0
and 'p approximately by an average or
"hybrid" energy
l
= (f' + 3'p )/4.
(30)
All electronic states below the Fermi energy ., of the solid are assumed to be occupied, and all states above ., to be empty. Our use of a cutoff R, makes this assumption physically reasonable: If an atom impinging on the surface is outside the range of interaction, there are
H. Henan and R.E. Allen
324
no forces and its occupancy is irrelevant. As Boon as it is within the interaction range, it is
interacting rather strongly, forming bonding and antibonding states with the solid, and its occupancy should be approximately correct. In Refs. 15 and 21, it is shown that the electronic force F:t, resulting from
Vorl
in (22)
and (26), can be simplified to F"•
F. '" ~Im {'" '" ~
= '" F· L..J'
(31)
;
aE;(E~) 10 az
[E! -
E;(E~)]}
(32)
gE!-E;(E~)
where the E~ are some set of !Sample energies in the complex-energy plane, along the contour C'J of Fig. 1 extending from (,.
+ ioo
down to the Fermi energy (,.. Also,
E!
and
E!
are
the lower and upper points of the interval along C, represented by E~. We have found that 5 such sample energies provide an accuracy of about 1% in the total energy and a few percent
in the electronic contribution to the force. The derivatives aE,jaz are calculated from the lIellman-Feynman theorem 1151
(33) Further details are given in Refs. 15 and 21.
ReE
Fig. 1. Complex energy plane. Integration over the rapidly varying structure between
f ..
(an
energy below the relevant electronic energy bands) a.nd t.,. (the Fermi energy of the solid) would require a very large number of sample energies for an accurate calculation. We
+ if mu and tr. Since the eigenvalues of the subspace Hamiltonian are 610wly varying functions of E for E complex, replace this interval by a. finite segment of C'J between t.,.
only a. few sample energies are required.
computer SImulatIons of ChemIsorptIon
325
Now let us turn to the results for various chemical species impinging on Ihe (UO) surface of GaAs. Figure 2 shows a top view of this surface, with the choice of coordinates indicated. Our unit of length is half the GaAs bond length, or 1.225
A. Then the As atoms at the lower
and upper right, respectively, have coordinates (1.63, -1.15, 0) and (1.63, 3.45, 0) before they are allowed to relax. Figure 2 also schematically shows a top view of 6 possible chemisorption sites on the surface, with bonding to a 'ingle As or Go. surface atom (lor 4), to a pair of As or Go. atoms (2 or 5), or in a Ga-As bridge .ite (3 or 6). These positions, which are observed as initial chemisorption sites in our simula.tions of time-dependent chemisorption, have also been observed as energy minima along the surface in calculations of energy surfaces associated with
static chemisorption [26,27J. . ,.) - - - - - (
GaAs(110)
Ga
0° ~
0
As,
0
CD
AS 2
Fig. 2. Top view or (110) surrace or GaAs, showing positions of Ga and As surface atoms and approxima.te positions of 6 initial chemisorption sites which have been observed in the present simulations.
In this paper, the unit oftime is the same as the time step,
t· =
t/~t.
~t
= 1.04 x 10-" sec; i.e.,
Also, in the simulations reported below, only two atoms are allowed to move-
the surface atom at the origin and the incoming atom. The other atoms in the semi-infinite solid affect the simulation, because they are electronically coupled to the moving atoms, and they are also represented by repulsive potentials. They are, however, motionless in these first simulations with the present technique. We have also performed some simulations with 4
atoms a.lIowed to move, and the differences are about what one expects. Namely, the results in most cases are fundamentally the same as those of the present paper, but the response of additional surface atoms does have some effect on the motion of the chemisorbing atom. One could also treat the propagation of vibrational energy into the solid with the time-
M. Menon and R. E. Allen
326
dependent phonon Green's function 117,28J. In the present work, however, we simply reduce each velocity by 0.5% at each time step. This estimate of the appropriate rate of energy dissipation is based on our study of vibrational relaxation in a semi-infinite one-dimensional chain of atoms, with a Budden impulse applied at the end 117,28J. In each case, the surface atoms are initially motionless at their relaxed positions. The initial position and velocity of
the incoming atom can be read off the graph.. Recall that the unit of length is r. /2, where r.
= 2.45A is the bulk GaA. bond length.
00 x -o . ~
-' 0
d--,,"oo"'-',"o'::o-::,o:::--:,oo o ::::--:,C: o"o '-',:o::o--::,oo :::--::!,oo
- 1.50
"
Fig. 3. Simulation of Cl chemisorbing on Go.Aa (110) . The coordinate system io defined in Fig. 2; also, t'
= t/ At, where At = 1.04 X 10-" sec.
Figure 3 shows a simulation for a Cl atom released with zero velocity one GaA. bond length above the surface, at the Ga-As bridge site 3 of Fig. 2. This atom moves away from the surface As, and bonds to a single Ga at site 4 of Fig. 2. Since CI is electronegative, this is the kind of behavior one expects, Notice that the z vibrations of the CI have large amplitude and a long period, because there is only a small, a.ngle-bending restoring force for motion in
that direction. We thus predict a low-frequency vibrational mode for this site. (Bonding to more distant neighbors will increase the frequency of this mode somewhat, but it should still be relatively low.) Notice that the present technique automatically gives angle-bending as well as bond-stretching forces, with both resulting from change. in the total electronic energy
Computer Slmulatlons of Chemlsorptlon
327
as the atomic positions are varied. Also notice that one can study the equilibrium vibra.tions
of surface and adsorbed atoms, even when such vibrations are highly anharmonic. Figure 4 shows & simula.tion for an 0 atom released with exactly the same initial conditions
as the CI above. Since 0 is light, it vibrates with high frequency in the II and
%
directions.
Again, however, it vibrates with a long period and large amplitude in the:z: direction, for which there is only a small angle-bending restoring force (in a nearest-neighbor approximation). 2.
0
2. N
,.
,.
•
0
-, - , ~o1--;:'OO;;;--;';;;OO;-;J;;;OO;-;"';;;;-"";;;OO~';;;O~O--;;'0-0-"""'" t'
....
Fig. 4. 0 on GaM (110). Figure 5 show a simulation for As, which bonds at the Ga-Ga bridge site 5. The initial conditions in the present simulations - zero velocity and the Ga-As bridge aite 3 - were chosen in order not to bias the motion of the chemisorbing atom; i.e .• it can choose to
bond either to the cation (Ga) or the anion (As). As one expects, the electronegative Group
V, VI, and VII atoms tend to bond to the cation. Furthermore, it will be seen below that the electropositive Group I, II, and III atoms tend to bond to the anion.
H. Menon and R.E . Allen
328 2.
....
1.
N
1. 1.
o.
,.,
o. o. -0.
-1. 1.
2.
"
1.
O.
a
l Oa
200
300
400
500
600
t'
700 8 00 GaAs
Fig. 5. A.!. For Group IV elements, one expects that Ga and As will be about equally attractive. This is confirmed by our simulations . In Fig. 6, Si is found to bond at the As-As bridge site 2. 51
2. N
1. 1.~=======:::=:::=====i
-o.s,
,.,
-1.
-1. -2.~============~ LOr o.
"
O. - 0. 0
Fig. 6. Si.
100
200
300
400
t·
500
600
700 Boe GaAa
Computer SImulatIons of ChemIsorptIon
329
In Figure 7, an AI atom is released 1.5 GaAs bond lengths above the surface, with an initial velocity downward. In this run alone, an As (rather than Ga) was placed at the origin. The AI bonds at site 4 with As at the origin, which is equivalent to site 1 with Ga at the origin. AI
·].i---:__
---J
o
100
200
300
400
500
600
t'
700
800
GoA. (Aol
Fig. 7. AI, with As at the origin.
2." . . . . - - - - - - - - - - - - ,
z.
o."~~========~
-o.5t;: -1 .
-0 ~o
Fig. 8. Zn.
100
200
100
400
t·
)00
600
100
800
CaA,
H. Henon and R.E. Allen
330
Figure 8 shows a simulation for Zn. As one expects, this electropositive metal atom is attracted by the anion, bonding at the As-As bridge .ite 2. Finally, Fig. 9 shows a simulation for Cu. Like Zn, this electropositive metal "tom bonds at the As-As bridge site 2.
2.
c.
2. N
,. 1. -0. -1.
>-
-, . -2.
o. x
o.
-O · 5!0--;-'07.0:-:: 20'-0-:'C:O-:: O -:.7. 00:-::'O'-O-:6C:O-::O~77. 00:-.::!OO t'
G.A.
Fig. 9. Cu. These results for various species chemisorbing on the (110) surface of GaA. are of conaiderable interest. Even more interesting, however, is the simple conclusion that the present method works, so that it is now possible to simulate atomic motion in real materials and at their surfa.ces.
Acknowledgment This work WM supported by the Office of N"val Research (NOOO14-82-K-0447) and by the Robert A. Welch Foundation.
1. 2. 3. 4.
J. J. J. J.
References Kouteay Phys. Rev. 108 (1957) 13. Koutecky, Trans. Faraday Soc. 54 (1958) 1038. Koutecky, Advan. Chern. Phys. !! (1965) 85. Koutecky, Progress in Surface "nd Membrane Science 11 (1976) 1.
Computer Simulations of Chemisorption
331
5. D. J. Chadi, Phy• . Rev. Lett. 41 (1978) 1062; Phy•. Rev. Lett. 52 (1984) 1911. 6. S. Y. Tong, A. R. Lubinsky, B. J. Mr.tik and M. A. van Hove, Phy•. Rev. Ill! (1978) 3303. 7. A. Kahn, E. So, P. Mark and C. B. Duke, J. Vac. Sci. Tech. 15 (1978) 580. 8. S. Y. Tong, G. Xu and W. N. Mei, Phys. Rev. Lett. 52 (1984) 1693. 9. R. P. Beres, R. E.Allen and J. D. Dow, Solid State Commun. 45 (1983) 13; Phys. Rev. ~ (1982) 769; Phys. Rev. B26 (1982) 5702. 10. G. P. Williams, R. J. Smith and G. J . Lapeyre, J. Vac. Sci. Tech. 15 (1978) 1249. 11. A. 11uijser, J. van Laar and T. L. van Rooy, Phys. Lett. 65A (1978) 337. 12. A. Ebina, T . Unno, Y. Suda, 11. Koinuma and T . Takahashi, J. Vac. Sci. Tech. 19 (1981) 301. 13. G. P. Srivastava, I. Singh, V. Montgomery and R . H. Williams, J. Phys. .Qlll (1983) 3627. 14. H. Carstensen, R. Manzke, I. Schafer and M. Skibowski, in Pree. Eighteenth Int. ConI. on the PhV.i.. 0/ Semicondudorl, Vol. 1edited by O. Engstrom (World Scientific, Singapore, 1987). 15. R. E. Allen and M. Menon, Phys. Rev. B33 (1986) 5611. 16. M. Menon and R. E. Allen, Phy•. Rev. B33 (1986) 7099. 17. M. Menon and R. E. Allen, Super lattices and Microstructures ~ (1987) 295. 18. M. Menon &nd R. E. Allen, Solid State Commun. M (1987) 353. 19. M. Menon and R. E. Allen, Proceedings 0/ the Second International Con,.,.nce on the Structure a' Sur'ace. (in press) . 20. M. Menon and R. E. Allen, Proceedings a' the 94th Annual Svmpo.ium the American
0'
21. 22. 23. 24. 25. 26. 27. 28.
Vacuum Societv, J. Vat. Sci. Tech. (in pres.) . M. Menon and R. E. Allen, to be published. R. Wang, J. Math. Phy•. 28 (1987) 2325,2329 and 2333. H. Lim and R. E. Allen, J. Vat. Sci. Tech. ~ (1985) 1221; 1M (1985) 2328. W. A. Harrison, Electroni. Strudu,. and the Properti.. 0/ Solids (W. H. Freeman, San Francisco, 1980). E. A. Kraut and W. A. Harrison, J. Vac. Sci. Tech. ID (1984) 409; R. Enderlein and W. A. Harrison, Phy •. Rev. B30 (1984) 1967. D. J. Chadi and R. Z. Bachrach, J . Vac. Sci. Tech. 16 (1979) 1159. J. lhm and J. D. Joannopoulos, Phy•. Rev. B26 (1982) 4429 &nd references therein. M. Menon &nd C. W. Myles, J. Phy•. Chem. Solids 48 (1987) 621.