Structure and oscillatory multilayer relaxation of the bismuth (100) surface

We present a combined experimental and theoretical study of the surface structure of single crystal Bi(100) via scanning tunneling microscopy (STM), low-energy electron diffraction intensity versus energy (LEED-IV) analysis and density functional theory (DFT). We find that the surface is unreconstructed and shows an unusually large oscillatory multilayer relaxation down to the sixth layer. This unexpected behavior will be explained by a novel mechanism related to the deeply penetrating electronic surface states. STM reveals wide (100) terraces, which are separated by two-layer high steps in which the shorter of the two interlayer spacings is terminating this surface, consistent with the LEED structural analysis and DFT.


Introduction
The surfaces of bismuth have recently attracted considerable attention because of the remarkable difference between bulk and surface properties: while bulk Bi is a semimetal with a very low density of states (DOS) at the Fermi level, all the surfaces studied so far are good reduceddimensional metals due to the presence of electronic surface states [1]- [4]. Moreover, these states have been shown to be non-spin-degenerate due to the strong spin-orbit interaction (SOI) in Bi [5]. This, and the semimetallic bulk character, leads to interesting dynamic properties [6,7]. It was suggested that the SOI was an important ingredient for finding metallic surface states on Bi [8]. While symmetry requires the states to be degenerate at some k-points, the strong splitting gives rise to a steep dispersion between these points. Some types of surface reconstructions, which are often observed to remove metallic states on semiconductors, are not permitted on Bi because of the SOI [9] (for a recent review on Bi surfaces see [8]).
Moreover, there is a fundamental interest in metallic surface states on Bi due to the electronic similarity of Bi to the topological insulator Bi 1−x Sb x (x ≈ 0.1) [10]- [12]. The surface of a topological insulator can always be expected to support metallic states because of topological arguments. These states should be stable against small perturbations and correspond to the edge states in the quantum spin Hall effect [13]- [17]. Similar metallic states should hence be found on the surfaces of pure Bi. Indeed, the surface states on Bi 0.9 Sb 0.1 (111) [18,19] have been shown to be very similar to those observed on Bi(111) [2,5,9].
Bi(100) stands out from the other surfaces because some of the metallic surface states penetrate very deeply into the bulk crystal, down to about 20 layers [3]. For Bi(110), on the other hand, surface states are rather localized in the surface [20] and for Bi(111) they are strongly localized in the first two layers. The difference between Bi(100) and Bi(111) is particularly striking because of the similar geometric structure, as both are pseudo-cubic (111) surfaces.
A deep penetration of surface states into the bulk, as found on Bi(100), is significant for Bi nanostructures such as thin films [21]- [23], clusters [24,25], and nanowires [26]- [28] because the penetration depth of the states is of the same scale as the size of the entire nanostructure. It has, in fact, been realized that Bi nanostructures cannot be described as a semimetallic core with a metallic surface [25], even if the core has the bulk structure, but that deeply penetrating surface states can have a major impact on the behavior of the whole nano-object [28].
The surface structures of bismuth have been largely unexplored except for the pioneering qualitative work of Jona [29] and our recent work on Bi(111) [30] and (110) [31]. In the present paper, we report on a first quantitative determination of the Bi(100) surface structure by a complementary combination of scanning tunneling microscopy (STM), low-energy electron diffraction intensity versus energy (LEED-I V ) analysis and first-principles calculations. From the two possible surface terminations, the one with the shorter interlayer spacing is shown to prevail. We find an oscillatory interlayer relaxation penetrating deeply into the crystal by both LEED analysis and first-principles calculations. Experiment and theory give consistent results as to the magnitude of the relaxation.
Oscillatory relaxation is an important surface phenomenon first predicted by Landman et al [32] with an electrostatic model for simple metal surfaces. This phenomenon has been observed on fcc(110), bcc(100), hcp(1010) [33]- [36], and some high-index surfaces. An important concept in models predicting the multilayer relaxation is that the ions follow the forces set up by the Friedel oscillations of the itinerant electron charge density near the surface [32,36]. However, in the case of Bi surfaces, the physical cause for the multilayer relaxations must be different: the bulk carriers in Bi have a very low Fermi energy and effective mass, leading to small Fermi wave vectors [37,38]. The periodicity for screening oscillations can be expected to be of the order of 10 nm [39], irrelevant for an oscillation in the interlayer spacing. The simultaneous appearance of deeply penetrating surface states and multilayer relaxations suggests that in the case of Bi(100) the multilayer relaxations are related to the surface electronic structure, and not to the bulk electronic structure, presenting a novel mechanism. We show that the dangling bonds, occurring on the bulk-truncated (100) surface, are the driving force for long-range relaxations that transform the surface region. The resulting surface structure is even more similar to Bi(111) than the truncated bulk structure, on which no nearest-neighbor bonds are broken. Thus, instead of a surface-localized dangling bond, we observe on the (100) surface a surface state that is spread over several layers.
In the following section, we visualize the bulk truncated surface structure of Bi(100). Section 3 describes experimental and computational methods, and section 4 presents results and discussion. Finally, a summary and the main conclusions are given in section 5.

The bulk-truncated Bi(100) surface structure
As a typical group V semimetal, bismuth crystallizes in the rhombohedral (A7) structure, which can be viewed as two interpenetrating FCC sublattices slightly stretched in the [111] direction [29,40]. There are two atoms per unit cell and three nearest neighbors for each Bi atom. The bonding to these three neighbors is comparatively strong and can be viewed as pseudo-covalent. The Bi crystal is most naturally viewed as hexagonal bilayers, stacked in the (111) direction. In this picture, there are no pseudo-covalent bonds linking the bilayers and, indeed, the crystal is easily cleaved in this plane.
Here we view the bismuth crystal as a stack of quasi-hexagonal Bi(100) layers. The top view and side view of Bi(100) are depicted in figure 1. Note that the surface index (100) is given in the rhombohedral structure; it is (0111) in the hexagonal system and (111) in the pseudo-cubic system. In each layer, an atom has two nearest neighbors at the same distance as on the Bi(111) surface of 4.54 Å and four next-nearest neighbors at a slightly larger distance  of 4.72 Å. The three nearest neighbors of any Bi atom in the bulk structure, however, do not lie in the same quasi-hexagonal layer but they connect these layers. The connections to the nearest neighbors are indicated as 'bonds' in figure 1. The side view shows that there are two different possible interlayer spacings. The termination of the surface is not known, but it seems likely that the shorter interlayer spacing prevails, not least because it requires the breaking of only one covalent bond per unit cell instead of two. The broken bonds at the surface are also indicated as dashed lines in figure 1. It can be seen that the structure is close to an ABCABC stacking sequence of the quasi-hexagonal layers, except that the fourth layer atoms are not quite in registry with the first layer atoms. They are actually shifted by a distance of 0.57 Å. Therefore, the stacked layers have a registry that does not repeat itself or, in other words, the stacking sequence has an infinite repeat distance due to its incommensurate translation with the in-plane lattice parallel to the mirror plane. This results in a very low symmetry of the surface. The only symmetry element is a mirror plane as indicated in figure 1(a), in contrast to the higher symmetry of a threefold axis and three mirror planes in Bi(111). Consequently there are some important differences between Bi(100) and Bi(111): Bi(111) is the natural cleavage plane and its creation does not require breaking any pseudo-covalent bonds; and Bi(111) has a higher symmetry than Bi(100).
Due to the low symmetry, and similar to the case of Bi(110) [31], there are fewer constraints for a re-arrangement of near-surface atoms. It would be possible to shift the atoms of Bi(100) along the mirror line without breaking the translational symmetry parallel to the surface, retaining a (1 × 1) LEED pattern. Such displacements have been considered in the LEED analysis and first-principles calculations.

Surface morphology observation by STM
The surface atomic morphology was investigated using a SPECS 150 Aarhus STM. The Bi(100) sample was introduced into the ultrahigh vacuum (UHV) chamber with a base pressure P = 1 × 10 −10 mbar. To clean the Bi(100) surface, cycles of 0.5 keV argon ion sputtering (20 min) and 120 • C annealing (45 min) were carried out. STM images were acquired both at room temperature and low temperature achieved by liquid nitrogen cooling. Two Zener diodes in the STM cradle, one for the STM tip and the other for the sample stage, were used to control the temperature during cooling.

Atomic structure determination by LEED
The LEED experiment was performed in a UHV chamber with a base pressure of 2 × 10 −10 mbar. The cleaning procedure was the same as that used in the STM study. The cleanliness of the surface was monitored by Auger electron spectroscopy. No surface contamination was detected. The sample was carefully aligned to achieve normal incidence of the primary electron beam by matching the kinematically calculated pattern with the experimental pattern and by comparing symmetry-equivalent beams. Note that the latter procedure can only be used to a limited extent due to the low surface symmetry [31]. A series of LEED-I V spectra were recorded with a CCD camera at the sample temperature of −160 • C. The sample cooling was realized via a liquid nitrogen flow. The temperature was stable within ±3 • C during data acquisition.
A quantitative analysis of the experimental spectra was performed with the standard package of ATLMLEED by Barbieri and van Hove [41]. The package is capable of handling multiple angles of incidence and/or multiple surface structures present on the surface. We have taken advantage of these features to explore the possible coexistence of multiple domains and to account for possible errors in the alignment of normal incidence. The optical potential, V im , was taken to be constant at −4 eV. Optimized phase shifts were calculated based on the energydependent inner potential suggested by Rundgren [42,43] in the form of where a i (i = 1, . . . , 4) are parameters to be optimized and E is the kinetic energy of the incident electron. The Pendry reliability factor (R p ) [44] was used to measure the agreement between experimental and simulated I V curves. It is assumed that, due to random fluctuations, the residual value of R P at the minimum value R P min is R P = R P min × √ 8|V im |/ E, where E is 6 the total energy range used in the intensity analysis. From this number and the curvature of R P near the minimum, we can estimate the uncertainty in the optimized parameters.

Surface structure by first-principles calculations
We have performed first-principles calculations of the surface crystal structure of Bi(100). The full-potential linearized augmented plane wave method in film geometry [45,46] as implemented in the FLEUR code is used and the local density approximation [47] to the density functional theory (DFT) is employed. Spin-orbit coupling (SOC) is included in the self-consistent calculations [48]. Evaluation of the surface relaxation has been carried out for a symmetric 14-layer film, both with the inclusion of the SOC term and without. Force calculations have been performed for the first four layers without SOC, while relaxations have been carried out only for the first two interlayer spacings with the inclusion of SOC. In the latter evaluations, we keep the interlayer spacings d 34 and d 45 equal to those obtained from the force calculation without SOC. The geometry is chosen such that on both sides of the film the first interlayer distance is a short one (see figure 1(b)). A wavefunction cutoff of 3.8 au −1 is chosen and the Brillouin zone is sampled with 32 k-points. For the calculation of the DOS and the charge density of the surface state, we used a 22-layer film with H termination on one side, similar to the calculations in [3], but with the relaxations as reported below. For the (111) surface a 20-layer film with one-sided H termination was chosen.

Surface morphology observation by STM
An atomically resolved room temperature STM image of Bi(100) is shown in figure 2(a). The small rhombus denotes the primitive unit cell. No indication of a surface reconstruction is observed. Figure 2(b) shows a typical STM image of the Bi(100) step morphology at room temperature. Wide terraces are observed with a width of a few hundred Ångstroms. Figure 2(c) depicts the height profile along the line shown in figure 2(b), showing the step structure. We have measured 634 random step positions of Bi(100) to get a histogram plot of the step heights shown in figure 2(d). A Gaussian fit (dotted curve) is applied to find the best estimated value of 3.76 Å, which is very close to the double interlayer spacing of 3.72 Å in the bulk. One interesting observation in figure 2(b) is the hump-like structure running along the left half of the image passing point B. It is believed to be a twin crystal structure similar to the twin microlayers observed on the Bi(111) [49] and Sb(111) surfaces [50]. The STM experiment has shown that only double-layer steps exist on the surface. These results indicate that all the terraces are geometrically equal and can be viewed as equivalent but incoherent domains in the LEED experiment. Although the above results suggest that single layer steps do not occur on Bi(100), it is nevertheless possible to prepare such metastable steps by slightly sputtering the cooled surface with 300 eV Ar + ions. Taking STM images immediately after this preparation reveals the presence of some single-layer steps; see figure 3(a). A height profile corresponding to the line in (a) is plotted in (b). It shows the following step height sequence: a double layer, a short single layer, a long single layer and two double layers, which are indicated by D, S, L, D and D in figure 3. A quantitative measurement gives the short and long single-layer spacing as

Dynamical LEED intensity analysis
The (1 × 1) LEED pattern for the surface was previously observed [3,29] but no quantitative LEED study has yet been reported. In our LEED-I V study we analyzed 16 symmetry inequivalent beams with a cumulative energy range of 2902 eV. The real part of the inner potential was obtained as V 0 = max{−8.57, 0.20 − 54.08/ √ E + 15.36} (eV). A rigid shift ( V 0 ) of this optimized energy-dependent inner potential was implemented and the optimum value of (−1.1 ± 1.5) eV was obtained. The small value of V 0 proves that V 0 obtained above is a very good model for Bi(100). Our results show a favorable short termination, which yields a best-fit between experimental and calculated data with R P = 0.27. The long spacing termination has been tested and ruled out by a worse level of agreement (R P = 0.45). The statistical errors for all parameters are evaluated based on the minimum R factor deviation of R P = 0.03.
The vertical interlayer spacings and lateral displacements along the mirror plane were optimized. Large oscillatory relaxations penetrate deeply into the sixth layer, reminiscent of the deeply penetrating surface electronic states [3]. The relaxation magnitudes are summarized in table 1. The first, third and fifth interlayers contract by (−10.2 ± 1.7), (−8.6 ± 2.0) and (−11.9 ± 2.4)%, respectively. The second and fourth interlayers expand by (+15.0 ± 1.5)% and (+15.3 ± 2.2)%, respectively. The sixth interlayer spacing is identical to the bulk value within the precision of the analysis, (−1.6 ± 2.1)%. The R P factors versus the deviation from the optimized values for each of the first six interlayer spacings, while keeping the other five fixed at the optimized values, are plotted in figure 4; they show a parabolic shape and the uncertainty values range from about ±0.03 to ±0.05 Å.
The atomic in-plane displacements in the direction of the mirror plane are (0.05 ± 0.09), (−0.32 ± 0.11), (0.09 ± 0.13), (−0.10 ± 0.15), (0.26 ± 0.19) and (0.09 ± 0.40) Å for the six topmost layers, respectively. Positive (negative) values mean displacements from the left (right) to the right (left) in figure 1(a). The alternating signs of the displacement from layer to layer   can be understood by considering the surface bonding configuration. When one dangling bond is broken for each atom, the non-vertical back-bonds between the first and second layers pull the first layer atoms in both the vertical and lateral directions. This effect can penetrate into the deeper layers and leads to an alternating sign of the displacements.
Some other non-structural parameters were also tested or optimized. The layer-resolved surface atomic Debye temperatures were found to be close to the bulk value of 119 K. So for all layers this value was used. Normal incidence has been confirmed by a variation of the polar incidence angle.
A comparison between experimental and calculated I V curves for the optimized parameters is displayed in figure 5. Almost all of the peak positions match and the overall agreement is good. The final R p factor of 0.27 for this system is acceptable when considering the low symmetry of Bi(100) and the complex pseudo-covalent structure.

Comparison with first-principles calculations
The energy-minimized geometric structure has been calculated for this surface from first principles. The calculations performed for bulk Bi without the inclusion of SOI give bulk short and long interlayer spacings of 1.763 and 1.925 Å, respectively. Evaluations that include the SOC term lead to a very slight modification of approximately 0.01 Å of these results. Our scalar relativistic force calculations give the first four interlayer spacing relaxations at 0 K. The results agree reasonably well with those obtained by the LEED-I V analysis at −160 • C (see table 1), considering that also the finite thickness of the film introduces constraints, especially on the deeper layer relaxations. Both the experiment and calculation lead to the oscillatory relaxation near the surface. The experimental results show a strong relaxation even in the deeper fourth and fifth layers. The relaxation of the fifth interlayer distance was not considered in the first-principles calculations, since this would require a thicker film for the simulation of the unperturbed bulk region. The effect of SOI on the relaxations is also small here and does not change the results by more than 2%. Since this is smaller than the expected effects of the relaxations of the deeper layers, we will not consider the modifications due to SOI further.
In our force calculations, we also optimized the position of the surface atoms in a plane parallel to the surface. By symmetry, this movement is then confined to the mirror plane shown in figure 1. We note that the displacements of +0.08, −0.09, +0.05, −0.08 Å for the first four layers are consistent with the experimental findings in their alternating signs. These numbers are approximately half as large as the interlayer relaxations that amount to −0.14, +0.31, −0.11 and +0.18 Å for the uppermost layers. To visualize the effect of these relaxations, we plotted the upper layers of the Bi(100) surface in figure 6, connecting all Bi atoms that are closer than 3.26 Å (about two times the covalent radius). It can be seen that the atoms of the layers S and S-1, as well as S-2 and S-3, are now connected by bonds that are very similar to the double layers found for the Bi(111) surface. Of course, while the (111) unit cell is hexagonal, the (100) surface unit cell is quasi-hexagonal and the interatomic distances in the double layer are larger than on the (111) surface. But as compared to the bulk of the crystal, where the double layers are oriented almost perpendicular to the surface, the orientation of the nearest-neighbor bonds has changed on the surface. The existence of double-layers is also suggested by the STM experiments, as shown in figure 2, where the stability of these entities can be seen from the step height of the line scans.
In a bulk-truncated structure, i.e. without relaxations, a dangling bond would be formed when the surface cuts through these double layers. Relaxations help to avoid this situation, by a partial reorientation of the preferred bonding plane. The surface state on this (100) surface is now distributed over several layers, accompanying the relaxations. In this way, the relaxation and the deep penetration of the surface state that was already observed in [3] are coupled. This can also be seen from a comparison with the Bi(111) surface and the DOS for both the (100) and the (111) orientation shown in figure 6. The (111) surface, where no nearest-neighbor bond is broken, shows just a weakly pronounced surface state at the Fermi level, which can be identified by a small shoulder in the layer-resolved density of states (LDOS). This shoulder was also experimentally verified by scanning tunneling spectroscopy of (111)-oriented Bi films [51]. In contrast, the Bi(100) LDOS shows a prominent peak at the Fermi level, which is still visible Figure 6. Comparison between the relaxed Bi(100) surface (left of both panels) and the Bi(111) surface (right of both panels). The left panels show the charge density resulting from states ±0.2 eV from the Fermi level for the first 13 layers. The structure is indicated by the Bi atoms, which are connected by lines if they are closer than 3.26 Å. The color coding for the charge density is the same in both plots. The right panels show the DOS for the topmost layers and two deeper-lying ones. Also from these data, a prominent surface state at the Fermi level can be observed in the upper layers of the (100) surface. several layers deeper in the bulk. Charge density plots from states ±0.2 eV from the Fermi level confirm the different decay of these states in the two surfaces.
The attribution of the surface relaxation to the existing surface states can be clearly made by considering the long Fermi wavelength of the bulk states, which essentially rules out any influence of these states on the short-range relaxations. For example, the projected Fermi wave vector perpendicular to the surface of the electron pocket around L is 0.03 Å −1 , which corresponds to a wavelength of about 100 Å, far longer than our observed oscillation period. A strong influence of the surface states on the relaxation has also been found and discussed for Al(100) [52]. One fundamental difference is that for Al(100) the projected bulk Fermi wavelength is similar to the spacing between the atomic layers, while for Bi(100) the surface relaxation is clearly dominated by the deeply penetrating surface states, due to small values of the Fermi wave vector for the bulk states.

Summary and conclusions
The structure of Bi(100) was determined by STM, quantitative LEED-I V analysis and firstprinciples calculations. The combination of these three approaches gives a consistent picture of a surface which is terminated by short interlayer spacing and a few hundred Ångstromwide terraces separated by double-layer high steps. Deeply penetrating oscillatory multilayer relaxations are observed that agree well with theoretical predictions.
Bismuth is a material on the borderline between covalent bonding, common for semiconductors, and metallic behavior. Whereas in the former class of materials the formation of a surface can result in dangling bonds that lead to extended reconstructions, in metals unreconstructed surfaces with surface states (usually accompanied by moderate relaxations) are often found. For the Bi(100) surface, we observe a deep impact of the surface on the electronic and geometric structures without significant changes of the in-plane periodicity. We find strong relaxations and a prominent surface state that might be characteristic of this class of materials.