# Ab initio lattice dynamics and electron-phonon coupling of Bi(111)

## Abstract

We present a comprehensive ab initio study of structural, electronic, lattice dynamical and electron-phonon coupling properties of the Bi(111) surface within density functional perturbation theory. Relativistic corrections due to spin-orbit coupling are consistently taken into account. As calculations are carried out in a periodic slab geometry, special attention is given to the convergence with respect to the slab thickness. Although the electronic structure of Bi(111) thin films varies significantly with thickness, we found that the lattice dynamics of Bi(111) is quite robust and appears converged already for slabs as thin as 6 bilayers. Changes of interatomic couplings are confined mostly to the first two bilayers, resulting in super-bulk modes with frequencies higher than the optic bulk spectrum, and in an enhanced density of states at lower frequencies for atoms in the first bilayer. Electronic states of the surface band related to the outer part of the hole Fermi surfaces exhibit a moderate electron-phonon coupling of about 0.45, which is larger than the coupling constant of bulk Bi. States at the inner part of the hole surface as well as those forming the electron pocket close to the zone center show much increased couplings due to transitions into bulk projected states near . For these cases, the state dependent Eliashberg functions exhibit pronounced peaks at low energy and strongly deviate in shape from a Debye-like spectrum, indicating that an extraction of the coupling strength from measured electronic self-energies based on this simple model is likely to fail.

###### pacs:

73.20.At,68.35.Ja,63.20.kd,71.15.Mb## I Introduction

Looking back since the advent of modern solid state physics, bismuth (Bi) has been perhaps one of the most intriguing element for research. Its singular properties make Bi the metal with the lowest thermal conductivity after mercury and the element with largest diamagnetism. Bi allowed an early discovery of the Seebeck, the de Haas-van Alphen, the Shubnikov-de Haas, and the Nernst effects, all of which are inherently present in metals but were more challenging to observe in general. Bismuth was in fact the first metal whose Fermi surface was experimentally identified (2) and provided the basis to determine that of other metals. Moreover, bulk Bi was the first non-superconducting material (at least in the ordered phase) that was found to be superconducting in structurally bulk-like nanoparticles. (3)

Many of the outstanding properties of bulk Bi are linked to its peculiar electronic structure. Bi is a semimetal with tiny electron and hole pockets and a very small density of states (DOS) at the Fermi level. Metallic screening is therefore much weaker than in typical metals, and the interatomic bonding has a stronger directional component. This favors a layered crystal structure with alternating weak and strong interlayer bonds, which can be described as a stacking of bilayers.

Being a semimetal with low DOS at the Fermi energy, bulk Bi is expected to have a small average electron-phonon (-ph) coupling constant, . Simple estimates based on a tight-binding model(4) as well as a very recent ab initio calculation(5) yielded rather small values of =0.13 and 0.09, respectively, consistent with the finding that Bi does not display superconductivity down to 50 mK.(6) However, since long ago amorphous bulk Bi was found to have a quite strong electron-phonon coupling of =2.46 (the value is larger than that of amorphous lead, (7) which is known to have a strong electron-phonon coupling also in the crystalline form (8)). This indicates that the intrinsic coupling strength of electronic states to the phonons is comparable to those of lead, but superconductivity is only absent due to the very low DOS at .

Interest in surface properties of Bi arises on the one hand from the large difference of the surface electronic structure from the bulk one. Typically, an enhanced metalicity at the surface is observed which is related to the appearance of surface electronic states with an increased DOS at the Fermi energy ().(9) On the other hand, spin-orbit interaction is strong for Bi, and is a necessary ingredient for even a qualitative description of the electronic structure already for the bulk. On surfaces it invokes large splittings of surface states with deep implications for the Fermi surface topology. Among the low-index surfaces, Bi(111) is a kind of model surface for investigating these properties, because it does not break the bilayer structure. This enables the study of modifications in the electronic structure originating in the loss of translational symmetry in a pure form without complications due to additional breaking of chemical bonds. The surface electronic structure of Bi(111) has been elucidated in a series of photoemission experiments,(10); (11); (12); (13); (14); (15) which established the enhanced density of charge carriers and the presence of surface states split by spin-orbit interaction. These states create two types of Fermi surfaces, known as hole and electron pockets.

The enhanced DOS at Bi surfaces has stimulated experimental investigations of the coupling strength of surface electronic states. For the Bi(111) surface Ast and Höchst have estimated the electron-phonon coupling from spectral functions of angular-resolved photoemission spectroscopy (ARPES) measurements assuming a Debye model to account for the spectral shape.(13) The results turned out to be strongly dependent on the cutoff frequency of the assumed Debye spectrum. Fitting the data with a cutoff Debye frequency of 10 meV, appropriate for bulk Bi, yields =0.6, while assuming a surface Debye frequency of 5 meV yields =2.3. This strong difference was later ascribed mainly to the limited accuracy of the experiment.(16) An alternative road was pursued by Gayone et al., in which the imaginary part of the self-energy was extracted from the temperature dependence of the linewidth of momentum distribution curves.(17) They found smaller values of =0.40(5) for the surface electronic states, albeit still significantly larger than the bulk value.

A prerequisite for a detailed analysis of the electron-phonon coupling at surfaces is the knowledge of the surface vibrational spectrum. In particular, surface localized vibrational modes can couple more strongly to surface electronic states due to their larger spatial overlap.(18) In thin films, it has been shown that important contributions to the coupling may come from optical vibrations.(19) Experimentally the surface phonon spectrum of Bi(111) has been investigated very recently via inelastic Helium atom scattering (HAS) by Tamtögl et al.(20); (21) They found a rich spectrum of localized modes, including super-bulk modes with frequencies above the maximum bulk frequency.

In view of the experimental uncertainties in extracting reliable numbers for the electron-phonon coupling, theoretical calculations on an ab initio basis are highly desirable. The density functional perturbation theory (DFPT) provides a unified scheme to predict electronic, phonon and electron-phonon coupling properties. However, applications of this technique to bulk Bi and in particular to Bi surfaces are still a challenge, because the inclusion of spin-orbit interaction increases the computational effort substantially, especially for calculations of lattice dynamical quantities. Therefore, work along this line was first devoted to the electronic structure. Large splittings of surface bands due to spin-orbit interaction were found for all low-index Bi surfaces by Koroteev et al.(15) In a later work, it was shown that the surface electronic structure of Bi(111) converges very slowly with increasing thickness of the slab, in particular near the zone boundary at .(22)

Applications of DFPT to lattice dynamical properties including spin-orbit coupling (SOC) have been devoted first to bulk Bi.(23); (24) Díaz-Sánchez et al. showed that inclusion of SOC is crucial for an improved agreement between calculation and measurement,(24) but it did not provide such a good description as it was the case, for example, for the phonon dispersion of lead.(8). Very recently, DFPT has been used to investigate the lattice dynamics of Bi(111). Due to the numerical effort to include spin-orbit interaction, only slabs up to 6 bilayers could be treated,(21); (25); (5); (26); (27), raising the question to what extent the obtained results represent the surface dynamics in the limit of macroscopically thick slabs. For the 6-bilayer calculation, good agreement with the HAS data was found.(21); (27)

No attempt has been made so far to address the question of electron-phonon coupling of surface states within such an ab initio approach. The aim of this paper is to provide a comprehensive ab initio analysis of the electronic structure, lattice dynamics, and electron-phonon coupling of surface electronic states for the Bi(111) surface. Towards this goal, we will also discuss the bulk lattice dynamics, which enters the evaluation of the surface vibrations. As we want to focus on the properties of a semi-infinite surface, we carefully investigate the convergence of the various quantities with increasing thickness of the slabs.

The paper is organized as follows. The computational details are described in the following Section. Results are presented and discussed in Sec. III. First we briefly recapitulate the structure and lattice dynamics of bulk Bi in Sec. III.1 which is a prerequisite for the surface dynamics study, and present results for the electron-phonon coupling constant. Secs. III.2 to III.4 are devoted to the structural, electronic, and lattice dynamical properties of the Bi(111) surface, respectively, with emphasis on convergence with respect to the slab thickness. The coupling of surface localized electronic states is then discussed in Sec. III.5. Finally the results are summarized in Sec. IV.

## Ii Computational details

We performed density functional theory (DFT) calculations of bulk Bi and Bi(111) within the local density approximation (LDA) in the parameterization of Hedin and Lundqvist.(28) For the bulk geometry we also report our results obtained with the PBE variant of the generalized gradient approximation (GGA).(29) The electron-ion interaction was represented by norm-conserving pseudopotentials in the form proposed by Vanderbilt,(30) treating 6, 6, and 6 as valence states. The Kohn-Sham orbitals were expanded in a mixed basis (MB) of local functions and plane waves.(31); (32) Spin-orbit coupling was incorporated within the pseudopotential scheme via Kleinman’s formulation,(33) and used to obtain full charge-self-consistency.(8)

For convenience, we have chosen the rhombohedral representation of the structure of bulk Bi (which contains two atoms per unit cell) and the hexagonal one (9) to describe the Bi(111) surface. The Bi(111) surface was modeled with slabs having one atom per layer (11 in-plane periodicity) and thicknesses of 6 and 12 bilayers. In the calculations involving the surface, a vacuum of 12–14 Å separates the periodic images of the slab to avoid interaction between them. We have used the relaxed system to perform the lattice dynamics calculations.

Integrations over the bulk and surface Brillouin zones (BZ) were performed by sampling 121212 and 12121 meshes corresponding to 189 and 19 irreducible points for the bulk and surface calculations, respectively, combined with a Gaussian broadening with a smearing parameter of 0.1 eV. The Kohn-Sham orbitals were expanded in a basis set consisting of local and -type functions augmented by plane-waves with a kinetic energy cut-off of 12 Ry (163 eV). The Fourier expansion of the crystal potential and charge density has been found to be safely truncated at 50 Ry. The positions of all atoms in the slab were optimized until the forces on each atom and each direction was smaller than Ry (210 eV/Å). For this purpose, the Broyden-Fletcher-Goldfarb-Shanno algorithm (34) has been applied. The choice of the parameters has been verified through preliminary tests on the lattice parameters, band structure, volume, bulk modulus, and phonon frequencies at some high-symmetry points of bulk Bi. These tests also showed that the LDA is more appropriate than GGA.

The calculation of the lattice dynamical matrices at specific points of the surface BZ (SBZ) was performed using linear response theory embodied within DFPT. (35); (36) Its implementation in the MB scheme is described in Ref. (37). The dynamical matrices for bulk Bi and Bi(111) slabs were calculated at 189 and 7 irreducible points, respectively. Real-space force constants were obtained by taking the standard Fourier transform of the corresponding dynamical matrices.(38) The force constants calculated for the Bi(111) slabs were then combined with those of bulk Bi to model the dynamics of much thicker slabs.(39) The force constants for these thick slabs were used to derive the eigenvectors and frequencies of the vibrational eigenmodes of the slab at arbitrary points of the SBZ. Finally, electron-phonon coupling matrix elements were calculated directly from quantities obtained within DFPT including the contribution from spin-orbit interaction.(8)

## Iii Results and Discussion

Before discussing the properties of the Bi(111) surface, we first present our results for the structural, lattice dynamical and electron-phonon coupling properties of bulk Bi.

### iii.1 Bulk Bi

(Å) | (Å) | (Å) | |||
---|---|---|---|---|---|

LDA | 4.67 | 57.94 | 0.236 | 1.606 | 2.267 |

GGA | 4.93 | 56.18 | 0.232 | 1.621 | 2.513 |

LDA [(24)] | 4.69 | 57.57 | 0.234 | 1.576 | 2.326 |

Exp. [(41)] | 4.7236(5) | 57.35(1) | 0.23407(4) |

The A7 structure that characterizes Bi is a rhombohedral structure containing two atoms per primitive cell. In describing the crystal, one of them (R1) can be considered to lie at the origin and the other one (R2) along the trigonal axis (z-axis), as shown in Figure 1. The rhombohedral representation of the lattice is completely determined by three parameters, the length of its basis vectors, , and the two parameters gauging its departure from the simple cubic structure, (40) and . denotes the angle between two rhombohedral basis vectors, while is the internal parameter defining the position of R2 in the primitive unit cell. The crystal can also be viewed as being built from hexagonal planes stacked along the trigonal axis with two alternating distances, and . In terms of the rhombohedral parameters they are given by and , respectively, where is the length of the basis vectors projected onto the trigonal axis (see Fig. 1). Their ratio is solely determined by via . Our results for the lattice optimization are displayed in Table 1 together with those from measurement and previous calculations.

The A7 BZ is a distorted truncated octahedron (see Fig. 2). The distortions are relative to the trigonal direction (T). Six of the eight hexagonal facets are irregular and the squared facets are rather isosceles trapezoids. The regular hexagons contain the T-points and the irregular ones the L-points. It is customary to consider a binary-bisectrix-trigonal system as indicated in Fig. 2. The high-symmetry directions T , L, X, and the bisectrix direction lie in the mirror plane of the A7 structure, while the binary direction is perpendicular to it.

Figure 3 shows our calculated phonon dispersion curves along the bisectrix, trigonal, X, binary, and L directions together with the most comprehensive experimental data set by MacFarlane (42) for comparison. The agreement with experimental data is quite good at most points in the BZ, except for the frequencies of the optical modes around the zone center. This also leads to some discrepancies for the lower optic branches throughout the zone, while the acoustic branches are well reproduced.

These findings are in line with the previous theoretical study of the lattice dynamics of Bi by Díaz-Sánchez et al.(24) Comparing calculations without and with SOC, they found that SOC clearly improves the description of the phonon dispersion, but with the exception of the optical branches near the point. Overall the performance of DFT is poorer when compared with calculations for other heavy metals that are also largely influenced by the SOC.(8)

A prominent feature of both experimental and theoretical spectra is the pronounced dip in the dispersion of the two lower optical branches at , with a minimum frequency of about 9 meV found by inelastic neutron scattering and Raman experiments at low temperatures.(42); (43) Such dips also occur in other group V elements (As, Sb) as well as in IV-VI compounds like SnTe, and have been interpreted as caused by long-ranged interactions due to a resonant bonding effect in materials with pseudo-rocksalt structure.(44) Our calculation exaggerates this dip, a phenomenon well known from previous theoretical work.(24) We have checked that including the 5 semicore states in the valence space does not alter this feature. Furthermore, it was not changed by varying the Gaussian broadening, which indicates that Fermi surface related renormalization is not responsible for it. We observed, however, a high sensitivity on changes of the lattice structure, in particular on the internal structural parameter . Using the experimental value of instead of the optimized value raises the lower and upper optical frequencies at by 1.2 meV and 0.8 meV, respectively, thereby reducing the difference to the measured frequencies to less than 0.6 meV. Similar improvements are found for all optic branches.

We have performed calculations of the electron-phonon coupling constant in bulk Bi, given by

(1) |

Here, denotes the screened electron-phonon coupling matrix elements, and the sum runs over all phonon modes (with momentum , branch index , and frequency ) and over all electronic states (with momentum , band index , and energy ). denotes the electronic density of states per spin at the Fermi energy . Since the bulk Fermi surface consists of tiny Fermi pockets, care must be taken in the BZ sampling to assure a proper convergence. We have used meshes of up to 363636 points for the electronic states and approximated the delta functions by Gaussians with widths of 0.1–0.2 eV. For the sum over phonon modes, the 121212 -point mesh was applied. In all cases, we found , which is consistent with the absence of superconductivity in bulk Bi down to 50 mK.(6) In a recent linear-response calculation using GGA a very similar value of 0.09 was found.(5)

### iii.2 Bi(111) geometry

The (111) surface of Bi is that perpendicular to the trigonal axis. The surface can, in principle, be formed by breaking either the three bonds of each atom with its 2nd nearest neighbors (NN) or those with its 1st NNs. The latter case, however, would require much more energy than the former. Its low reactivity(45) also speaks in favor of breaking the weaker bonds with 2nd NNs. As such, Bi(111) can be considered as formed by a set of bilayers (see Fig. 4) where intra-bilayer bonds (between 1st NNs) are rather covalent (strong) and the vertical intra-bilayer spacing is relatively short (, , etc.) whereas inter-bilayer bonds (between 2nd NNs) are rather metallic (weak) and the vertical inter-bilayer spacing is relatively long (, , etc.). They correspond to the interlayer distances and of the bulk, respectively. Figure 4 shows that the stacking of planes in an ideally bulk-terminated Bi(111) surface (before relaxation) is such that every six planes the structure repeats itself. Each atom in any layer has three 1st NNs and three 2nd NNs in the two adjacent layers. There are three kinds of layers A(A’), B(B’), and C(C’). A and A’ layers have atoms lying on the line, B and B’ layers have atoms lying on , and C and C’ have atoms lying on . The stacking makes Bi(111) a quite open surface. The closest distance between atoms within a layer is 4.54 Å such that the surface leaves exposed even the third layer.

Theory | Experiment^{1} |
||||||
---|---|---|---|---|---|---|---|

6 BL | 10 BL | 12 BL | 6 BL^{2} |
7 BL^{3} |
T=140 K | T0 K | |

-1.62 | -1.56 | -1.56 | -0.83 | +0.6 | +0.51.1 | +1.22.3 | |

1.37 | 1.32 | 1.32 | +3.13 | +6.2 | +1.90.8 | +2.61.7 | |

-0.87 | -0.93 | -0.87 | 0.01.1 | ||||

0.13 | 0.18 | 0.13 | |||||

-0.25 | -0.31 | -0.31 |

In bulk Bi, the intra-bilayer distance is given by Å, while the vertical inter-bilayer distance is Å. Upon the creation of the surface the interlayer distances relax. Table 2 shows the calculated relaxation for three different thicknesses of the slab (6, 10 and 12 bilayers) together with previous calculations(21); (46) and low-energy electron diffraction (LEED) measurements.(46) Our results show a fast convergence of the relaxation with the slab thickness. Already for a 10-bilayer slab the interlayer distances are sufficiently converged, and the bulk spacing at the center of the slab is recovered. We found a shrinking of the first intra-bilayer distance () of 1.56% and an expansion of the first inter-bilayer distance () of 1.32%. The first result is at variance with the LEED experiment, where an expansion of was found, although the experimental error bars are large. Our value for , however, agrees with the LEED result. Our pseudopotential results contrasts a previous calculation using the full-potential linearized augmented plane wave method, where both distances were found to expand,(22) but a very large expansion of 6-7% was predicted for at variance with experiment. A recent pseudopotential calculation employing GGA found relaxations much closer to the present ones.(21)

### iii.3 Bi(111) electronic structure

As indicated in Fig. 5, the tiny electron and hole pockets of bulk Bi are projected onto regions around the and points of the (111) surface BZ, respectively. In between, there are larger regions where the bulk spectrum possesses an energy gap near the Fermi level. As shown in Fig. 6, the electronic bandstructure of Bi(111) exhibits two surface localized bands along and part of the line, which fall into this energy gap region. They are derived from two semi-relativistic bands, one from each side of the slab. Since the surface electronic bands show a rather poor convergence with slab thickness, we performed calculations of the electronic bandstructure for 12, 24, and 36 bilayers of Bi(111). As shown in the previous Subsection, these thicknesses are sufficient to converge the interlayer distances at the surface and to recover the bulk spacing in the center of the slabs. In contrast, the surface electronic bands are not converged, but exhibit a sizable dependence on the film thickness in particular when approaching the point. It was argued that for a semi-infinite surface, the two split bands should be degenerate at on symmetry grounds.(47); (22) However, in the case of finite slabs such a degeneracy is not found. There is a strong hybridization between the upper and lower surface that does not allow convergence of the energy of these bands around even for slabs of 36 bilayers or more.(22) In contrast, this inter-surface interaction is much weaker away from , such that both electron and hole pockets are well converged already for a 12-bilayer Bi(111) slab. A more detailed analysis of the electronic structure and its dependence on the slab thickness will be published elsewhere.(48) Due to such a slow convergence, we will restrict the analysis of the electron-phonon induced self-energy effects (Sec. III.5) to those parts of the lines, where the surface band energies are sufficiently converged for a 12-bilayer slab.

### iii.4 Bi(111) surface phonons

Point | This work | Experiment^{4} |
Theory | |
---|---|---|---|---|

5 BL ^{5} |
6 BL ^{6} |
|||

12.3 (L) | 11.8, 12.2 | |||

12.3 (SH) | ||||

13.3 (SV) | 14.2 | |||

13.6 (SV) | 14-15 | 13.6 | 15.3 | |

2.9 (RW, SV+L) | ||||

11.5 (SV+L) | ||||

12.7 (SV+L) | ||||

12.9 (SV+L) | ||||

3.4 (RW) | ||||

4.1, 4.6, 4.8, 5.9, 6.4 | ||||

11.1, 11.4, 11.8, 12.2 |

We now address the lattice dynamics of the Bi(111) surface. Recently, three theoretical studies were devoted to surface phonons in the framework of the density functional theory. (25); (5); (26) Due to the heavy numerical work involved in dealing with SOC, these calculations were performed only for thin slabs up to 6 bilayers, where the vibrational spectrum exhibited significant variations as a function of thickness. As we are primarily interested in the properties of a semi-infinity Bi(111) surface we checked the convergence of the surface localized vibrations by performing lattice dynamics calculations for 6- and 12-bilayer slabs, respectively. In both cases, the slabs were relaxed and SOC was taken into account. Surface phonon spectra were then obtained by combining the real-space force constants of the slab with those from the bulk calculation to simulate the dynamics of an asymmetric slab of 50 bilayers. An asymmetric slab contains one surface taken from the fully relaxed slab, while the other surface corresponds to the ideal bulk-truncated case. The phonon spectrum of such a slab thus contains surface modes related to both types of surfaces. From the components of each eigenvector, one can determine the vibrational weight at each atomic site and at each direction. Fig. 7 shows the phonon dispersions based on the 6 and 12-bilayer calculations, respectively. The optic and acoustic parts of the projected bulk bands (grey lines) slightly touch each other at , because the lower optic bulk mode at the zone center has such a small calculated frequency. Surface localized modes of the true (relaxed) surface are identified by their vibrational weight in the first two bilayers at the surface. Modes with a weight of more the 20% are highlighted in red, while those with a weight of more the 10% are indicated by black dots. A more detailed characterization of surface modes at the high-symmetry points , , and is given in Table 3.

One can note that after the slab filling procedure, both calculations give essentially the same surface localized spectrum. This indicates that changes in the real-space force constants due to the presence of the surface are essentially confined to the outer two or three bilayers and are reasonably well converged already for the 6-bilayer slab. The fast convergence of the dynamical properties contrasts the very slow convergence of the surface electronic band structure. This indicates that the uncertainty in the surface band structure near the point has little influence on the electronic response to atomic perturbations and hardly affects the frequencies of the surface phonons.

A prominent feature of the Bi(111) surface phonon dispersion is the appearance of modes above the optic bulk band almost everywhere along the and directions. The upper branch consists of predominately vertical vibrations of atoms in the first bilayer, and reflects a 12% stiffening of the force constants between the two layers forming the first bilayer and the contraction of the corresponding interlayer spacing (see Table 2). The second highest branch above the optic bulk spectrum also has a shear vertical polarization but is mainly localized in the 2nd and 3rd bilayer. A third surface branch with in-plane (shear-horizontal and longitudinal) polarization is present within the optic bulk spectrum and is localized over the first three bilayers.

In the acoustic spectrum, the Rayleigh wave (RW) falls slightly below the bulk band along , while it remains within the bulk continuum along . Along , above the RW, there appear several modes localized primarily in the first bilayer. The surface mode at about 6 meV stretching from halfway to involves vibrations of atoms in the 2nd bilayer with mainly shear vertical and to a lesser extend in-plane polarization. Along the directions, the longitudinal 1st-bilayer mode lying above the RW shows a weaker localization near , but when approaching it becomes more localized and acquires a stronger vertical component.

The surface phonons have been recently investigated via inelastic helium atom scattering (HAS) by Tamtögl et al.(21) In Fig. 7(b) we have added the reported frequencies taken at lower temperatures ( K). They clearly observed surface modes above the optic bulk band and interpreted their data as evidence for two branches with shear vertical polarization, in agreement with our calculation. Their frequencies are about 1 meV higher than our theoretical result, but also lies significantly above the experimental bulk maximum of 13.5 meV (see Fig. 3). In the acoustic part, good agreement between our calculations and the HAS results are found for the RW and the longitudinal 1st-bilayer mode as well as the 2nd-bilayer mode at about 6 meV.

In the HAS spectra, a variety of additional peeks were observed in the range of 4–5 meV, most prominent near . They were interpreted as caused by longitudinally polarized vibrations in the 3rd layer (2nd bilayer). This interpretation was based on the theoretical vibrational spectrum of a 6-bilayer slab and supported by calculations of the charge density variations induced by such vibrations, which were found to be large enough at the scattering point of the helium atom to explain the observed HAS peaks.

When we consider the 12-bilayer slab, we do find similar vibrational modes in this energy range with enhanced longitudinal polarization in the 3rd and 4th layer. However, this enhanced vibrational weight is a property of thin slabs only. When simulating the dynamics of thicker slabs by the slab-filling procedure, the weight becomes successively smaller with increasing thickness. It remains an open question, if the density oscillations related to these modes remain strong enough in the limit of thick slabs to explain the HAS peaks, which have been observed for Bi(111) samples of macroscopic thickness.

The presence of low-frequency optical modes causes the surface atoms to have an enhanced vibrational density of states at lower frequencies. Fig. 8 shows the layer-projected phonon density of states (LDOS) for the first, second, and twelfth layer of the 12-bilayer slab. The twelfth layer is deep enough so that it undergoes no relaxation (see Table 2) and thus its LDOS is representative of the phonon DOS of bulk Bi. Fig. 8 shows that the LDOS of the first two layers have an additional peak at the upper end of the optical spectrum as compared to the bulk due to the presence of the super-bulk surface modes. In contrast, the acoustic part of their LDOS is softened with respect to the bulk, which can be traced back to a softening of the nearest-neighbor coupling within the first layer. For each LDOS curve the Debye-like behavior applies in the beginning of the spectrum. From the initial quadratic behavior we deduce a Debye cutoff frequency of 8.9 meV for the bulk and of 7.3 meV for the first and second surface layer, respectively. These values are faily consistent with experimental Debye frequencies of about 10 meV for the bulk(13); (49) and 6.5 meV for the Bi(111) surface layer. The latter value was deduced from measurements of the mean-square atomic vibrational amplitudes of the surface atoms.(46) The Debye model is often used as an approximation for the Eliashberg function describing the coupling of an electronic state to the phonon system. We will show in the following Section, however, that the spectral shape of the coupling function can deviate significantly from the simple Debye form for the surface electronic states of Bi(111).

### iii.5 Electron-phonon interaction

The metallic character of the Bi(111) surface gives the opportunity of studying the -ph interaction for states which are well confined at the surface and subject to a strong spin-orbit interaction. To quantify the strength of the electron-phonon interaction, the dimensionless -ph coupling parameter is used:

(2) |

Here denotes an electron (hole) state momentum and band index, is the maximum phonon frequency and is the electronic state dependent Eliashberg spectral function corresponding to phonon emission (E) and absorption (A) processes:(50)

(3) |

The and signs in the delta function with electron energies correspond to phonon emission and absorption, respectively. The sum is carried out over final electronic states () and phonon modes ().

In the calculation of the -ph coupling a 12-bilayer Bi(111) film is used, and only those states in the surface electronic bands that do not differ in energy from those of the 24- and 36-bilayer calculation are considered, i.e. we restrict our analysis to momenta in the range ). To see how the surface electronic states couple to phonons we have calculated for different electron (hole) energies and momenta. The summation over phonons in Eq. 3 was carried out over 1296 wave vectors () in the surface BZ. The delta function with electronic energies was approximated by a first-order Hermite–Gaussian function with a smearing width in the range of 0.07–0.4 eV.

The electron-phonon parameter as a function of momentum is shown in Fig. 9 for the two surface bands. The electron-phonon coupling near the Fermi level is of intermediate strength (0.45, the lower state) unless the electronic state lies close to . It is markedly different from the very weak coupling constant of found at the Fermi level of bulk Bi (see Sec. III.1). The strength of -ph coupling for excited electrons on Bi(111) is also higher than the values of obtained for the surface electronic states near the Fermi level on the Bi(100) and Bi(110) surfaces.(4); (16)

Due to a small energy range covered by the surface states (0.15 eV around ) the coupling strength is not expected to vary considerably with electron momentum if the scattering within the surface bands dominates. That is the case at large electron momentum. As is evident from Fig. 9, for the calculated shows a weak dependence on electron momentum and the strength of the -ph interaction varies between 0.3 and 0.5 for both electronic bands, i.e. it is not large for either of the two. There is also an indication of a rather weak dependence on the electron binding energy, namely, the difference in the coupling strength between the two electronic bands at the same is small and does not exceed 0.15. In order to illustrate this as a function of binding energy is shown in Fig. 10.

However, when approaches the point the strength of -ph interaction increases significantly. It reaches up to in both surface bands. Furthermore the -ph coupling becomes very sensitive to the energy position of a hole (electron) state in the surface band. At the same electron momentum, in the lower surface band can be more than twice as large as in the higher band.

In the case of Bi(100) surface, the strong energy dependence of was explained by simple phase space arguments with the assumption of a weak energy dependence of the -ph matrix elements.(4) Here this argument is not valid. To illustrate that Fig. 10 shows the calculated density of electronic states (DOS). In the surface state energy range, the density of states is a smooth curve which does not reproduce the variation of . Only for hole states with large momentum in the upper surface band the strength of -ph coupling is determined to a certain extent by the available phase space.

The scattering processes leading to such a strong -ph interaction at small electron momenta can be derived from the corresponding spectral functions. Fig. 11 shows the Eliashberg functions calculated for surface states S1, S2, and S3 (see Figs. 9 and 10), which have similar binding energies but different electron momenta. Only the average of the emission and adsorption spectral function is shown because both parts nearly coincide.

A distinctive feature of all shown in Fig. 11 is that the lattice vibrations with small energies (and small wave vectors) get involved deeply into the scattering processes of electrons even though their contribution to the phonon density of states is negligible (Fig. 8). The type of transitions they are connected to depends on the momentum position of the electronic state in the surface bands. In the case of state S1, the peak in the spectral function at 1 meV is related to the nearest intraband transitions. With decreasing electron momentum (S1S2 S3) the small energy peak caused by the intraband scattering remains, but in addition the possibility of scattering to electronic states near the point increases. Such transitions via long wavelength phonons with small energies give a large contribution to the coupling constant according to its definition (Eq. 2). The latter contribution increases rapidly as the electron momentum approaches the point and is responsible for the sharp increase of at small electron momentum and at energies, where electronic states around the point are available for phonon-mediated transitions. These states are essentially bulk-like. Specifically, while the dispersion of the two bands of the surface states (Fig. 9) suggests that they rise above the Fermi level when approaching , in reality both bands quickly acquire a bulk-like character, and a true surface localized state at appears at much lower energies ( meV below ). Thus, states near involved in the small-momentum transitions are bulk-like. This low-frequency contribution clearly increases the coupling strength, but it is also associated with a larger uncertainty, because for the present 12 bilayer slab, the bulk-projected part is represented only by a few bands in this energy region. Calculations with much thicker slabs would be required for proper convergence of this contribution, which is currently not feasible.

Attempts to measure the coupling strength of Bi(111) surface states focused on the part of the hole pocket closest to , which corresponds approximately to our state S2. From measurements of the self-energy, Ast and Höchst extracted coupling constants of and when approximating the Eliashberg function by a Debye spectrum with either bulk ( meV) or surface Debye frequencies ( meV), respectively.(13) Kirkegaard et al. later argued that by taking into account the finite spectrometer energy resolution these values should be corrected down to couplings of the order 0.4.(16) Such a value was obtained by Gayone et al. from data on temperature dependent momentum distribution curves near the Fermi level crossing.(17) This analysis was also based on a Debye model for the Eliashberg function.

The unusual shape of the calculated spectral functions, in particular for the S2 and S3 states, renders the applicability of the Debye model rather questionable. As shown by the dashed lines in Fig. 11, the Debye approximation(51) misses the strong enhancement of at meV. On contrast, the Debye model might be more appropriate for an analysis in the momentum region , where this low-energy peak is not very pronounced. For example, the initial slope of the Eliashberg function for the S1 state is reasonably well described by the Debye model if the surface Debye frequency of 5 meV is used (Fig. 11).

## Iv Summary

We have performed a comprehensive investigation of the structural, electronic, lattice dynamical and electron-phonon coupling properties of the Bi(111) surface within density functional perturbation theory, taking into account the spin-orbit coupling consistently. Although some details of the electronic structure of Bi(111) depend on the slab thickness, we found that the lattice dynamics of Bi(111) is practically converged already for slabs of 6 bilayers. Changes in the dynamical couplings are confined essentially to interatomic bonds in the first two bilayers, where also the main structural relaxation occurs. The surface phonon spectrum exhibits super-bulk modes which lie above the bulk spectrum. In addition, the layer-projected vibrational density of states is enhanced at lower frequencies for the first bilayer, consistent with the observed enhanced average vibrational amplitudes of the surface atoms. A calculation of electronic state dependent coupling to phonons gave moderate coupling strengths of 0.45 for surface states with larger momenta. Surface states close to are predicted to have significantly higher couplings of the order of 1. This increase is connected to an enhanced scattering via phonons with long wavelengths and small energies into bulk-like electronic states near . The calculated coupling for states of the hole Fermi surface closest to is about twice as large as those deduced from experiment. This discrepancy may be partly due to insufficient thickness of the slab, but also may result from an inappropriate assumption in the experimental analysis. We found that the state dependent Eliashberg functions strongly deviate in shape from that of a Debye model, indicating that the use of the latter in an analysis of electronic self-energies from photoemission spectra is not justified and likely fails to give correct results for the coupling strength.

## V Acknowledgments

This work of MAO and TSR was partially funded by DOE Grant DE-FG02-07ER46354.

### Footnotes

- low-energy electron diffraction, Ref. (46)
- GGA, pseudopotential plane-wave method, Ref. (21)
- LDA, full-potential linearized augmented plane wave method, Ref. (46)
- Helium atom scattering data from Tamtögl et al., Ref. (21), taken from their low-temperature data of Fig. 4(a)
- LDA pseudopotential calculation by Chis et al., Ref. (25). Values shown for are estimates based on their Fig. 3
- GGA pseudopotential calculation without spin-orbit interaction, Ref. (21)

### References

- Corresponding author
- D. Shoenberg, Proc. Roy. Soc. A. Math. and Phys. Sci. 170, 341 (1939).
- B. Weitzel and H. Micklitz, Phys. Rev. Lett. 66, 385 (1991).
- J. E. Gayone, S. V. Hoffmann, Z. Li, and Ph. Hofmann, Phys. Rev. Lett. 91, 127601 (2003).
- G. Q. Huang and J. Yang, J. Phys.: Condens. Matter 25, 175004 (2013).
- M. Tian, N. Kumar, M. H. W. Chan, and T. E. Mallouk, Phys. Rev. B 78, 045417 (2008).
- T. T. Chen, J. D. Leslie, and H. J. T. Smith, Physica 55, 439 (1971).
- R. Heid, K. P. Bohnen, I. Y. Sklyadneva, and E. V. Chulkov, Phys. Rev. B 81, 174527 (2010).
- Ph. Hofmann, Prog. Surf. Sci. 81, 191 (2006).
- G. Jezequel, Y. Petroff, R. Pinchaux, and F. Yndurain, Phys. Rev. B 33, 4352 (1986).
- M. Hengsberger, P. Segovia, M. Garnier, D. Purdie, and Y. Baer, Eur. Phys. J. B 17, 603 (2000).
- Ch. R. Ast and H. Höchst, Phys. Rev. Lett. 87, 177602 (2001).
- Ch. R. Ast and H. Höchst, Phys. Rev. B 66, 125103 (2002).
- Ch. R. Ast and H. Höchst, Phys. Rev. B 70, 245122 (2004).
- Yu. M. Koroteev, G. Bihlmayer, J. E. Gayone, E. V. Chulkov, S. Blügel, P. M. Echenique, and Ph. Hofmann, Phys. Rev. Lett. 93, 046403 (2004).
- C. Kirkegaard, T. K. Kim, and Ph. Hofmann, New J. Phys. 7, 99 (2005).
- J. E. Gayone, C. Kirkegaard, J. W. Wells, S. V. Hoffmann, Z. Li, and Ph. Hofmann, Appl. Phys. A 80, 943 (2005).
- Ph. Hofmann, I.Yu. Sklyadneva, E.D.L. Rienks, and E.V. Chulkov, New J. Phys. 11, 125005 (2009).
- I. Y. Sklyadneva, G. Benedek, E. V. Chulkov, P. M. Echenique, R. Heid, K. P. Bohnen, and J. P. Toennies, Phys. Rev. Lett. 107, 095502 (2011).
- A. Tamtögl, M. Mayrhofer-Reinhartshuber, N. Balak, W. E. Ernst, and K. H. Rieder, J. Phys.: Condens. Matter 22, 304019 (2010).
- A. Tamtögl, P. Kraus, M. Mayrhofer-Reinhartshuber, D. Campi, M. Bernasconi, G. Benedek, and W. E. Ernst, Phys. Rev. B 87, 035410 (2013).
- Y. M. Koroteev, G. Bihlmayer, E. V. Chulkov, and S. Blügel, Phys. Rev. B 77, 045428 (2008).
- É. D. Murray, S. Fahy, D. Prendergast, T. Ogitsu, D. M. Fritz, and D. A. Reis, Phys. Rev. B 75, 184301 (2007).
- L. E. Díaz-Sánchez, A. H. Romero, and X. Gonze, Phys. Rev. B 76, 104302 (2007).
- V. Chis, G. Benedek, P. M. Echenique, and E. V. Chulkov, Phys. Rev. B 87, 075412 (2013).
- J. Yang, G. Q. Huang, and X. F. Zhu, Phys. Stat. Solidi B 250, 1937 (2013).
- G. Benedek, M. Bernasconi, K.-P. Bohnen, D. Campi, E. V. Chulkov, P. M. Echenique, R. Heid, I. Yu. Sklyadneva, and J. P. Toennies, Phys. Chem. Chem. Phys. 16, 7159 (2014).
- L. Hedin and B. I. Lundqvist, J. Phys. C 4, 2064 (1971).
- J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
- D. Vanderbilt, Phys. Rev. B 32, 8412 (1985).
- S. G. Louie, K. M. Ho, and M. L. Cohen, Phys. Rev. B 19, 1774 (1979).
- B. Meyer, C. Elsässer, M. Fähnle, FORTRAN90 Program for Mixed-Basis Pseudopotential Calculations for Crystals, Max-Planck-Institut für Metallforschung, Stuttgart (unpublished).
- L. Kleinman, Phys. Rev. B 21, 2630 (1980).
- W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in Fortran, The Art of Scientific Computing, vol. 2nd Ed. (Cambridge Cambridge University Press,1992).
- N. E. Zein, Fiz. Tverd. Tela (Leningrad) 26, 3028 (1984) [Sov. Phys. Solid State 26, 1825 (1984)].
- S. Baroni, P. Giannozzi, and A. Testa, Phys. Rev. Lett. 58, 1861 (1987).
- R. Heid and K. P. Bohnen, Phys. Rev. B 60, R3709 (1999).
- P. Giannozzi, S. de Gironcoli, P. Pavone, and S. Baroni, Phys. Rev. B 43, 7231 (1991).
- R. Heid and K. P. Bohnen, Phys. Rep. 387, 151 (2003).
- M. S. Dresselhaus, The physics of semimetals and narrow-gap semiconductors: Proceedings (Supplement No. 1 to the Journal of physics and chemistry of solids, v. 32, Pergamon Press, 1971).
- D. Schiferl and C. S. Barrett, J. Appl. Cryst. 2, 30 (1969).
- R. E. MacFarlane, The physics of semimetals and narrow-gap semiconductors: Proceedings (Supplement No. 1 to the Journal of physics and chemistry of solids, v. 32, Pergamon Press, 1971).
- J. Höhne, U. Wenning, H. Schultz, and S. Hüfner, Z. Physik B 27, 297 (1977).
- S. Lee, K. Esfarjani, T. Luo, J. Zhou, Z. Tian, and G. Chen, Nat. Commun. 5, 3525 (2014).
- F. Jona, Surf. Sci. 8, 57 (1967).
- H. Mönig, J. Sun, Y. M. Koroteev, G. Bihlmayer, J. Wells, E. V. Chulkov, K. Pohl, and P. Hofmann, Phys. Rev. B 72, 085410 (2005).
- T. Hirahara, T. Nagao, I. Matsuda, G. Bihlmayer, E.V. Chulkov, Yu.M. Koroteev, P. M. Echenique, M. Saito, and S. Hasegawa, Phys. Rev. Lett. 97, 146803 (2006).
- M. Alcántara Ortigoza et al., to be published.
- K.G. Ramanathan and T.M. Srinivasan, Phys. Rev. 99, 442 (1955).
- G. Grimvall, The Electron–Phonon Interaction in Metals (North-Holland, New York, 1981).
- we used the 2D model of: B. Hellsing, A. Eiguren, and E V Chulkov, J. Phys.: Condens. Matter 14, 5959 (2002).