Structure of the ( 111 ) surface of bismuth : LEED analysis and first-principles calculations

H. Mönig,1 J. Sun,2 Yu. M. Koroteev,3,4 G. Bihlmayer,5 J. Wells,1 E. V. Chulkov,3,6 K. Pohl,2 and Ph. Hofmann1,* 1Institute for Storage Ring Facilities, University of Aarhus, 8000 Aarhus C, Denmark 2Department of Physics and Materials Science Program, University of New Hampshire, Durham, New Hampshire 03824, USA 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, Apdo. 1072, 20080 San Sebastián, Basque Country, Spain !Received 6 February 2005; published 3 August 2005"


I. INTRODUCTION
Bismuth is a semimetal with a very small overlap between valence and conduction band. The density of states at the Fermi level is therefore small and the transport properties are quite different from those of normal metals. 1 Bi crystallizes in the rhombohedral A7 structure which is characteristic for the group V semimetals. In bulk Bi, each atom has three equidistant nearest-neighbor atoms and three equidistant next-nearest neighbors slightly further away. This results in puckered bilayers of atoms perpendicular to the ͓111͔ direction in which each atom is covalently bonded to its three nearest neighbors. The atoms' next-nearest neighbors are in the adjacent bilayers and the bonding within each bilayer is much stronger than the interbilayer bonding. This explains why Bi crystals easily cleave along the ͑111͒ plane. The A7 structure has two atoms per bulk unit cell, corresponding to the two atoms in the bilayers.
The surfaces of Bi have recently attracted considerable attention because their properties are very different from those of the bulk. Notably, the Bi͑110͒, 2 Bi͑100͒, 3 and Bi͑111͒ 4,5 surfaces are much better metals than the bulk due to the presence of metallic surface states. A particularly interesting feature of the surface electronic structure is the role of spin-orbit splitting. The spin-orbit interaction in the heavy element Bi is very strong. In the bulk band structure, the spin-orbit interaction is important 6 but it does not lead to a splitting of the bands with respect to the spin ͑i.e., every band still contains two electrons with opposite spin per k-point͒ because of the inversion symmetry of the bulk structure. On the surfaces, however, no inversion symmetry is present and the surface state bands are split with respect to their spin direction. 2,7 So far, no direct measurements of the spin-dependent surface band structure have been reported and the evidence for the spin-orbit splitting is the excellent agreement between the measured surface state dispersion and that calculated from first-principles with the inclusion of the spin-orbit interaction. 7 Additional direct evidence is found in the quasiparticle interference pattern on Bi͑110͒. 8 A quite different picture of the surface electronic structure of Bi͑111͒ was recently presented by Ast and Höchst. 9 They have shown that very good agreement between the measured surface state dispersion and a tight-binding calculation can be obtained for a very weak spin-orbit coupling parameter. However, in order to achieve this, a strong outward relaxation of the first interlayer spacing had to be assumed. This relaxation is ⌬d 12 / d 12 b = 71% where ⌬d 12 ͑=1.13 Å͒ is the change in the first to second interlayer distance and d 12 b ͑=1.59 Å͒ is the truncated bulk value for this distance. Such a big interlayer relaxation would be highly unusual. Moreover, the first interlayer spacing, d 12 , would be larger than the interbilayer distance d 23 , leading to a completely different type of bonding at the surface. In the first-principles calculations presented in Ref. 7 and here, the situation is quite different: they also give very good agreement with the experimental surface band structure but the spin-orbit splitting of the bands is pronounced and the interlayer relaxations are small. A direct structure determination of Bi͑111͒ is therefore required.
The amount of structural information on the Bi surfaces in general and Bi͑111͒ in particular is very limited. In a pioneering work, Jona has analyzed the LEED patterns for different Sb and Bi surfaces. 10 For Bi͑111͒ he found a hexagonal ͑1 ϫ 1͒ pattern and assumed that the surface was terminated with a short first interlayer spacing, i.e., with an intact bilayer. A later scanning tunneling microscopy investigation by Edelman et al. 11 supported the presence of only one termination. However, the detailed surface atomic configuration is so far unknown.
Here we present an experimental study of the surface structure of clean Bi͑111͒ by LEED intensity-vs-voltage ͑IV͒ analysis. Experimental diffraction intensities taken at sample temperatures of 140, 171, 218, 268, and 313 K at normal incidence have been analyzed by comparison to dynamical LEED calculations. The optimization of the structural and nonstructural parameters was carried out by minimizing a reliability ͑R͒ factor with the quadratic tensor model algorithm of Hanson and Krogh. 12 Two different surface models were tried: a termination with an intact bilayer and a termination with a broken bilayer. The agreement between measured and calculated intensities is excellent for the intact bilayer. The most important results are the structural parameters, i.e., the layer relaxations, and the mean-square vibrational amplitude of the surface atoms. In addition to the LEED study, we have performed first-principles calculations for the atomic structure of Bi͑111͒. The results of these are in good agreement with the experimental relaxations extrapolated to zero temperature.
The paper is organized as follows. Following this Introduction, the truncated bulk crystal structure for Bi͑111͒ is presented in Sec. II. Experimental and computational details are given in Sec. III followed by the results in Sec. IV. The paper concludes with a summarizing discussion in Sec. V.

II. Bi(111) CRYSTAL STRUCTURE
The truncated bulk crystal structure of Bi͑111͒ is shown in Fig. 1. For the figure, we have assumed that the crystal is terminated with a Bi bilayer. By considering only the ͑1 ϫ 1͒ LEED diffraction pattern reported by Jona, 10 a different termination would also be possible: the surface could be terminated with a long first interlayer spacing, i.e., with a bilayer cut upon surface creation. This is utterly inconceivable because of the covalent bonds that hold the bilayer together. We have indicated the locations of these bonds in the figure. To create a surface terminated with an intact bilayer, no covalent bonds have to be broken. The formation of a split bilayer, on the other hand, requires the breaking of three covalent bonds per surface unit cell. This picture is confirmed by the overall mechanical properties of Bi which is brittle and easily cleaves along the ͑111͒ plane, presumably between two bilayers. The bilayer-type structure gives rise to alternating interlayer distances. For the truncated bulk at 140 K we have 347 Å, and so on.
The rhombohedral A7 structure with two atoms per unit cell can also be described as a hexagonal structure with six atoms per unit cell. The hexagonal c axis lies in the rhombohedral ͓111͔ direction. In the hexagonal system, the lattice constants for 78 and 298 K are given in Table I. 13 Based on these lattice constants, we get the linear thermal expansion coefficients parallel and perpendicular to the ͑111͒ plane as 11.0ϫ 10 −6 K −1 and 18.5ϫ 10 −6 K −1 , respectively, which are consistent with thermal expansion coefficients reported by Erfling. 14 We obtain the lattice constants for 140, 171, 218, 268, and 313 K by interpolation.

A. Experiment
The experiment was performed in a -metal ultrahigh vacuum chamber with a base pressure of 7 ϫ 10 −9 Pa. It is equipped with a four-grid LEED optics from Omicron. Surface contamination was measured by Auger electron spectroscopy ͑AES͒ using a hemispherical electron analyzer and the LEED electron gun as an 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 sample temperature was adjusted via liquid nitrogen flow.
The surface was cleaned by cycles of Ar + sputtering ͑1 kV,2 A͒ 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 chargecoupled device ͑CCD͒ camera. A back-illuminated and Peltier cooled ͑−40°C͒ CCD chip guaranteed an extraordinarily high quantum efficiency. The camera was mounted on a base, which allowed rotation around all three axes. Great care was taken to align the camera with respect to the electron gun and the Bi crystal.
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 350 eV, while the energy was increased in steps of 1 eV after every recorded image. The integrated spot intensity of every single diffracted beam hk was extracted from these images. Normal incidence was found by minimizing the R factor between symmetry-equivalent beams. In the final data set, the intensities of the symmetry-equivalent beams were averaged.

B. LEED calculations and R factor analysis
The dynamical LEED calculations were performed using the computer code of Adams 15,16 which was developed from the programs of Pendry 17 and of Van Hove and Tong. 18 Atomic scattering phase shifts have been calculated using a muffin-tin potential model and the standard Barbieri/Van Hove phase shift package available online. 19 A muffin-tin radius of 2.87 atomic units has been adopted for Bi. Variations in the muffin-tin radius have only resulted in nonsignificant changes of the optimum R factor and structural parameters. Furthermore, the phase shifts have been renormalized by the thermal effects of root-mean-square ͑rms͒ isotropic vibrational amplitudes. Up to 17 phase shifts have been used because of the strong scattering of the heavy Bi atom ͑Z =83͒, especially at energies above 300 eV. The muffin-tin constant V 0 and the imaginary damping potential V im were both taken to be energy-independent. The surface potential step of height V 0 was located half a short bulk interlayer spacing away from the topmost layer. The reflection and transmission matrices for composite layers were calculated using Van Hove and Tong's combined space method. 18 The reflection matrices for the bulk crystal were calculated by Pendry's layer-doubling method. 17 The agreement between experimental and calculated LEED intensities is quantified by an R factor defined [20][21][22] as where I hk,i ex , I hk,i cal are the experimental and calculated intensities, respectively. The index i represents the data points spanning over the electron energy range. The experimental uncertainty of the beam hk, hk , was obtained by comparing symmetry-equivalent beams. 22 The global scaling constant c used for all beams can be determined by the requirement of ‫ץ‬R / ‫ץ‬c =0 as

͑2͒
In order to realize simultaneous optimization of all the structural and nonstructural parameters, Adams introduced the quadratic tensor model ͑QTM͒ algorithm for the R factor minimization; see Refs. 12 and 23.
Estimates of the parameter uncertainties i have been taken to be were n is the number of optimized parameters and R min is the local minimum value. The error matrix ⑀ can be obtained in terms of curvature matrix ␣ via ⑀ = ␣ −1 , where 24 N free is taken to be where ⌬k Ќ,hk ͑E͒ is the surface-normal component of the diffraction vector for the hk beam at energy E, E 2 − E 1 is the energy range of the measurements for the hk beam, and d is the long bulk interlayer spacing.

C. Ab initio calculations
We have also performed ab initio calculations of the surface crystal structure of Bi͑111͒. The full-potential linearized augmented plane wave method in film geometry 25,26 as implemented in the FLEUR code was used and the local density approximation 27 to the density functional theory was employed. Spin-orbit coupling was included in the selfconsistent calculations. 28 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 interlayer spacings without spin-orbit coupling while the total energy computations have been carried out for the first two interlayer spacings with the inclusion of SOC. As the latter evaluations are computationally more demanding we kept the d 34 and d 45 interlayer spacings equal to those obtained from the force calculation. The geometry was chosen such that both sides of the film were terminated with an intact bilayer. A wave-function cutoff of 3.6 a.u. −1 was chosen and the irreducible part of the Brillouin zone was sampled with 21 k-points.

A. LEED structural determination
After some initial tests with both structural models for the truncated bulk, the model involving a cut bilayer, i.e., a long first interlayer spacing, was discarded because of the poor agreement with the experimental data. The structural and nonstructural parameters for the surface terminated with an intact bilayer were refined using the procedure outlined above. The optimized parameters for a sample temperature of 140 K are listed in Table II. These are the interlayer spac- ings d ij , the rms atomic vibrational amplitudes u in the surface layers, the real part of inner potential V 0 , the damping potential V im , and the global scaling constant c. The Debye temperatures are calculated according to the relation between Debye temperature ⌰ D and mean-square atomic vibrational displacement ͗u 2 ͘ T at sample temperature T ͑Ref. 17͒ where m a is the atomic mass, ប Planck's constant, and k B the Boltzmann constant. At a sample temperature of 140 K ten symmetry in-equivalent beams with a total energy range of 2079 eV have been used to minimize the R factor. Experimental and calculated intensity-energy curves are shown in Fig. 2. The excellent agreement between the experimental and calculated results is evident. For comparison, the Pendry R factor for the agreement between calculated and experimental spectra is found to be R p = 0.25. It should be emphasized that a proper comparison requires a full and independent optimization of R p , which was not performed in the present work. The sensitivity of the R factor to the optimization parameters of interlayer spacings and atomic vibrational amplitudes has been investigated in the case of T = 140 K. Figure 3 shows the R factor as a function of these parameters individually. All the sensitivity curves take on a parabolic shape as discussed in Refs. 16 and 29. As expected, we can see a larger sensitivity of the R factor to the interlayer spacings than to the atomic vibrational amplitudes. From the error analysis outlined above, we determine that one standard deviation corresponds to a 4% increase of the R factor with respect to its minimum, i.e., R = 1.04R min . Hence the intersec-tions between R = 1.04R min and the sensitivity curves correspond to the error bars of the individual parameters.
The optimized parameters for the complete set of temperatures are given in Table III together with the number of beams and the total energy range. The first interlayer spacing shows no significant relaxation in the temperature range studied while the second interlayer spacing shows a small expansion.
The changes of rms atomic vibrational amplitudes with temperature for the first two layers are shown in Fig. 4. The corresponding Debye temperatures as calculated from Eq. ͑6͒ are also given. It is found that the first two layers exhibit lower Debye temperatures than the bulk for all sample temperatures investigated. Note that the Debye temperature for the first two layers increases slightly with sample temperature. There are two reasons for this. The first is the experimental uncertainty in the determination of u and the second is that the Debye model is only a very simple approximation.

B. Comparison to first-principles calculations
The first-principle calculations performed with the inclusion of spin-orbit interaction give bulk short and long interlayer spacings of 1.667 and 2.228 Å, respectively. Scalar relativistic evaluations that do not include the SOC term lead to a very slight modification ͑of approximately 1%͒ of these results. To check the sensitivity of the surface relaxation to the number of atomic layers incorporated into the relaxation, we have computed the 14-layer Bi͑111͒ film without the SOC term for two geometries. In the first geometry only the two surface layers were allowed to relax to the equilibrium and in the second geometry the four surface layers have been relaxed. Both these calculations result in a small contraction of the first interlayer spacing ͑⌬d 12 / d 12 b Ӎ −0.6% ͒ and in an expansion of the second one ͑⌬d 23 / d 23 b Ӎ + 6.6% ͒. The third interlayer spacing d 34 ͑d 34 b = d 12 b in bulk Bi͒ also experiences a small contraction similar to ⌬d 12 while the fourth one expands by ⌬d 45 / d 45 b Ӎ + 2.3%, i.e., significantly less than the expansion of the d 23 distance. The inclusion of spin-orbit interaction into the calculation results in a small expansion of the d 12 distance, ⌬d 12 / d 12 b = + 0.6%, and in an expansion of the second interlayer spacing, ⌬d 23 / d 23 b = + 6.2%. These data can be compared with the experimental results linearly fitted and extrapolated to 0 K which are ⌬d 12 / d 12 b = +͑1.2± 2.3͒% and ⌬d 23 / d 23 b = +͑2.6± 1.7͒%. For the ⌬d 12 / d 12 b relaxation a very good agreement is obtained with the computed value when the SOC term is included.
For ⌬d 23 / d 23 b the theoretical value clearly lies outside the experimental error but the agreement can still be judged to be fair in view of the remaining uncertainties. In Fig. 5 we show the calculated difference between total energies per surface atom for the films, as a function of the relaxation for the first and the second interlayer spacings as well as the whole energy landscape around the obtained energy minimum. The ⌬d 12 / d 12 b relaxation curve ͑open squares͒ was computed keeping ⌬d 23 / d 23 b = 0%, while the d 23 was relaxed  keeping ⌬d 12 / d 12 b =1% ͑open circles͒. As one can see in the figure, ⌬d 23 shows an extremely weak energy dependence. For instance, the energy difference between ⌬d 23 / d 23 b =2% and ⌬d 23 / d 23 b = 6% is about 4.3 meV/ surface atom. This very small quantity can be influenced by the choice of the approximation for the exchange-correlation potential, by the finite thickness of the slab, by small inconsistencies in the k-point sampling between bulk and film calculations, as well as by the fact that the linear extrapolation to 0 K may be oversimplified.

V. DISCUSSION AND SUMMARY
Our results give a consistent picture of the surface geometric structure of Bi͑111͒. The most important result is the fact that the relaxations are quite small, unlike the values assumed in the tight-binding fit to the electronic structure in Ref. 31. Indeed, moderate relaxations would be expected for a system like Bi͑111͒ that can be viewed, drastically simpli-fied, as a stack of covalent bilayers bound by van der Waalslike bonds. It is even conceivable that the relaxation in d 23 should be bigger than that in d 12 . The argument is again that the much stronger bonds within a bilayer would be modified less upon a change in the environment.
Another interesting result is the reduced surface Debye temperature for the first two atomic layers. This finding is consistent with a very early study of Goodman and Somorjai. 32 Reduced surface Debye temperatures are a common phenomenon on many surfaces. The actual numerical values of the surface Debye temperature are an important ingredient for the determination of the electron-phonon coupling strength from angle-resolved photoemission data. 33,34 Finally, there are some interesting technical issues worth mentioning. The QTM algorithm introduced for the R factor minimization by Adams has been tested 23 for the surface structure determination of clean Al͑110͒, Al͑100͒-c͑2 ϫ 2͒ − Li and Al͑111͒ − ͑2 ϫ 2͒ − Na. The results obtained by this new procedure showed an excellent agreement with those from a simple grid-search procedure. 24,35 Some advantages have been shown and discussed, including ͑i͒ the rapid convergence to the optimum values, ͑ii͒ the ability to define linear correlations between search parameters and perform parameter optimization within defined intervals, and ͑iii͒ the capability of optimizing both the structural and nonstructural parameters simultaneously. These characteristics have been confirmed by our calculations. Tests for some more complicated systems, such as surface alloy structures, are under way. An important feature in our optimization is the employment of a global scaling constant c for all beams, which makes the R factor retain the information of relative intensities from beam to beam.

ACKNOWLEDGMENTS
Valuable discussions with Anders Mikkelsen are gratefully acknowledged. This work was supported by the American National Science Foundation, Grant No. DMR-0134933, the Danish National Science Foundation, the Basque Country Government, and by the University of the Basque Country.