Structural determination of the Bi(110) semimetal surface by LEED analysis and ab initio calculations

J. Sun,1,* A. Mikkelsen,2,† M. Fuglsang Jensen,2 Y. M. Koroteev,3,4 G. Bihlmayer,5 E. V. Chulkov,3,6 D. L. Adams,7 Ph. Hofmann,2 and K. Pohl1 1Department of Physics and Materials Science Program, University of New Hampshire, Durham, New Hampshire 03824, USA 2Institute for Storage Ring Facilities and Interdisciplinary Nanoscience Center (iNano), University of Aarhus, 8000 Aarhus C, Denmark 3Donostia International Physics Center (DIPC), 20018 San Sebastián, Basque Country, Spain 4Institute of Strength Physics and Materials Science, Russian Academy of Sciences, 634021, Tomsk, Russia 5Institut für Festkörperforschung, Forschungszentrum Jülich, D-52425 Jülich, Germany 6Departamento de Física de Materiales and Centro Mixto CSIC-UPV/EHU, Facultad de Ciencias Químicas, UPV/EHU, Apartado 1072, 20080 San Sebastián, Basque Country, Spain 7Institute for Physics and Astronomy, University of Aarhus, 8000 Aarhus C, Denmark Received 16 June 2006; revised manuscript received 5 October 2006; published 6 December 2006


I. INTRODUCTION
Characteristic of group-V elements, bismuth crystallizes in the rhombohedral A7 structure as a semimetal with a small density of states at the Fermi level. 1 But interestingly, the surfaces of Bi show very different electronic properties from the bulk. Studies on the Bi͑110͒, 2 Bi͑100͒, 3 and Bi͑111͒ ͑Refs. 4 and 5͒ surfaces have shown that they are much more metallic than the bulk due to a significantly higher density of states at the Fermi level at the surface. 6 It has been found that one significant contribution is from a strong spin-orbit coupling at the surface due to broken inversion symmetry. 2,[6][7][8] From a chemical point of view, the creation of a surface requires the breaking of atomic bonds. Covalent bonding plays only a minor role in most metals. Thus the effect of bond breaking is small and surface properties are similar to those of the bulk, although localized electronic surface states may be present. On semiconductors, creating a surface leaves so-called dangling bonds which should give rise to half-filled and therefore metallic bands. However, it turns out that on most semiconductor surfaces the atoms rearrange their positions such that the dangling bonds are removed and the surface is again a semiconductor and not a metal. 9 Semimetals such as bismuth lie in between these two cases. On one hand, a semimetal is close to being a semiconductor since directional bonding is important and the valence and conduction bands are almost separated by a gap. On the other hand, there is a very small overlap between both bands such that the material is formally a metal. This delicate balance between being a metal and a semiconductor depends crucially on the atomic structure 10 and it can be expected to be severely disturbed at the surface.
Detailed structural information on Bi surfaces is so far limited to a recent low-energy electron diffraction ͑LEED͒ intensity vs voltage ͑IV͒ and first-principles study of the Bi͑111͒ surface. 11 One important difference between bulkterminated Bi͑110͒ and Bi͑111͒ is that the Bi͑110͒ surface exhibits dangling bonds, while Bi͑111͒ does not. In a pioneering study by Jona,12 oxygen adsorption experiments suggest that Bi͑110͒ is noticeably more active than Bi͑111͒. A qualitative analysis of LEED patterns in Jona's study shows an unreconstructed ͑1 ϫ 1͒ Bi͑110͒ surface structure. From the bulk structure Jona erroneously concluded that the unit cell ͑and hence the LEED pattern͒ should not be exactly rectangular but that the lattice vectors should include an angle slightly different from 90°. This is not correct, as will become apparent below. The unit cell is rectangular and almost square. A recent scanning tunneling microscopy study by Pascual et al., 8 revealed images of the Bi͑110͒ surface that are consistent with a near-square surface unit cell.
In contrast to most metal surfaces, Bi͑110͒ has a very low symmetry-the only symmetry element being a mirror plane while it has no translational symmetry normal to the surface. This makes the LEED IV analysis of this surface challenging. Along the ͓110͔ direction the Bi has a close stacking of atomic layers, i.e., the buckled bilayer as described below that requires the combined space method 13 to calculate diffraction matrices. 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 nonsymmetrical translation parallel to the mirror plane. More importantly, the dangling bonds at the surface may complicate the surface electronic and geometrical structures, which makes this open surface quite similar to many semiconductors and binary compounds.
In this paper, we report on a study of the surface structure of clean Bi͑110͒ by quantitative LEED intensity vs voltage analysis and ab initio calculations. Experimental diffraction intensities taken at a sample temperature of 110 K under normal incidence have been analyzed by comparison to dynamical LEED calculations. Great care was taken to align the sample considering the low-symmetry diffraction pattern. The main structural parameters that were optimized in the LEED IV analysis include the first four interlayer spacings and the Debye temperatures for the first two surface layers. Furthermore, we have performed first-principles calculations for the atomic structure of Bi͑110͒. The results are in good agreement with the experimental relaxations.
In the following section, we introduce the bulk truncated surface structure of Bi͑110͒. Experimental and computational details are described in Sec. III, followed by the results and discussion in Sec. IV. Conclusions are given in Sec. V.

II. BULK TRUNCATED Bi(110) SURFACE STRUCTURE
The A7 ͑␣-arsenic͒ structure of bulk bismuth has a rhombohedral unit cell with a two-atom basis. This structure can be obtained by a slight stretching of two penetrating facecentered cubic structures along the body diagonal. 14 The bulk truncated surface structure of Bi͑110͒ is shown in Fig. 1.
Each atom has three nearest neighbors to which it is connected by quasicovalent bonds. The side views show the stacked bilayers loosely bound by a single bond between every other atom in neighboring bilayers. Within one bilayer, each atom in one layer closely bonds with two nearestneighbor atoms in the other layer, forming a buckled structure. The covalent bonds have been drawn by solid lines and the dangling bonds at the surface layer by dashed lines. The bilayer-type structure gives rise to alternating interlayer distances. For the truncated bulk at 110 K we have d 12 b = 0.208 Å, d 23 b = 3.064 Å, d 34 b = 0.208 Å, d 45 b = 3.064 Å, and so on. Interlayer spacings between the ith and jth bulk layers are indicated as d ij b . Noticeably, the Bi͑110͒ surface has very low symmetry: the only symmetry element is a mirror plane as indicated in Fig. 1. The lengths of unit vectors at 110 K are taken as 4.731 and 4.538 Å; see Refs. 12, 14, and 15. If the rhombohedral structure is treated as a pseudocubic structure as in Ref. 12, Bi͑110͒ will be denoted as Bi͑100͒. The pseudosquare character of the surface unit cell is evident: for a cubic Bi structure all the atoms in the first bilayer would have the same height, the unit cell would be rotated by about 45°, and contains only one atom.

A. Experiment
The experiment was performed in a -metal ultrahighvacuum chamber equipped with a four-grid LEED optic with a base pressure of 7 ϫ 10 −9 Pa. Surface contamination was measured by Auger electron spectroscopy ͑AES͒ using a hemispherical electron analyzer and the LEED electron gun as electron source. The sample was mounted on a manipulator, allowing positioning to within 0.1°around all three axes of the crystal. The sample was cooled by liquid nitrogen. The surface was cleaned by cycles of 1 keV Ar + sputtering and annealing to 150°C. With AES no surface contamination could be detected. The maximum possible oxygen contamination was determined to be 0.02 monolayers. Spot intensities were measured using a 16-bit charge-coupled device ͑CCD͒ camera. A back-illuminated and Peltier-cooled ͑−40°C͒ CCD chip guaranteed a high quantum efficiency and low dark current. The camera was mounted on a base that allowed rotation around all three axes. Great care was taken to align the camera with respect to the electron gun and the Bi crystal, as described below.
To obtain intensities of the diffracted beams as a function of electron energy, the following procedure was employed. A series of images was recorded within the energy range from 30 to 300 eV, while the energy was increased in steps of 1 eV after every recorded image. The integrated spot intensity of every single diffracted beam ͑h , k͒ was extracted from these images. The presence of only one mirror line symmetry for Bi͑110͒ leads to technical challenges for the LEED experiment. These difficulties are illustrated in Fig. 2 that shows two measured LEED patterns taken at different incident energies. The pseudosquare pattern of the reciprocal lattice and the missing left ͑right͒ symmetry are clearly evident. The up ͑down͒ symmetry is given by the mirror plane in the crystal ͑the horizontal plane in Fig. 2͒. It is necessary to align FIG. 1. ͑Color online͒ Truncated bulk structure of Bi͑110͒. The solid and dotted lines mark covalent and dangling bonds, respectively. ͑a͒ Top view of the first two atomic layers. Each layer consists of a two-dimensional rectangular lattice and the lattice constants at 110 K are given. The mirror planes of the structure are also shown as dashed lines. ͑b͒ and ͑c͒ Side views of the first eight layers perpendicular and parallel to the mirror plane, respectively. The bilayerlike structure with alternating short and long interlayer spacings is evident. the sample surface perpendicular to the incoming electron beam, and this is usually done by comparing the IV curves of the symmetry-equivalent beams. Here this procedure can only be applied for the up ͑down͒ angle. In order to align the left ͑right͒ angle we optimized the diffraction spot position on the LEED screen until they agreed with the kinematically calculated positions. We estimate that this approach leads to an error of less than 1°in the angle of incidence. In the final data set, the intensities of the symmetry-equivalent beams were averaged.

B. Dynamical LEED calculations
The dynamical LEED intensity calculations were performed using the standard package SATLEED ͑symmetrized automated tensor LEED͒ by Barbieri and Van Hove 16 within the renormalized forward-scattering perturbation formalism. Atomic scattering phase shifts have been calculated using a muffin-tin potential model and the standard Barbieri and Van Hove phase shift package. 16 The bulk diffraction matrices for the closely spaced bilayers were calculated with the combined space method. 13 The same muffin-tin radius of 2.87 a.u. and phase shifts have been used as in Ref. 11. Phase shifts have been renormalized by the thermal effects of rootmean-square ͑rms͒ isotropic vibrational amplitudes. Up to 15 ͑L =14͒ phase shifts have been used because of the strong scattering of the heavy Bi atom ͑Z =83͒. The muffin-tin constant V 0 is taken to be energy independent and is optimized. V im , the imaginary part of the inner potential, also referred to as the damping or optical potential, is taken as 4 eV for the bulk and 4.2 eV for the first two overlayers. The slightly larger value at the surface was chosen to model the presence of dangling bonds, which increases the electron damping. The surface potential step of height V 0 is located half a long bulk interlayer spacing away from the topmost layer nuclei. The bulk Debye temperature is fixed at 119 K, 17 while the Debye temperatures for the first two layers are optimized. Mean-square atomic vibrational amplitudes ͗u 2 ͘ T at temperature T for the Debye-Waller factor calculation are derived from Debye temperatures ⌰ D according to the following equation: 18 where m a is the atomic mass, ប Planck's constant, and k B the Boltzmann constant.
In the LEED intensity analysis, agreement between experimental and calculated LEED intensities is quantified by the widely used Pendry R factor R P , which is particularly sensitive to relative peak position and the existence of small peaks. 19 The uncertainties in the optimized structural parameters are estimated from the variation around the minimum R P min , where ⌬E is the total energy range compared in the IV analysis. 19

C. Ab initio calculations
We have also performed ab initio calculations of the surface crystal structure of Bi͑110͒. The full-potential linearized augmented plane wave method in film geometry 20,21 as implemented in the FLEUR code was used and the local density approximation 22 to the density functional theory was employed. Spin-orbit coupling was included in the selfconsistent calculations. 23 The evaluation of the surface relaxation has been carried out for the symmetric 14-layer film, both with the inclusion of the spin-orbit coupling ͑SOC͒ term and without this term. Force calculations have been performed for the first four layers without spin-orbit coupling while relaxations have been carried out only for the first two interlayer spacings with the inclusion of SOC. In the latter evaluations we kept the interlayer spacings d 34 and d 45 equal to those obtained from the force calculation without SOC. The geometry was chosen such that both sides of the film were terminated with an intact bilayer. A wave-function cutoff of 3.8 a.u. −1 was chosen and the Brillouin zone was sampled with 32 k points.

A. LEED structure determination
The LEED pattern of Bi͑110͒ has previously been discussed by Jona. 12 He defined a pseudocubic bulk unit cell and concluded that the unit cell ͑and hence the LEED pat-tern͒ should not be exactly rectangular but that the lattice vectors should include an angle slightly different from 90°. Our study does not confirm this conclusion. Our LEED patterns as presented in Fig. 2 show an exact rectangular net from careful measurements of the diffraction spots positions and, indeed, such an exact rectangle can also be expected from a projection of the bulk reciprocal lattice onto the surface. 6 The measured ratio of the two reciprocal unit cell vectors is 0.96͑2͒ in good agreement with the expected value of 0.959. Moreover, the observed patterns show no indication of any reconstruction of the Bi͑110͒ surface, despite the existence of active dangling bond at the surface. Apparently, Bi͑110͒ is found to be very different from typical semiconductor surfaces, such as Si͑100͒ and Ge͑100͒ which both exhibit 2 ϫ 1 reconstructions.
The structural and nonstructural parameters were optimized for a Bi͑110͒ surface terminated by an intact bilayer. A termination with a split bilayer was immediately excluded due to lack of agreement with the experimental IV curves  spacings correspond to the small separation ͑0.21 Å͒ between the two layers making up the bilayer in the bulk. Their seemingly dramatic relative relaxations are very small in absolute terms. Also, the fourth layer appears to move above the third layer by 0.01 Å. However, this very small value is clearly below our detection limit. Overall no significant relaxation for the Bi͑110͒ surface is found. We have tried many possible displacement patterns allowed due to the low surface symmetry. However, we found no significant improvement in R P when changing the relative distance between the two basis atoms in the first and second layers parallel to the mirror line. The Debye temperature for the first layer is found to be lower than that of the bulk, which is consistent with an early study of Goodman and Somorjai. 24 Reduced surface Debye temperatures are a common phenomenon reflecting the weaker bonding of surface atoms compared to the bulk. 25 The actual numerical values of the surface Debye temperature are an important ingredient for the determination of the electron-phonon coupling strength from angleresolved photoemission data. [26][27][28][29] Meanwhile the second layer shows a Debye temperature close to the bulk value. The LEED IV analysis gives a relatively high R P factor of about 0.455 compared to typical values of 0.1 to 0.3 for clean unreconstructed metal surfaces. Many efforts have been made to find out the possible causes. We simulated non-normal-incidence conditions extensively in the LEED IV calculations and found that an increase in the incident angle gave a dramatic rise in the R factor from its minimum at zero or normal incidence. This suggests that the sample is properly aligned and the relatively high value is not caused by deviations from normal incidence during the IV measure-ment. The influence of the muffin-tin radius on the structure has been studied and results show essentially the same geometry with a minimum R factor at the muffin-tin radius of 2.87 a.u. We also tried the atomic potential derived from our ab initio calculations with no significance changes in the optimized structural parameters nor an improvement in the R factor. So it might be the structural complexity and low symmetry of the Bi͑110͒ surface itself that complicates the LEED process. As seen on open semiconductor surfaces, the presence of dangling bonds and the presence of voids in the open surface structure is a real challenge for the muffin-tin approximation of the crystal potential and could also contribute to the relatively high R P for this surface. However, the low surface symmetry of Bi͑110͒ gives rise to the large number of nonequivalent beams. Here we present an accumulative energy range of 3591 eV ͑20 beams͒ which is larger than about 1000 eV used in typical LEED studies. When an overall range of only 1071 eV ͑the first six beams͒ was analyzed in our work, the R factor decreased to 0.36 without changes in the optimized structural parameters. This value is comparable to LEED results for many semiconductors and metal oxide compounds with similar geometric and electronic structures. Interestingly, the reduction of R factor due to fewer beams indicates more intricate scatterings in the higher ordered beams. Furthermore, with horizontal surface atomic displacement in the mirror plane allowed, a smaller R has been obtained but, considering that the displacement values are within error limit and no significant geometry change occurred, we are not including these displacements in the report. The agreement between this large experimental data set and the calculated intensity-energy curves, as shown in Figs. 3 and 4, gives us confidence in the reliability of our results.
The error bars of the optimized parameters were analyzed based on the variation of the R factor around R P min , ⌬R P = 0.043 according to Eq. ͑2͒. The dependence of R P on a change of the interlayer spacings away from their optimized values is shown in Fig. 5. In this analysis, all other parameters were fixed at their optimized values. We can see that all the sensitivity curves take on a parabolic shape. The errors for the individual parameters are also listed in Table I.

B. Comparison to first-principles calculations
The first-principles calculations performed for bulk Bi without the inclusion of spin-orbit interaction give bulk short and long interlayer spacings of 0.142 and 3.087 Å, respectively. Evaluations that include the SOC term lead to a very slight modification of approximately 0.01 Å of these results. As shown in Table I These results agree reasonably well with those obtained by the LEED IV analysis at 110 K considering the fact that the absolute distance difference between the experimental and calculated first interlayer relaxations of −13% and −62% is only 0.06 Å. Both the experiment and theory lead to the contraction of the first interlayer spacing. For the second interlayer spacing the theory gives a small expansion while the experiment shows a small contraction of the spacing. However, the theoretical result is within the experimental error bar. The absolute distance difference between the experimental and calculated second interlayer relaxations of 0.015 Å is even smaller than that for the first interlayer spacing. For the third and fourth interlayer spacings the theory and experiment are in excellent agreement. The first-principles calculations that include the spin-orbit interaction term lead to ⌬d 12 / d 12 b = −43% and ⌬d 23 / d 23 b = + 0.4% for the first and second interlayer spacings, respectively. These values have been obtained by keeping the interlayer spacings ⌬d 34 / d 34 b and ⌬d 45 / d 45 b equal to those found in the scalar relativistic calculations. This shows that the influence of spin-orbit interaction on the relaxation is small and probably will not change the values of d 34 and d 45 significantly. Notice that in the relaxed geometry a change of 6% or 0.01 Å in ⌬d 12 / d 12 b corresponds to an energy change of only 0.5 meV per surface atom, which is certainly at the limit of our accuracy.
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 Fig. 1͑a͒. We notice that these relaxations are small and do not exceed 1.0% in the top four layers, consistent with the experimental findings.

V. CONCLUSIONS
Our results give a consistent picture of the very lowsymmetry surface geometric structure of Bi͑110͒ by LEED intensity analysis and first-principles calculations. Good agreement is reached between experimental LEED and theoretical IV curves. No structural reconstruction occurs despite of dangling bonds present at the surface. No significant absolute value of relaxation is found for the first four interlayer spacings. The reduced top-layer Debye temperature suggests essentially larger vibrational atomic amplitudes at the surface. Experimentally, the approach of sample alignment by calculating the diffraction spot positions on the LEED screen is very efficient and can be used for surfaces with low symmetry as well as for in-situ cleaved surfaces.