- Email: [email protected]
Cluster calculations of the surface dimer structure on Si(100) surfaces
Surface Science 0 North-Holland
99 (1980) 581-597 Publishing Company
CLUSTER CALCULATIONS Si( 100) SURFACES W.S. VERWOERD
OF THE SURFACE DIMER STRUCTURE
ON
*
Institut fiir Theoretische Physik, Technische Universitli’t Clausthal, D-3392 Clausthal-Zellerfeld, Fed. Rep. of Germany
Received
17 March
1980
The MIND0/3 quantum chemistry procedure is used to calculate the structure of clusters having 9-15 silicon atoms, representing the first four layers and l-2 surface pairs of the Si(100) surface. We find an asymmetric dimer as the minimum energy configuration, and use a localised orbital analysis to show that it may be interpreted as the superposition of dimer bond formation and a secondary Jahn-Teller reconstruction. There is considerable interaction between neighbouring dimers, leading to a lowest energy surface arrangement having a 2 X 2 periodic structure. This is characterised by an overall relaxation of -0.29 A, a buckling of 0.24 A and a slightly stretched dimer bond.
1. Introduction The structure of the 2 X 1 reconstructed (100) surface of Si has aroused considerable interest in recent years, and apart from the previously well-established surface dimer model [ 1,2 ] and vacancy models [2,3] some new structural models have been proposed [4-61. The Levine dimer model [ 11, according to which row pairing of the first layer atoms takes place by a pair-wise linking up of dangling bonds (fig. lb), is intuitively appealing because partial bond saturation takes place while geometrical distortion is limited to a rotation of first to second layer bonds. It was also found by Appelbaum et al. [7] in their band structure calculations to give a good account of the surface density of states compared to ultraviolet photoemission (UPS) measurements, except that the splitting between the gap surface state bands are too small, giving rise to too much surface metallic character. Harrison has found [S] that the dangling bond electrons on the unrelaxed surface occupy a single doubly filled s-like state, instead of the separate “rabbit’s ears” dangling bonds suggested by a naive cleaved crystal model. This is supported by the
* Alexander von Humboldt research fellow. Permanent P.O. Box 392, Pretoria 0001, Rep. of South Africa.
581
address:
Department
of Physics,
Unisa,
W.S. Verwoerd / Cluster calculations
582
R R (b) 2x1
t
(cl 2 x 2
(d) 3x1
Fig. 1. The linking up of dangling bonds on the (100) surface terns with various periodicities: (a) unreconstructed surface, hypothetical bondlinking patterns.
of silicon, to form (b) Levine model,
surface pat(c) and (d)
charge density found by Appelbaum et al. [9] on the unrelaxed surface, as well as our own results discussed below. From this fact he concluded that a linking up of dangling bonds is unlikely to be energetically favourable because it would require promotion of electrons to form bonding hybrids. He therefore favours vacancy models which have atoms missing from the first layer. However, a direct comparison of the electronic spectra calculated for the vacancy and pairing models by Appelbaum et al. [7] shows that the pairing model fits the UPS data much better than the vacancy model. A more serious objection to the pairing model is that, at least in the simple form described above, it does not fit the results of LEED analyses [4]. This has led Jona et al. [4] to propose a complex chain-like structure, which may be seen as a linking up of the dimer pairs into a zig-zag chain while there are also substantial lateral shifts in the second layer, and which is found in a dynamical analysis to give a reasonable correspondence with the LEED data. The electronic structure of the chain model has been studied by Kerker et al. [lo], but is found to disagree with the UPS data, as a result of the large number of broken bonds and bond distortions involved. It was therefore suggested by Appelbaum and Hamman (AH) [ 1 l] on the grounds of a Keating elastic model calcula-
W.S. Verwoerd / Cluster calculations
583
tion that the pairing model is in fact correct, but that subsurface reconstruction of up to five layers deep must be included to obtain good correspondence with LEED data. Their kinematic analysis was followed up by a quasidynamical LEED analysis by Mitchell and Van Hove [12] which compared the chain model, the AH deep layer pairing model and some hexagonal close packed overlayer models due to Mitchell [S]. The conclusion reached is that only the pairing model and an overlayer model involving slight puckering, seems promising for further study. An even more satisfactory agreement with LEED 1-P’ curves was recently obtained by Tong and Maldonado [ 131 using the AH model with slightly modified overlayer spacings. While the available evidence therefore still seems to favour a pairing model of the reconstruction, its detailed structure remains uncertain. A direct determination of the atomic geometry by energy optimisation may therefore be useful to establish the input structure for LEED and band structure calculations. The results of such a calculation has recently been published by Chadi [6]. In his calculation, a Slater-Koster’ tight binding calculation of the electronic band energy is combined with a Keating-type elastic model of interatomic forces. This leads to a new variation of the pairing model, with subsurface displacements similar to those of the AH model but having an asymmetric dimer pair. The first layer thus becomes buckled in a similar way as is accepted for the (111) surface, and the metallic character inherent in the symmetric dimer is thereby removed. A less explicit suggestion of such weak secondary reconstruction superimposed on the dimerisation, was also previously made by Appelbaum, Baraff and Hamann [ 141. The approach presented in this article is also to determine the structure by energy minimisation, but this time an SCF quantum chemistry calculation is applied to an atomic cluster of Si atoms with the external bonds saturated with hydrogen. As will become clear below, the typical cluster size needed to study the dimerisation is too large for ab initio methods. However, our previous work on the Si( 111) surface [ 15 ,161 has shown that the semi-empirical MIND0/3 program [ 171 is quite suitable for such work and here we discuss its application to the Si(100) surface. The calculation of total energy is also appropriate for studying the question why a particular reconstruction pattern is observed instead of others which may also seem chemically plausible. For example, the 2 X 2 structure in fig. lc has the same number of bonds and seems as feasible as the standard Levine model of fig. lb, while the 3 X 1 triplet configuration of fig. Id has fewer dangling bonds and may seem even more favourable. A more detailed energy calculation is therefore needed to settle this question. Another feature of the work presented here is the use of localised molecular orbitals (LO’s) to interpret the reconstruction in terms of bond changes. The basic idea behind this analysis is to contruct spatially localised orbitals by making suitable linear combinations of the occupied molecular orbitals (MO’s). The two descriptions are equivalent in the sense that due to its determinantal structure, the many-electron wave function is left intact by such a linear transformation. How-
584
W.S. Verwoerd / Cluster calculations
ever, unlike the delocalised MO’s, the LO’s can each be associated with a specific bond or atom. A detailed discussion of the procedure used here to calculate the LO’s, as well as a demonstration of the one-to-one correspondence between intuitive chemical bonds and LO’s, has been published elsewhere [18]. The calculation itself is embodied in a Fortran program POPLOC, which has been deposited with the QCPE [19]. The same method has also been used in ref. [16] to analyse reconstruction on the (111) surface, where its relation to bond analysis by means of band structure calculations is discussed, and it is demonstrated that the definite hybridisation value obtained for each bond from an LO analysis enables one to make a rather detailed interpretation of the bond changes which take place. After a description of the atomic clusters and their optimisation in section 2, we apply such an analysis to the symmetric dimer model in section 3. It is found that bond formation between surface pairs does take place, but is also strongly influenced by interaction with neighbouring pairs. However, the main conclusion of this paper is discussed in section 4, namely that the minimum energy configuration is that of an asymmetric dimer, formed by superposition of a Jahn-Teller distortion and pair bonding. Different arrangements of the asymmetric pairs on the surface are possible, and we find the one corresponding to a 2 X 2 periodicity (fig. 7c) to be energetically most favourable. Finally we show that the electronic structure of all bonds in the first two atomic layers are considerably disturbed by the reconstruction. Since the geometric results of our calculation are quite similar to those of Chadi [6], it is also of interest to compare the methods of the two calculations in more detail. Based in part on a band structure calculation, the method of ref. [6] has the advantage of including the crystal periodicity and of avoiding the boundary condition problems of a cluster calculation. However, we believe our method to be more basic and more accurate for the following reasons. Firstly, the assumption of elastic behaviour with bulk parameter values by Chadi neglects the fact that the bonds near the surface may be substantially different, and our results in section 4 show that this does in fact happen. The effect of this difference may also be seen in a comparison of the two methods as applied to the Si( 111) surface. Cluster methods generally predict an inward relaxation of the 1 X 1 “unreconstructed” surface, our own calculation [ 161 yielding a value of 0.17 A which is in good agreement with experimental values, while Chadi [20] finds no relaxation. The same effect shows up on the buckled Si( 111) 2 X 1 surface: while both methods agree on the amount of buckling involved, our method [ 161 predicts a substantially larger overall relaxation but smaller second layer shifts (i.e. larger bond distortion) than found in ref. [20]. Secondly, Chadi does not include interelectron Coulomb repulsion, while intraatomic repulsion integrals are explicitly included [17] in the MIND0/3 program. This is particularly relevant here in view of the charge transfer that accompanies the buckling, and on the Si( 111) surface we have found indications that it is a major factor in determining the extent of buckling. Its importance is also supported
W.S. Verwoerd / Cluster calculations
585
by Chadi’s estimate [6] that Coulombic repulsion may amount to 20% of the total relaxation energy. Thirdly, unlike the band structure calculations of refs. [7] and [lo], the calculation of ref. [6] is not done self-consistently, and fixed bulk-derived Slater-Koster parameters are used. However, the very fact that the calculation predicts a substantial non-bulk-like interatomic charge transfer indicates the need to adapt the parameters self-consistently to the surface configuration. In the MIND0/3 method, on the other hand, the wave functions and Hamiltonian matrix are calculated fully self-consistently. This also implies that our calculation includes some substrate polarisation effects, at least within the cluster region, while it is completely neglected in ref. [6]. Finally, our method is also more fundamental in the sense that the establishment of a dimer bond is a result and not an a priori assumption of the calculation. The results of the present calculation and that of Chadi [6] differ in their predictions of the final agreement of dimers on the surface, and on the details of the charge distribution and dimer geometry. However, the agreement of two quite different calculations on the main result that dimer pairs on the Si( 100) surface are asymmetric, gives strong support to this basic idea.
2. Choice of atomic clusters For a realistic description of the dimerisation, we need at least the surface Si pair and their bonded neighbours in the second layer. A substitution of hydrogen atoms for all remaining bulk silicons at this stage would place H atoms in bridge positions where their valency would be insufficient to saturate both Si bonds. To avoid this problem we continue the Si subcluster down to the 4th layer, and then saturate all external bonds with hydrogen. Thus we arrive at the Si9Hrz cluster shown in fig. 2 to model the Si(100) surface. According to the philosophy also followed in our previous work, the first step in the calculation is to establish a corresponding bulk structure to which the surface relaxation may be compared. This is needed firstly, because as discussed extensively elsewhere [ 181, there are noticeable boundary effects associated with the proximity of H atoms to Si-Si bonds (at least in MIND0/3 calculations). Specifically, the bond length depends to some extent on the size of the cluster. Secondly, a bulk calculation is necessary to fix the optimal positions for the “bulk” hydrogen atoms. Our bulk cluster here is obtained by adding an SiHz group to saturate the “bridge” dangling bonds of the surface pair in fig. 2, and two H atoms to saturate their external bonds. The resulting cluster is simply the Si equivalent of the adamantane molecule. The bulk structure is obtained by restricting the Si subcluster to tetragonal symmetry, but with variable bond length, while the positions of all hydrogen atoms are freely optimised.
586
W.S. Verwoerd / Cluster calculations
Fig. 2. The SigHI cluster used to model layer to which they belong.
the surface.
Si atoms
are numbered
according
to the
All the hydrogen atoms representing bulk silicons, are kept fixed at these optimal positions for the surface calculation, while “1st layer” hydrogens are allowed to move freely together with the surface silicons. All subsurface silicons have also been kept fixed at bulk positions for the cluster under discussion since without all their crystal neighbours present realistic restrictions on their motion cannot otherwise be guaranteed. For the (111) surface [ 161 we had found it necessary to perform a secondary bulk calculation by allowing the bulk cluster to relax under the same constraints as those for the surface cluster. Here we find that (probably due to the larger cluster size) even if allowed, such “bulk relaxation” is negligible. For a study of the interaction between a surface pair and its neighbours on the surface, we also study two larger clusters. The first is an Si15H18 cluster with 3 surface atoms in a row, obtained by adding to the one in fig. 2 its mirror image in a plane parallel to the XZ plane but containing one of the surface atoms and its second layer neighbours. The second is a Si15H16 cluster with two adjacent 1st layer pairs, found by adding to fig. 2 its mirror image in an YZ plane which passes through two 2nd layer Si atoms. The calculations for these two larger surface clusters were preceded by that of their common bulk cluster.
3. The symmetric dimer model As a first step, we consider a symmetric relaxation in which one member of the surface pair in the Si9H,, cluster is allowed to relax freely, while the other is restricted to move symmetrically with respect to its (unrelaxed) bulk position. Since this is not the final minimum energy configuration, we limit ourselves in this section to a mainly qualitative discussion.
W.S. Verwoerd / Cluster calculations
5.07
I-
3.87
II-
587
I )
2.75
Fig. 3. Inwardly and outwardly relaxed configurations of a single dimer pair. All lengths shown in A. Inset shows the relative energies of the unrelaxed and two relaxed positions.
are
Fig. 3 illustrates the somewhat surprising result obtained, that there are two different relaxed positions: either an inwardly or outwardly relaxed configuration is attained, depending on the initial disturbances of the metastable unrelaxed position. The reason for this behaviour becomes clear from the localised orbital analysis shown in fig. 4. For the unrelaxed surface (fig. 4a) we find each surface atom to have a doubly occupied, largely s-like dangling bond orbital accompanied by two very p-like hybrids in the back bonds (8% s-content, not shown on the figure). These atoms are thus approximately in the atomic ground state, and do not have “rabbit’s ears” dangling bonds. This agrees with the results of Harrison [8], and his conjecture that such orbitals will not interact by bond formation because this requires promotion into bonding hybrids, agrees with the electronic configuration of the outwardly relaxed pair. It shows (fig. 4b) a repulsive interaction between the two orbitals, but their composition, as well as that of the backbonds, is essentially the same as for the unrelaxed case. In contrast to this, for inward relaxation (fig. 4c) we do find the establishment of a pair bond between the atoms, with some electronic promotion indicated by an s-content decrease from 84% to 70% in the bonding hybrids. The orbitals of the non-bonding electrons stay unchanged, but are now singly occupied on each atom since they jointly make up a doubly occupied LO. We note that the pair bond in this case is not a typical covalent Si-Si bond like that shown for a bulk bond in fig. 4a, as seen both from the overlap geometry and the relatively high s-content still present, and at 2.75 a is also much longer than the bulk length. However, for modelling the surface we must also take into account that mutual
588
W.S. Verwoerd / Cluster calculations
b
I
Fig. 4. Localised orbitals (LO’s) calculated for the single pair cluster: (a) unrelaxed, (b) outwardly relaxed and (c) inwardly relaxed. Outlines are curves of constant wave function amplitude $ = 0.21 on each atom, a value arbitrarily chosen to give a good overlap for typical bulklike bonds like that shown in fig. 4a. Each LO has a total occupation of two electrons and is represented by a different kind of shading. The percentage s content of each atomic hybrid is shown and a reduced scale has been used perpendicular to the hybrid direction to clarify the directionality.
W.S. Verwoerd / Cluster calculations
589
approach of any two surface atoms moves them away from their outer neighbours along the dimer axis, and according to the results discussed above this will generate a repulsive interaction with compressive effects on the pair. To study such effects, we therefore add another surface atom, arriving at the SirsHI cluster discussed in section 2. The first calculation performed with this cluster is to let the two outer surface atoms approach the centre one symmetrically, as expected for the hypothetical 3 X 1 pattern of fig. Id. We find that in the optimum configuration under this restriction, the energy gain per atom is considerably less than for pairing, there is little evidence for bond formation (see fig. Sa) and if the symmetry restriction is lifted the system is unstable against the formation of a pair while the remaining atom is repelled. We believe that this result is, also of relevance for a consideration of the chain model of 2 X 1 reconstruction mentioned in the introduction. The chain model may also be interpreted as a linking of the external dangling bonds of a pair to other atoms, but involves much larger distortions of the back bonds than the 3 X 1 pattern. Since we find the latter to be energetically unfavourable, it is difficult to believe that the chain model would be feasible. The final geometry attained when one surface atom is allowed to move freely while the other two are coupled to it by the restrictions of the symmetric dimer model, is shown in fig. 5b, and the corresponding electronic structure in fig. SC. Comparison with the case of a single pair shows the additional neighbour to have a substantialleffect. There is a large decrease in the dimer bond length, making it even shorter than the bulk bonds. The electronic structure now shows a welldeveloped bulk-like dimer bond, while the dangling bonds have become largely p-like and their orientation suggests that some n-bonding between them takes place. This may explain the short bond length, and evidence for such n-bonding in the symmetric dimer model has also been found by Appelbaum, Baraff and Hamann [7,14:] from their band structure calculations, leading them to a bond length proposal of 2.22 a very similar to our optimised value. The result above is certainly not conclusive about the dimer state on a surface, since only one neighbour is included and furthermore its electronic configuration is different from that of a pair member. Nevertheless, it does show that there can be a surprisingly large interaction between surface pairs despite their separation of about 5.6 A. Also shown on fig. Sb are some 2nd and 3rd layer relaxations. This was obtained by allowing the indicated atoms freedom of motion subject to the 2 X 1 symmetry of the first layer. The displacements that we find agree in direction with those found from elastic models by Appelbaum and Hamann [ 1 l] and Chadi [6], but especially our second layer movements are much smaller. Part of the difference may be ascribed to the fact that our cluster does not include the full nearest neighbour complement for subsurface atoms, and that it is too small to include 4th layer relaxation. However, we believe that it is also partly due to an overestimate of the
590
W.S. Verwoerd / Cluster calculations
back bond strain by elastic models because of the assumption that these are identical to bulk bonds. This point will be discussed further in section 4. The subsurface relaxation also plays a role in distinguishing between the 2 X 1 and 2 X 2 reconstruction patterns of figs. lb and lc. If we define a dimer plane (for each dimer pair) containing the pair bond and the [loo] axis, it is easily seen that the subsurface relaxations induced by pairs in neighbouring dimer planes will be mutually reinforcing for the 2 X 1 pattern but opposing for the 2 X 2 structure, so that the former will be energetically more favourable. Seen dynamically, the pairing of two atoms on the unreconstructed surface will propagate longitudinally
a Ay=O.832 AZ =-0.356
W.S. Verwoerd / Cluster calculations
591
Fig. 5. Localised orbitals and geometry of the relaxed triplet cluster: (a) 3 X 1 configuration, (b) and (c) symmetric 2 x 1. Lengths are measured in A and LO’s are drawn to the same conventions as in fig. 4.
(in t.he dimer plye) by the repulsive interaction discussed previously, and transversely (perpendicular to the dimer plane) by the fact that a neighbouring surface pair will tend to be pushed together by the relaxation that the first pair induces in the second layer.
4. Asymmetric dimerisation We now turn to a free optimisation of the single pair SigHI cluster, allowing both surface atoms to move independently. This results in an asymmetric dimer configuration, shown in fig. 6a. Whereas a symmetric dimer can be specified by the dimer bond length d and the overall surface layer relaxation 6, the asymmetry gives rise in addition to a lateral shift h of the centre of gravity of each pair, and a buckling B of the surface layer. Parameter values are listed in table 1 and will be discussed below. First, however, we consider the reconstruction mechanism. A comparison of the localised orbitals of the symmetric pair (figs. 4c and SC) and the asymmetric case (figs. 6b and 6c) shows that an electron has been transferred from one of the degenerate singly occupied dangling bond orbitals of the symmetric pair to the other, to give a doubly occupied orbital on the “raised” member of the asymmetric dimer. This suggests a Jahn?Teller mechanism. To substantiate it, one must confirm that the degeneracy was lifted by the geometric distortion, and this is done by noting
592
WS. Verwoerd / Clusrer calculations
b
I
Fig. 6. Asymmetric dimerisations: (a) geometric parameters characterising the reconstruction; (b) LO’s of a single pair, drawn as in fig. 4; Cc) LO’s of the triplet, dangling bonds on the outer atoms.
the large difference in hybridisation of the two dangling bonds in fig. 6c. Almost all of the electron occupation is carried by the s-like and therefore low-energy orbital while the p-like orbital is &most empty. The same applies to fig, 6b, but here the occupation of the p-like orbital is too small to be seen in the Lo’s, which by
W.S. Verwoerd / Cluster calculations
593
Table 1 The reconstructed surface layer atomic geometry, all entries in A; 6 represents an overall relaxation and h a lateral shift of the dimer pair centre of gravity, the buckling 6 is the deviation along [ 1001 from this mean position, and d is the dimer bond length Coord.
6
h
b
d
Ref.
AH Levine Chadi Pair 2x1 2x2
-0.09 -0.22 -0.20 -0.25 -0.31 -0.29
0.0 0.0 0.31 0.20 0.23 0.21
0.0 0.0 0.24 0.24 0.24 0.24
2.45 2.35 2.35 2.46 2.40 2.42
[l&13] [ll
161 Present work Present work Present work
definition are constructed from occupied orbitals. The change in dangling bond hybridisation on each atom is accompanied by an opposite change in its other bonds, as clearly seen in the dimer bond in fig. 6, and this is also true of the back bonds (see table 2). According to a well known argument that bonds formed by s-like hybrids tend to be shorter than those that are p-like [ 181 one would therefore expect differences in the back bond lengths, and indeed the “donor” atom has a back bond length of only 2.32 a compared to the bulk cluster value of 2.37 a and the “acceptor” atom’s 2.42 a. While we have previously found evidence [ 181 that the MIND0/3 program exaggerates the effect of hybridisation on bond lengths, it seems reasonable to conclude that the connection between geometrical distortions and the electron transfer that make up the Jahn-Teller effect, is brought about here by accompanying hybridisation and bond length changes in the dimer and back bonds. Table 2 also illustrates the penetration of surface effects on the electronic structure into the cluster. It is seen that even the considerable charge redistribution
Table 2 The hybridisation @OScontent) on each atom of bonds in the SisH 12 cluster Layer
Bulk
Unrelaxed
Symmetric
Asymmetric
-
-
_
-
70 70
79 3
45 47
8 51
54 47
1‘3 43
56 44
48 45
4.5 46
40 44
52 47
38 48
45 47
47 47
46 47
45 47
45 47
594
W.S. Verwoerd / Cluster calculations
associated with the asymmetric reconstruction stage, only seriously affects bonds in which the first layer atoms participate. This justifies the use of elastic models for deeper bonds as done in refs. [6] and [ 111, but we find that the differences in the surface bonds lead to smaller deep layer relaxations for the asymmetric dimer. The increase in dimer and back bond s content of the donor atom may be expected to increase its electronegativity [2 11, thus acquiring more electrons (while the reverse is true for the acceptor) and so counteracting the effect of the dangling bond charge transfer. This effect is visible in the dimer bond of fig. 6 which shows the s-like hybrid of the donor atom to be slightly larger than its partner, and is reflected in the fact that we find the nett charge transfer between the two atoms to be considerably smaller than a full electron, namely 0.221el. The latter figure is derived from the calculated nett atomic charges on the lowered and raised atoms, namely 0.32lel and -0.12lel respectively. Taking into account also the atomic charges on the unrelaxed surface (0.191el) and on the symmetric pair (-0.05lel) the picture that emerges is that 0.24 electrons per atom is transferred to the surface layer during formation of the dimer bond, while in the Jahn-Teller stage the donor loses 0.37 electrons of which -0.22lel is transferred to its pair partner, while 0.15 electrons is gained from each atom by the substrate. Our value for the pair charge transfer is considerably smaller than the 0.361el value of Chadi [6], as is to be expected from our inclusion of Coulombic effects. On the other hand, it is quite similar to the -0.20lel that we have found for Jahn-Teller reconstruction of the (111) surface [ 161. The geometrical results of our single pair configuration is shown in table 1 to be quite similar to those of Chadi, with an identical value predicted for the buckling. However, our calculation gives more “vertical” relaxation and less lateral displacement of the pair, implying bigger back bond changes. Also notable is that we find the dimer bond to be somewhat longer than the bulk bond length, which agrees with the value used in recent LEED analyses [ 1 l-l 31. The asymmetry of dimer pairs gives rise to various possible arrangements on the surface with different periodicities, even assuming the basic pairing to follow the 2 X 1 pattern. In the dimer plane (containing the dimer bond and the [loo] axis) it can be either two-fold or four-fold in terms of the unreconstructed mesh, and perpendicular to this plane have either single or double period. Some of the resulting configurations are shown in fig. 7. In order to assess which of these have the lowest energy, we next have to evaluate the interaction between dimer pairs. First, we consider the interaction between a dimer and its nearest neighbour in the direction perpendicular to the dimer plane. As described in section 2, we use a 15-atom cluster containing two such pairs for this purpose, and find that the reconstruction energy decrease of the double period configuration (crossed dimers) is at A,!? = -0.77 eV significantly more favourable than the single period (parallel dimer axes) case with AZ? = -0.70 eV (all energies compared with the unreconstructed surface). As a basis of comparison, it may be noted that for a single pair we find
W.S. Verwoerd /Cluster calculations
595
(b) 4 xl
(cl
2 x2
Fig. 7. Different periodicities resulting from the superposition the basic 2 X 1 dimerisation pattern.
AE = -0.61 eV. By comparison
of Jahn-Teller
reconstruction
on
of the symmetric and asymmetric pairs, we find that roughly -0.2 eV may be ascribed to the Jahr-Teller phase and -0.4 eV to the initial symmetric dimer bond formation. The discrepancy between the latter figure and the typical Si-Si bond energy of 2.0 eV [15] is mainly due to the angular distortion of the back bonds, as may be justified by noting that the strain energy calculated from the empirical bond force constants used by Harrison [3] amounts to about 1.7 eV. Secondly, the interaction of a dimer with its neighbours in the dimer plane has to be considered. Unfortunately the cluster size for a 4-layer double pair cluster becomes excessive in this case. However, inspection of fig. 6c shows that the dangling bond orbital of the third atom in the previously discussed triplet cluster is actually quite similar to that of the dimer acceptor atom. The interaction between doubly occupied dangling bonds was shown above to be substantial, and we try to model at least this part of the interaction between dimers by using the third atom of the triplet clusteras representative of a neighbouring dimer. It should be noted, however, that the nett atomic charge on this third atom is +0.2lel instead of
596
W.S. Verwoerd 1 Cluster calculations
-0.12lel (mainly because it lacks participation in a dimer bond) for which an electrostatic energy correction must be added to the triplet results. We find that the energy of the triplet with occupied dangling bonds on the outer atoms (as for periodicity of 2) is -0.97 eV compared to -0.95 eV when they are on adjacent atoms (periodicity of 4). From the energy of the outwardly relaxed pair (fig. 3b) we derive an energy contribution for the relaxation of the third atom only, and subtracting this and the electrostatic energy correction, we find reconstruction energies of -1.0 eV and -0.90 eV per dimer respectively for the two configurations. It should be noted that this does not include the electrostatic interaction between dimer pairs, which we find to increase the difference between the two configurations by an additional 0.09 eV. The point charge electrostatic interaction thus contributes roughly half of the total interaction energy between neighbouring dimers in this plane, and this is also found to apply perpendicularly to the dimer plane. Combining these results, one may now compare the dimer energies of the three reconstruction patterns in fig. 7 and of the 4 X 2 pattern formed by combining the periodicities in figs. 7b and 7c. Including the interaction of a dimer with its four nearest neighbour pairs (one on each side, both in the dimer plane and perpendicular to it) we find an energy of -1.25 eV for the 2 X 2 pattern of fig. 7c, -1 .l 1 eV for both the 2 X 1 and the 4 X 2 patterns and -0.97 eV for the 4 X 1 pattern. Structural parameters for the 2 X 1 and 2 X 2 arrangements (arrived at by interpolation of the triplet and double pair results) are shown in the last two columns of table 1. The actual energy values mentioned here are probably not very reliable, especially as a result of the rather indirect approach that we had to use for interactions between dimers in the same plane. However, we believe that the energy margin favouring the 2 X 2 structure is sufficient to qualify it as a serious contender for the interpretation of LEED and other experimental results. The various arrangements of asymmetric dimers on the surface have also been considered by Chadi [6]. While he mentions a variety of experimental indications in favour of a 2 X 2 pattern, his final conclusion is that the 4 X 2 structure has the lowest energy, mainly as the result of a difference in surface Madelung energy. Consequently we have also investigated the effect of including electrostatic interactions with the more remote dimer pairs, but find the differences between their contributions for different arrangements to be too small to affect our result in favour of the 2 X 2 pattern seriously. The discrepancy between our conclusion and that of Chadi, is probably due mainly to the smaller dimer dipole moment that we find, and to the fact that we find substantial nearest neighbour interactions which are not of a point charge electrostatic nature and which already favours the 2 X 2 structure. Another factor which might contribute, is the subsurface relaxation. Both from our own results and those of Chadi it appears that this contribution to the dimer relaxation energy is of the order of 0.1 eV, which is comparable to the energies separating the various patterns. However, our calculations indicate that this is
W.S. Verwoerd / Cluster calculations
597
related mainly to the dimer bond formation stage, whch is also supported by the very similar subsurface relaxations reported by AH [ 1 I] for the symmetric and Chadi [6] for the asymmetric dimers. Since the surface layer patterns under discussion only differ in the way the asymmetric part of the reconstruction is arranged, we therefore do not expect subsurface relaxation to play an important role in distinguishing between them.
Acknowledgements The author gratefully acknowledges financial support for computer time received from the Council for Scientific and Industrial Research, Pretoria. He is also indebted to the Alexander von Humboldt foundation for the receipt of a fellowship while the later stages of this work was completed.
[l] J.P. Levine, Surface Sci. 34 (1973) 90. [ 21 R.E. Schlier and H.E. Farnsworth, Semiconductor Surface Physics (Univ. of Pennsylvania Press, Philadelphia, PA, 1957). [3] W.A. Harrison, Surface Sci. 55 (1976) 1. [4] I’. Jona, H.D. Shih, A. Ignatiev, D.W. Jepsen and P.M. Marcus, J. Phys. Cl0 (1977) L67. [5] K.A.R. Mitchell, Surface Sci. 64 (1977) 797. [6] D.J. Chadi, Phys. Rev. Letters 43 (1979) 43. [7] J.A. Appelbaum, G.A. Baraff and D.R. Hamann, Phys. Rev. Letters 35 (1976) 729; Phys. Rev. B14 (1976) 588; B15 (1977) 2408. [8] W.A. Harrison, Surface Sci. 55 (1976) 1. [9] J.A. Appelbaum, G.A. Baraff and D.R. Hamann, Phys. Rev. Bll (1975) 3822. [lo] G.P. Kerker, S.G. Louie and M.L. Cohen, Phys. Rev. B17 (1978) 706. [ll] J.A. Appelbaum and D.R. Hamann, Surface Sci. 74 (1978) 21. [ 121 K.A.R. Mitchell and M.A. Van Hove, Surface Sci. 75 (1978) 147L. [13] S.V. Tong and A.L. Maldonado, Surface Sci. 78 (1978) 459. [14] J.A. Appelbaum, G.A. Baraff and D.R. Hamann, Phys. Rev. Letters 36 (1976) 450. [15] W.S. Verwoerd and F.J. Kok, Surface Sci. 80 (1979) 89. [ 161 El. Gaigher and W.S. Verwoerd, Surface Sci., submitted. [17] R.C. Bingham,M.J.S. Dewar and D.H. Lo, J. Am. Chem. Sot. 97 (1975) 1285; M.J.S. Dewar et al., QCPE 10 (1976) 309. [18] W.S. Verwoerd, Chem. Phys. 44 (1979) 151; W.S. Verwoerd, Intern. J. Quantum Chem., to be published. [ 191 Quantum Chemistry Program Exchange, Indiana University, Bloomington, IN 47401, USA. [20] D.J. Chadi, Phys. Rev. Letters 41 (1978) 1062. [21] J. Hinze and H.H. Jaff6, J. Am. Chem. Sot. 84 (1962) 540.
99 (1980) 581-597 Publishing Company
CLUSTER CALCULATIONS Si( 100) SURFACES W.S. VERWOERD
OF THE SURFACE DIMER STRUCTURE
ON
*
Institut fiir Theoretische Physik, Technische Universitli’t Clausthal, D-3392 Clausthal-Zellerfeld, Fed. Rep. of Germany
Received
17 March
1980
The MIND0/3 quantum chemistry procedure is used to calculate the structure of clusters having 9-15 silicon atoms, representing the first four layers and l-2 surface pairs of the Si(100) surface. We find an asymmetric dimer as the minimum energy configuration, and use a localised orbital analysis to show that it may be interpreted as the superposition of dimer bond formation and a secondary Jahn-Teller reconstruction. There is considerable interaction between neighbouring dimers, leading to a lowest energy surface arrangement having a 2 X 2 periodic structure. This is characterised by an overall relaxation of -0.29 A, a buckling of 0.24 A and a slightly stretched dimer bond.
1. Introduction The structure of the 2 X 1 reconstructed (100) surface of Si has aroused considerable interest in recent years, and apart from the previously well-established surface dimer model [ 1,2 ] and vacancy models [2,3] some new structural models have been proposed [4-61. The Levine dimer model [ 11, according to which row pairing of the first layer atoms takes place by a pair-wise linking up of dangling bonds (fig. lb), is intuitively appealing because partial bond saturation takes place while geometrical distortion is limited to a rotation of first to second layer bonds. It was also found by Appelbaum et al. [7] in their band structure calculations to give a good account of the surface density of states compared to ultraviolet photoemission (UPS) measurements, except that the splitting between the gap surface state bands are too small, giving rise to too much surface metallic character. Harrison has found [S] that the dangling bond electrons on the unrelaxed surface occupy a single doubly filled s-like state, instead of the separate “rabbit’s ears” dangling bonds suggested by a naive cleaved crystal model. This is supported by the
* Alexander von Humboldt research fellow. Permanent P.O. Box 392, Pretoria 0001, Rep. of South Africa.
581
address:
Department
of Physics,
Unisa,
W.S. Verwoerd / Cluster calculations
582
R R (b) 2x1
t
(cl 2 x 2
(d) 3x1
Fig. 1. The linking up of dangling bonds on the (100) surface terns with various periodicities: (a) unreconstructed surface, hypothetical bondlinking patterns.
of silicon, to form (b) Levine model,
surface pat(c) and (d)
charge density found by Appelbaum et al. [9] on the unrelaxed surface, as well as our own results discussed below. From this fact he concluded that a linking up of dangling bonds is unlikely to be energetically favourable because it would require promotion of electrons to form bonding hybrids. He therefore favours vacancy models which have atoms missing from the first layer. However, a direct comparison of the electronic spectra calculated for the vacancy and pairing models by Appelbaum et al. [7] shows that the pairing model fits the UPS data much better than the vacancy model. A more serious objection to the pairing model is that, at least in the simple form described above, it does not fit the results of LEED analyses [4]. This has led Jona et al. [4] to propose a complex chain-like structure, which may be seen as a linking up of the dimer pairs into a zig-zag chain while there are also substantial lateral shifts in the second layer, and which is found in a dynamical analysis to give a reasonable correspondence with the LEED data. The electronic structure of the chain model has been studied by Kerker et al. [lo], but is found to disagree with the UPS data, as a result of the large number of broken bonds and bond distortions involved. It was therefore suggested by Appelbaum and Hamman (AH) [ 1 l] on the grounds of a Keating elastic model calcula-
W.S. Verwoerd / Cluster calculations
583
tion that the pairing model is in fact correct, but that subsurface reconstruction of up to five layers deep must be included to obtain good correspondence with LEED data. Their kinematic analysis was followed up by a quasidynamical LEED analysis by Mitchell and Van Hove [12] which compared the chain model, the AH deep layer pairing model and some hexagonal close packed overlayer models due to Mitchell [S]. The conclusion reached is that only the pairing model and an overlayer model involving slight puckering, seems promising for further study. An even more satisfactory agreement with LEED 1-P’ curves was recently obtained by Tong and Maldonado [ 131 using the AH model with slightly modified overlayer spacings. While the available evidence therefore still seems to favour a pairing model of the reconstruction, its detailed structure remains uncertain. A direct determination of the atomic geometry by energy optimisation may therefore be useful to establish the input structure for LEED and band structure calculations. The results of such a calculation has recently been published by Chadi [6]. In his calculation, a Slater-Koster’ tight binding calculation of the electronic band energy is combined with a Keating-type elastic model of interatomic forces. This leads to a new variation of the pairing model, with subsurface displacements similar to those of the AH model but having an asymmetric dimer pair. The first layer thus becomes buckled in a similar way as is accepted for the (111) surface, and the metallic character inherent in the symmetric dimer is thereby removed. A less explicit suggestion of such weak secondary reconstruction superimposed on the dimerisation, was also previously made by Appelbaum, Baraff and Hamann [ 141. The approach presented in this article is also to determine the structure by energy minimisation, but this time an SCF quantum chemistry calculation is applied to an atomic cluster of Si atoms with the external bonds saturated with hydrogen. As will become clear below, the typical cluster size needed to study the dimerisation is too large for ab initio methods. However, our previous work on the Si( 111) surface [ 15 ,161 has shown that the semi-empirical MIND0/3 program [ 171 is quite suitable for such work and here we discuss its application to the Si(100) surface. The calculation of total energy is also appropriate for studying the question why a particular reconstruction pattern is observed instead of others which may also seem chemically plausible. For example, the 2 X 2 structure in fig. lc has the same number of bonds and seems as feasible as the standard Levine model of fig. lb, while the 3 X 1 triplet configuration of fig. Id has fewer dangling bonds and may seem even more favourable. A more detailed energy calculation is therefore needed to settle this question. Another feature of the work presented here is the use of localised molecular orbitals (LO’s) to interpret the reconstruction in terms of bond changes. The basic idea behind this analysis is to contruct spatially localised orbitals by making suitable linear combinations of the occupied molecular orbitals (MO’s). The two descriptions are equivalent in the sense that due to its determinantal structure, the many-electron wave function is left intact by such a linear transformation. How-
584
W.S. Verwoerd / Cluster calculations
ever, unlike the delocalised MO’s, the LO’s can each be associated with a specific bond or atom. A detailed discussion of the procedure used here to calculate the LO’s, as well as a demonstration of the one-to-one correspondence between intuitive chemical bonds and LO’s, has been published elsewhere [18]. The calculation itself is embodied in a Fortran program POPLOC, which has been deposited with the QCPE [19]. The same method has also been used in ref. [16] to analyse reconstruction on the (111) surface, where its relation to bond analysis by means of band structure calculations is discussed, and it is demonstrated that the definite hybridisation value obtained for each bond from an LO analysis enables one to make a rather detailed interpretation of the bond changes which take place. After a description of the atomic clusters and their optimisation in section 2, we apply such an analysis to the symmetric dimer model in section 3. It is found that bond formation between surface pairs does take place, but is also strongly influenced by interaction with neighbouring pairs. However, the main conclusion of this paper is discussed in section 4, namely that the minimum energy configuration is that of an asymmetric dimer, formed by superposition of a Jahn-Teller distortion and pair bonding. Different arrangements of the asymmetric pairs on the surface are possible, and we find the one corresponding to a 2 X 2 periodicity (fig. 7c) to be energetically most favourable. Finally we show that the electronic structure of all bonds in the first two atomic layers are considerably disturbed by the reconstruction. Since the geometric results of our calculation are quite similar to those of Chadi [6], it is also of interest to compare the methods of the two calculations in more detail. Based in part on a band structure calculation, the method of ref. [6] has the advantage of including the crystal periodicity and of avoiding the boundary condition problems of a cluster calculation. However, we believe our method to be more basic and more accurate for the following reasons. Firstly, the assumption of elastic behaviour with bulk parameter values by Chadi neglects the fact that the bonds near the surface may be substantially different, and our results in section 4 show that this does in fact happen. The effect of this difference may also be seen in a comparison of the two methods as applied to the Si( 111) surface. Cluster methods generally predict an inward relaxation of the 1 X 1 “unreconstructed” surface, our own calculation [ 161 yielding a value of 0.17 A which is in good agreement with experimental values, while Chadi [20] finds no relaxation. The same effect shows up on the buckled Si( 111) 2 X 1 surface: while both methods agree on the amount of buckling involved, our method [ 161 predicts a substantially larger overall relaxation but smaller second layer shifts (i.e. larger bond distortion) than found in ref. [20]. Secondly, Chadi does not include interelectron Coulomb repulsion, while intraatomic repulsion integrals are explicitly included [17] in the MIND0/3 program. This is particularly relevant here in view of the charge transfer that accompanies the buckling, and on the Si( 111) surface we have found indications that it is a major factor in determining the extent of buckling. Its importance is also supported
W.S. Verwoerd / Cluster calculations
585
by Chadi’s estimate [6] that Coulombic repulsion may amount to 20% of the total relaxation energy. Thirdly, unlike the band structure calculations of refs. [7] and [lo], the calculation of ref. [6] is not done self-consistently, and fixed bulk-derived Slater-Koster parameters are used. However, the very fact that the calculation predicts a substantial non-bulk-like interatomic charge transfer indicates the need to adapt the parameters self-consistently to the surface configuration. In the MIND0/3 method, on the other hand, the wave functions and Hamiltonian matrix are calculated fully self-consistently. This also implies that our calculation includes some substrate polarisation effects, at least within the cluster region, while it is completely neglected in ref. [6]. Finally, our method is also more fundamental in the sense that the establishment of a dimer bond is a result and not an a priori assumption of the calculation. The results of the present calculation and that of Chadi [6] differ in their predictions of the final agreement of dimers on the surface, and on the details of the charge distribution and dimer geometry. However, the agreement of two quite different calculations on the main result that dimer pairs on the Si( 100) surface are asymmetric, gives strong support to this basic idea.
2. Choice of atomic clusters For a realistic description of the dimerisation, we need at least the surface Si pair and their bonded neighbours in the second layer. A substitution of hydrogen atoms for all remaining bulk silicons at this stage would place H atoms in bridge positions where their valency would be insufficient to saturate both Si bonds. To avoid this problem we continue the Si subcluster down to the 4th layer, and then saturate all external bonds with hydrogen. Thus we arrive at the Si9Hrz cluster shown in fig. 2 to model the Si(100) surface. According to the philosophy also followed in our previous work, the first step in the calculation is to establish a corresponding bulk structure to which the surface relaxation may be compared. This is needed firstly, because as discussed extensively elsewhere [ 181, there are noticeable boundary effects associated with the proximity of H atoms to Si-Si bonds (at least in MIND0/3 calculations). Specifically, the bond length depends to some extent on the size of the cluster. Secondly, a bulk calculation is necessary to fix the optimal positions for the “bulk” hydrogen atoms. Our bulk cluster here is obtained by adding an SiHz group to saturate the “bridge” dangling bonds of the surface pair in fig. 2, and two H atoms to saturate their external bonds. The resulting cluster is simply the Si equivalent of the adamantane molecule. The bulk structure is obtained by restricting the Si subcluster to tetragonal symmetry, but with variable bond length, while the positions of all hydrogen atoms are freely optimised.
586
W.S. Verwoerd / Cluster calculations
Fig. 2. The SigHI cluster used to model layer to which they belong.
the surface.
Si atoms
are numbered
according
to the
All the hydrogen atoms representing bulk silicons, are kept fixed at these optimal positions for the surface calculation, while “1st layer” hydrogens are allowed to move freely together with the surface silicons. All subsurface silicons have also been kept fixed at bulk positions for the cluster under discussion since without all their crystal neighbours present realistic restrictions on their motion cannot otherwise be guaranteed. For the (111) surface [ 161 we had found it necessary to perform a secondary bulk calculation by allowing the bulk cluster to relax under the same constraints as those for the surface cluster. Here we find that (probably due to the larger cluster size) even if allowed, such “bulk relaxation” is negligible. For a study of the interaction between a surface pair and its neighbours on the surface, we also study two larger clusters. The first is an Si15H18 cluster with 3 surface atoms in a row, obtained by adding to the one in fig. 2 its mirror image in a plane parallel to the XZ plane but containing one of the surface atoms and its second layer neighbours. The second is a Si15H16 cluster with two adjacent 1st layer pairs, found by adding to fig. 2 its mirror image in an YZ plane which passes through two 2nd layer Si atoms. The calculations for these two larger surface clusters were preceded by that of their common bulk cluster.
3. The symmetric dimer model As a first step, we consider a symmetric relaxation in which one member of the surface pair in the Si9H,, cluster is allowed to relax freely, while the other is restricted to move symmetrically with respect to its (unrelaxed) bulk position. Since this is not the final minimum energy configuration, we limit ourselves in this section to a mainly qualitative discussion.
W.S. Verwoerd / Cluster calculations
5.07
I-
3.87
II-
587
I )
2.75
Fig. 3. Inwardly and outwardly relaxed configurations of a single dimer pair. All lengths shown in A. Inset shows the relative energies of the unrelaxed and two relaxed positions.
are
Fig. 3 illustrates the somewhat surprising result obtained, that there are two different relaxed positions: either an inwardly or outwardly relaxed configuration is attained, depending on the initial disturbances of the metastable unrelaxed position. The reason for this behaviour becomes clear from the localised orbital analysis shown in fig. 4. For the unrelaxed surface (fig. 4a) we find each surface atom to have a doubly occupied, largely s-like dangling bond orbital accompanied by two very p-like hybrids in the back bonds (8% s-content, not shown on the figure). These atoms are thus approximately in the atomic ground state, and do not have “rabbit’s ears” dangling bonds. This agrees with the results of Harrison [8], and his conjecture that such orbitals will not interact by bond formation because this requires promotion into bonding hybrids, agrees with the electronic configuration of the outwardly relaxed pair. It shows (fig. 4b) a repulsive interaction between the two orbitals, but their composition, as well as that of the backbonds, is essentially the same as for the unrelaxed case. In contrast to this, for inward relaxation (fig. 4c) we do find the establishment of a pair bond between the atoms, with some electronic promotion indicated by an s-content decrease from 84% to 70% in the bonding hybrids. The orbitals of the non-bonding electrons stay unchanged, but are now singly occupied on each atom since they jointly make up a doubly occupied LO. We note that the pair bond in this case is not a typical covalent Si-Si bond like that shown for a bulk bond in fig. 4a, as seen both from the overlap geometry and the relatively high s-content still present, and at 2.75 a is also much longer than the bulk length. However, for modelling the surface we must also take into account that mutual
588
W.S. Verwoerd / Cluster calculations
b
I
Fig. 4. Localised orbitals (LO’s) calculated for the single pair cluster: (a) unrelaxed, (b) outwardly relaxed and (c) inwardly relaxed. Outlines are curves of constant wave function amplitude $ = 0.21 on each atom, a value arbitrarily chosen to give a good overlap for typical bulklike bonds like that shown in fig. 4a. Each LO has a total occupation of two electrons and is represented by a different kind of shading. The percentage s content of each atomic hybrid is shown and a reduced scale has been used perpendicular to the hybrid direction to clarify the directionality.
W.S. Verwoerd / Cluster calculations
589
approach of any two surface atoms moves them away from their outer neighbours along the dimer axis, and according to the results discussed above this will generate a repulsive interaction with compressive effects on the pair. To study such effects, we therefore add another surface atom, arriving at the SirsHI cluster discussed in section 2. The first calculation performed with this cluster is to let the two outer surface atoms approach the centre one symmetrically, as expected for the hypothetical 3 X 1 pattern of fig. Id. We find that in the optimum configuration under this restriction, the energy gain per atom is considerably less than for pairing, there is little evidence for bond formation (see fig. Sa) and if the symmetry restriction is lifted the system is unstable against the formation of a pair while the remaining atom is repelled. We believe that this result is, also of relevance for a consideration of the chain model of 2 X 1 reconstruction mentioned in the introduction. The chain model may also be interpreted as a linking of the external dangling bonds of a pair to other atoms, but involves much larger distortions of the back bonds than the 3 X 1 pattern. Since we find the latter to be energetically unfavourable, it is difficult to believe that the chain model would be feasible. The final geometry attained when one surface atom is allowed to move freely while the other two are coupled to it by the restrictions of the symmetric dimer model, is shown in fig. 5b, and the corresponding electronic structure in fig. SC. Comparison with the case of a single pair shows the additional neighbour to have a substantialleffect. There is a large decrease in the dimer bond length, making it even shorter than the bulk bonds. The electronic structure now shows a welldeveloped bulk-like dimer bond, while the dangling bonds have become largely p-like and their orientation suggests that some n-bonding between them takes place. This may explain the short bond length, and evidence for such n-bonding in the symmetric dimer model has also been found by Appelbaum, Baraff and Hamann [7,14:] from their band structure calculations, leading them to a bond length proposal of 2.22 a very similar to our optimised value. The result above is certainly not conclusive about the dimer state on a surface, since only one neighbour is included and furthermore its electronic configuration is different from that of a pair member. Nevertheless, it does show that there can be a surprisingly large interaction between surface pairs despite their separation of about 5.6 A. Also shown on fig. Sb are some 2nd and 3rd layer relaxations. This was obtained by allowing the indicated atoms freedom of motion subject to the 2 X 1 symmetry of the first layer. The displacements that we find agree in direction with those found from elastic models by Appelbaum and Hamann [ 1 l] and Chadi [6], but especially our second layer movements are much smaller. Part of the difference may be ascribed to the fact that our cluster does not include the full nearest neighbour complement for subsurface atoms, and that it is too small to include 4th layer relaxation. However, we believe that it is also partly due to an overestimate of the
590
W.S. Verwoerd / Cluster calculations
back bond strain by elastic models because of the assumption that these are identical to bulk bonds. This point will be discussed further in section 4. The subsurface relaxation also plays a role in distinguishing between the 2 X 1 and 2 X 2 reconstruction patterns of figs. lb and lc. If we define a dimer plane (for each dimer pair) containing the pair bond and the [loo] axis, it is easily seen that the subsurface relaxations induced by pairs in neighbouring dimer planes will be mutually reinforcing for the 2 X 1 pattern but opposing for the 2 X 2 structure, so that the former will be energetically more favourable. Seen dynamically, the pairing of two atoms on the unreconstructed surface will propagate longitudinally
a Ay=O.832 AZ =-0.356
W.S. Verwoerd / Cluster calculations
591
Fig. 5. Localised orbitals and geometry of the relaxed triplet cluster: (a) 3 X 1 configuration, (b) and (c) symmetric 2 x 1. Lengths are measured in A and LO’s are drawn to the same conventions as in fig. 4.
(in t.he dimer plye) by the repulsive interaction discussed previously, and transversely (perpendicular to the dimer plane) by the fact that a neighbouring surface pair will tend to be pushed together by the relaxation that the first pair induces in the second layer.
4. Asymmetric dimerisation We now turn to a free optimisation of the single pair SigHI cluster, allowing both surface atoms to move independently. This results in an asymmetric dimer configuration, shown in fig. 6a. Whereas a symmetric dimer can be specified by the dimer bond length d and the overall surface layer relaxation 6, the asymmetry gives rise in addition to a lateral shift h of the centre of gravity of each pair, and a buckling B of the surface layer. Parameter values are listed in table 1 and will be discussed below. First, however, we consider the reconstruction mechanism. A comparison of the localised orbitals of the symmetric pair (figs. 4c and SC) and the asymmetric case (figs. 6b and 6c) shows that an electron has been transferred from one of the degenerate singly occupied dangling bond orbitals of the symmetric pair to the other, to give a doubly occupied orbital on the “raised” member of the asymmetric dimer. This suggests a Jahn?Teller mechanism. To substantiate it, one must confirm that the degeneracy was lifted by the geometric distortion, and this is done by noting
592
WS. Verwoerd / Clusrer calculations
b
I
Fig. 6. Asymmetric dimerisations: (a) geometric parameters characterising the reconstruction; (b) LO’s of a single pair, drawn as in fig. 4; Cc) LO’s of the triplet, dangling bonds on the outer atoms.
the large difference in hybridisation of the two dangling bonds in fig. 6c. Almost all of the electron occupation is carried by the s-like and therefore low-energy orbital while the p-like orbital is &most empty. The same applies to fig, 6b, but here the occupation of the p-like orbital is too small to be seen in the Lo’s, which by
W.S. Verwoerd / Cluster calculations
593
Table 1 The reconstructed surface layer atomic geometry, all entries in A; 6 represents an overall relaxation and h a lateral shift of the dimer pair centre of gravity, the buckling 6 is the deviation along [ 1001 from this mean position, and d is the dimer bond length Coord.
6
h
b
d
Ref.
AH Levine Chadi Pair 2x1 2x2
-0.09 -0.22 -0.20 -0.25 -0.31 -0.29
0.0 0.0 0.31 0.20 0.23 0.21
0.0 0.0 0.24 0.24 0.24 0.24
2.45 2.35 2.35 2.46 2.40 2.42
[l&13] [ll
161 Present work Present work Present work
definition are constructed from occupied orbitals. The change in dangling bond hybridisation on each atom is accompanied by an opposite change in its other bonds, as clearly seen in the dimer bond in fig. 6, and this is also true of the back bonds (see table 2). According to a well known argument that bonds formed by s-like hybrids tend to be shorter than those that are p-like [ 181 one would therefore expect differences in the back bond lengths, and indeed the “donor” atom has a back bond length of only 2.32 a compared to the bulk cluster value of 2.37 a and the “acceptor” atom’s 2.42 a. While we have previously found evidence [ 181 that the MIND0/3 program exaggerates the effect of hybridisation on bond lengths, it seems reasonable to conclude that the connection between geometrical distortions and the electron transfer that make up the Jahn-Teller effect, is brought about here by accompanying hybridisation and bond length changes in the dimer and back bonds. Table 2 also illustrates the penetration of surface effects on the electronic structure into the cluster. It is seen that even the considerable charge redistribution
Table 2 The hybridisation @OScontent) on each atom of bonds in the SisH 12 cluster Layer
Bulk
Unrelaxed
Symmetric
Asymmetric
-
-
_
-
70 70
79 3
45 47
8 51
54 47
1‘3 43
56 44
48 45
4.5 46
40 44
52 47
38 48
45 47
47 47
46 47
45 47
45 47
594
W.S. Verwoerd / Cluster calculations
associated with the asymmetric reconstruction stage, only seriously affects bonds in which the first layer atoms participate. This justifies the use of elastic models for deeper bonds as done in refs. [6] and [ 111, but we find that the differences in the surface bonds lead to smaller deep layer relaxations for the asymmetric dimer. The increase in dimer and back bond s content of the donor atom may be expected to increase its electronegativity [2 11, thus acquiring more electrons (while the reverse is true for the acceptor) and so counteracting the effect of the dangling bond charge transfer. This effect is visible in the dimer bond of fig. 6 which shows the s-like hybrid of the donor atom to be slightly larger than its partner, and is reflected in the fact that we find the nett charge transfer between the two atoms to be considerably smaller than a full electron, namely 0.221el. The latter figure is derived from the calculated nett atomic charges on the lowered and raised atoms, namely 0.32lel and -0.12lel respectively. Taking into account also the atomic charges on the unrelaxed surface (0.191el) and on the symmetric pair (-0.05lel) the picture that emerges is that 0.24 electrons per atom is transferred to the surface layer during formation of the dimer bond, while in the Jahn-Teller stage the donor loses 0.37 electrons of which -0.22lel is transferred to its pair partner, while 0.15 electrons is gained from each atom by the substrate. Our value for the pair charge transfer is considerably smaller than the 0.361el value of Chadi [6], as is to be expected from our inclusion of Coulombic effects. On the other hand, it is quite similar to the -0.20lel that we have found for Jahn-Teller reconstruction of the (111) surface [ 161. The geometrical results of our single pair configuration is shown in table 1 to be quite similar to those of Chadi, with an identical value predicted for the buckling. However, our calculation gives more “vertical” relaxation and less lateral displacement of the pair, implying bigger back bond changes. Also notable is that we find the dimer bond to be somewhat longer than the bulk bond length, which agrees with the value used in recent LEED analyses [ 1 l-l 31. The asymmetry of dimer pairs gives rise to various possible arrangements on the surface with different periodicities, even assuming the basic pairing to follow the 2 X 1 pattern. In the dimer plane (containing the dimer bond and the [loo] axis) it can be either two-fold or four-fold in terms of the unreconstructed mesh, and perpendicular to this plane have either single or double period. Some of the resulting configurations are shown in fig. 7. In order to assess which of these have the lowest energy, we next have to evaluate the interaction between dimer pairs. First, we consider the interaction between a dimer and its nearest neighbour in the direction perpendicular to the dimer plane. As described in section 2, we use a 15-atom cluster containing two such pairs for this purpose, and find that the reconstruction energy decrease of the double period configuration (crossed dimers) is at A,!? = -0.77 eV significantly more favourable than the single period (parallel dimer axes) case with AZ? = -0.70 eV (all energies compared with the unreconstructed surface). As a basis of comparison, it may be noted that for a single pair we find
W.S. Verwoerd /Cluster calculations
595
(b) 4 xl
(cl
2 x2
Fig. 7. Different periodicities resulting from the superposition the basic 2 X 1 dimerisation pattern.
AE = -0.61 eV. By comparison
of Jahn-Teller
reconstruction
on
of the symmetric and asymmetric pairs, we find that roughly -0.2 eV may be ascribed to the Jahr-Teller phase and -0.4 eV to the initial symmetric dimer bond formation. The discrepancy between the latter figure and the typical Si-Si bond energy of 2.0 eV [15] is mainly due to the angular distortion of the back bonds, as may be justified by noting that the strain energy calculated from the empirical bond force constants used by Harrison [3] amounts to about 1.7 eV. Secondly, the interaction of a dimer with its neighbours in the dimer plane has to be considered. Unfortunately the cluster size for a 4-layer double pair cluster becomes excessive in this case. However, inspection of fig. 6c shows that the dangling bond orbital of the third atom in the previously discussed triplet cluster is actually quite similar to that of the dimer acceptor atom. The interaction between doubly occupied dangling bonds was shown above to be substantial, and we try to model at least this part of the interaction between dimers by using the third atom of the triplet clusteras representative of a neighbouring dimer. It should be noted, however, that the nett atomic charge on this third atom is +0.2lel instead of
596
W.S. Verwoerd 1 Cluster calculations
-0.12lel (mainly because it lacks participation in a dimer bond) for which an electrostatic energy correction must be added to the triplet results. We find that the energy of the triplet with occupied dangling bonds on the outer atoms (as for periodicity of 2) is -0.97 eV compared to -0.95 eV when they are on adjacent atoms (periodicity of 4). From the energy of the outwardly relaxed pair (fig. 3b) we derive an energy contribution for the relaxation of the third atom only, and subtracting this and the electrostatic energy correction, we find reconstruction energies of -1.0 eV and -0.90 eV per dimer respectively for the two configurations. It should be noted that this does not include the electrostatic interaction between dimer pairs, which we find to increase the difference between the two configurations by an additional 0.09 eV. The point charge electrostatic interaction thus contributes roughly half of the total interaction energy between neighbouring dimers in this plane, and this is also found to apply perpendicularly to the dimer plane. Combining these results, one may now compare the dimer energies of the three reconstruction patterns in fig. 7 and of the 4 X 2 pattern formed by combining the periodicities in figs. 7b and 7c. Including the interaction of a dimer with its four nearest neighbour pairs (one on each side, both in the dimer plane and perpendicular to it) we find an energy of -1.25 eV for the 2 X 2 pattern of fig. 7c, -1 .l 1 eV for both the 2 X 1 and the 4 X 2 patterns and -0.97 eV for the 4 X 1 pattern. Structural parameters for the 2 X 1 and 2 X 2 arrangements (arrived at by interpolation of the triplet and double pair results) are shown in the last two columns of table 1. The actual energy values mentioned here are probably not very reliable, especially as a result of the rather indirect approach that we had to use for interactions between dimers in the same plane. However, we believe that the energy margin favouring the 2 X 2 structure is sufficient to qualify it as a serious contender for the interpretation of LEED and other experimental results. The various arrangements of asymmetric dimers on the surface have also been considered by Chadi [6]. While he mentions a variety of experimental indications in favour of a 2 X 2 pattern, his final conclusion is that the 4 X 2 structure has the lowest energy, mainly as the result of a difference in surface Madelung energy. Consequently we have also investigated the effect of including electrostatic interactions with the more remote dimer pairs, but find the differences between their contributions for different arrangements to be too small to affect our result in favour of the 2 X 2 pattern seriously. The discrepancy between our conclusion and that of Chadi, is probably due mainly to the smaller dimer dipole moment that we find, and to the fact that we find substantial nearest neighbour interactions which are not of a point charge electrostatic nature and which already favours the 2 X 2 structure. Another factor which might contribute, is the subsurface relaxation. Both from our own results and those of Chadi it appears that this contribution to the dimer relaxation energy is of the order of 0.1 eV, which is comparable to the energies separating the various patterns. However, our calculations indicate that this is
W.S. Verwoerd / Cluster calculations
597
related mainly to the dimer bond formation stage, whch is also supported by the very similar subsurface relaxations reported by AH [ 1 I] for the symmetric and Chadi [6] for the asymmetric dimers. Since the surface layer patterns under discussion only differ in the way the asymmetric part of the reconstruction is arranged, we therefore do not expect subsurface relaxation to play an important role in distinguishing between them.
Acknowledgements The author gratefully acknowledges financial support for computer time received from the Council for Scientific and Industrial Research, Pretoria. He is also indebted to the Alexander von Humboldt foundation for the receipt of a fellowship while the later stages of this work was completed.
[l] J.P. Levine, Surface Sci. 34 (1973) 90. [ 21 R.E. Schlier and H.E. Farnsworth, Semiconductor Surface Physics (Univ. of Pennsylvania Press, Philadelphia, PA, 1957). [3] W.A. Harrison, Surface Sci. 55 (1976) 1. [4] I’. Jona, H.D. Shih, A. Ignatiev, D.W. Jepsen and P.M. Marcus, J. Phys. Cl0 (1977) L67. [5] K.A.R. Mitchell, Surface Sci. 64 (1977) 797. [6] D.J. Chadi, Phys. Rev. Letters 43 (1979) 43. [7] J.A. Appelbaum, G.A. Baraff and D.R. Hamann, Phys. Rev. Letters 35 (1976) 729; Phys. Rev. B14 (1976) 588; B15 (1977) 2408. [8] W.A. Harrison, Surface Sci. 55 (1976) 1. [9] J.A. Appelbaum, G.A. Baraff and D.R. Hamann, Phys. Rev. Bll (1975) 3822. [lo] G.P. Kerker, S.G. Louie and M.L. Cohen, Phys. Rev. B17 (1978) 706. [ll] J.A. Appelbaum and D.R. Hamann, Surface Sci. 74 (1978) 21. [ 121 K.A.R. Mitchell and M.A. Van Hove, Surface Sci. 75 (1978) 147L. [13] S.V. Tong and A.L. Maldonado, Surface Sci. 78 (1978) 459. [14] J.A. Appelbaum, G.A. Baraff and D.R. Hamann, Phys. Rev. Letters 36 (1976) 450. [15] W.S. Verwoerd and F.J. Kok, Surface Sci. 80 (1979) 89. [ 161 El. Gaigher and W.S. Verwoerd, Surface Sci., submitted. [17] R.C. Bingham,M.J.S. Dewar and D.H. Lo, J. Am. Chem. Sot. 97 (1975) 1285; M.J.S. Dewar et al., QCPE 10 (1976) 309. [18] W.S. Verwoerd, Chem. Phys. 44 (1979) 151; W.S. Verwoerd, Intern. J. Quantum Chem., to be published. [ 191 Quantum Chemistry Program Exchange, Indiana University, Bloomington, IN 47401, USA. [20] D.J. Chadi, Phys. Rev. Letters 41 (1978) 1062. [21] J. Hinze and H.H. Jaff6, J. Am. Chem. Sot. 84 (1962) 540.