THz Surface Modes and Electron-Phonon Coupling in Bi$_2$Se$_3$(111)

We present a combined experimental and theoretical study of the surface vibrational modes of the topological insulator (TI) Bi$_2$Se$_3$ with particular emphasis on the low-energy region below 10 meV that has been difficult to resolve experimentally. By applying inelastic helium atom scattering (HAS), the entire phonon dispersion was determined and compared with density functional perturbation theory (DFPT) calculations. The intensity of the phonon modes is dominated by a strong Rayleigh mode, in contrast to previous experimental works. Moreover, also at variance with recent reports, no Kohn-anomaly is observed. These observations are in excellent agreement with DFPT calculations. Besides these results, the experimental data reveal$-$via bound-state resonance enhancement$-$two additional dispersion curves in the gap below the Rayleigh mode. They are possibly associated with an excitation of a surface electron density superstructure that we observe in HAS diffraction patterns. The electron-phonon coupling paramenter $\lambda$ = 0.23 derived from our temperature dependent Debye-Waller measurements compares well with values determined by angular resolved photoemission or Landau level spectroscopy. Our work opens up a new perspective for THz measurements on 2D materials as well as the investigation of subtle details (band bending, the presence of quantum well states) with respect to the electron-phonon coupling.


I. INTRODUCTION
Bi 2 Se 3 ( Figure 1) is categorised as a three-dimensional topological insulator, a new state of quantum matter with a bulk gap and spin-orbit split surface states forming a Dirac cone across the gap 1,2 . The interaction of electrons with surface phonons in Bi 2 Se 3 has been studied intensively due to its important role in transport properties and possible constraints for potential applications in a variety of nanotechnologies [3][4][5][6][7][8][9] . Bismuth selenide as well as telluride alloys are classical thermoelectric materials 10,11 with a large Seebeck coefficient and, as such, they have been used in thermoelectric refrigeration for a long time 12 . However, to fully understand the thermoelectric properties of Bi 2 Se 3 thin films and nanoscale devices 13,14 , information on the surface phonon dispersion curves and the electron-phonon (e-ph) interaction is crucial 13,15,16 . So far experimental information about the surface phonon dispersion curves of Bi 2 Se 3 (111) was limited to previous helium atom scattering (HAS) studies by Zhu et al. 7,17 , in the low energy part of the phonon spectrum. These studies suggested the presence of a deep Kohn anomaly in the 7.5 meV optical phonon branch (S2) coupled to the electronic (spin-forbidden) transition across the Dirac cone 17 . However, existing first-principle calculations of Bi 2 Se 3 (111) phonon dispersion curves, do not show any evidence of Kohn anomalies in the S2 branch 18 . A convenient parameter to characterise the e-ph coupling strength is the mass-enhancement λ 21 and in recent years it was demonstrated that HAS from conducting surfaces can directly measure the mode-selected e-ph cou-pling constants λ Qj 22,23 , besides the surface phonon dispersion curves 24 . Moreover the temperature-dependence of the HAS Debye-Waller factor was shown to provide the global e-ph coupling constant λ at the surface of thin metal films 24,25 and topological insulators 26 . Yet the large λ as estimated based on the Kohn anomaly 7 is at odds with theoretical findings that indicate that the major contribution to λ comes from the higher optical phonon branches 18 , especially when the Fermi level crosses the surface quantum-well states above the conduction band minimum (see Figure 1(c)). The situation is actually met in recent high-resolution 3 He-spin scattering studies on Bi 2 Te 3 (111), where the weak signature of a Kohn anomaly is detected in the surface longitudinal acoustic resonance 27 , also not found in adiabatic ab-initio calculations of the phonon branches. In order to elucidate these conflicting results, we have undertaken a HAS study of the surface phonon dispersion curves and the e-ph interaction of Bi 2 Se 3 (111). Supersonic neutral He atom beams with incident energies in the range ≤ 20 meV have been used to probe low-energy surface excitations with the best available resolution, while being an inert completely nondestructive probe 24,28 . The technique allows to measure most of the surface phonon branches in the acoustic and optical regions. Low-energy He atoms impinging on a conducting surface are exclusively scattered by the surface charge density 23,29 and inelastic scattering from surface phonons only occurs via the phonon-induced charge density oscillations, i.e., via the e-ph interaction. It is in this way that inelastic HAS provides a first-hand information on the e-ph interaction, with the neutral He atoms acting as a sort of local me-  19,20 ). Despite the Dirac cone, quantum-well states as illustrated by the orange lines exist on the surface as well. The quantum-well states are split due to Rashba coupling, thus making the previously degenerate spin states (illustrated by the red and blue arrows) separate into an inner and outer branch, each having a spin texture with opposite chirality. chanical probe on the electron density. Energy and momentum, inelastically exchanged by He atoms with the surface can, however, be retained by the electron system in the form of low-energy collective excitations. In principle, the HAS signal from this kind of excitations is expected to be quite small. Nevertheless, an increased e-ph interaction due to surface quantumwell states 30 in combination with an enhancement from HAS bound-state resonances 31 , suggests to assign two branches of low-energy modes in the gap well below the Rayleigh waves (RW) to some sort of collective electronic excitations. Actually anomalous acoustic plasmons have been recently reported in Bi 2 Se 3 (111) by Jia et al. 32 , from high-resolution electron energy-loss spectroscopy, although these modes turn out to be superimposed in the first Brillouin zone onto the RW branch. Plasmons in a two-dimensional electron gas (2DEG) with a √ Q dispersion (2D plasmons) have been predicted long ago by Frank Stern 33,34 . Later it was shown that the coupling of 2DEG plasmons arising from two different quantum-well minibands, as found in semiconductor sur-face accumulation layers, yield a surface plasmon pair: a 2D plasmon and an acoustic surface plasmon (ASP) with a linear dispersion above the upper edge of the singleparticle excitation spectrum 35,36 . Similarly the coupling of a 2DEG at a metal surface coupled to the underlying 3D electron gas yields an ASP in addition to the ordinary surface plasmon 30,[37][38][39] . As discussed below, the assignment of the two additional low-energy branches as collective polaron excitation recently suggested by Shvonski et al. 35 , although plausible in semimetals with a large dielectric constant, definitely requires further ad-hoc studies, possibly with even higher resolution as available, e.g., with 3 He spin-echo spectroscopy 27 .   Figure 2(b) the temperature dependence of the specular peak intensity is plotted, which can be used to extract the electron-phonon coupling constant λ (section ). The diffraction scans in Figure 2(a) have been measured at three different incident energies E i , at a sample temperature of 113 K. The intensity scale has been scaled to show additional features with smaller intensity. Besides some features assigned to bound-state resonances and kinematic focusing, which are easily recognised due to the strong dependence of their position on the incident energy 29 , there are features that occur at fixed values of ∆K, independently of E i with a distance of about 0.2Å −1 from the specular and first order diffraction peaks, as indicated by the vertical shaded regions (further diffraction scans, including also the ΓK azimuth can be found in the supplementary information). We recently observed with HAS a multivalley charge density wave (CDW) in Sb(111) originating from the M-point electron pockets giving rise to additional peaks in the diffraction pattern 41 . In Bi 2 Se 3 (111), however, no carrier pockets exist besides the Dirac cone and the quantum-well minibands occurring around the zone center (see Figure 1(c)). The latter provide nesting wavevectors of about 0.2Å −1 between states of equal spin, which correspond fairly well to the parallel momentum transfers of the satellites observed aside the (0, ±1) peaks in HAS diffraction spectra (and likely also aside the specular peak, despite the coincidence with bound-state resonances) (Figure 2(a)). It should be noted, however, that the observation of satellite peaks whose position is independent of the HAS incident energy is by itself indicative of a long-period superstructure of the electron density, possibly incommensurate or weakly commensurate with the surface atomic lattice.
Charge density oscillations as low as 10 −6 atomic units, presently accessible to HAS, can in principle sustain very low-energy collective phase and amplitude excitations in the meV spectral range 41 , and possibly suggest an assignment of the present low-energy branches.
A. Time-of-flight measurements and phonon dispersion curves The phonon energies were determined by performing time-of-flight (TOF) measurements over a wide range of incident angles between the first-order diffraction peaks and at various beam energies. The phonon dispersion was then obtained by calculating the parallel momentum transfer |∆K| for each extracted phonon energy from the conservation laws of energy and parallel momentum providing the so-called scan curve for planar scattering (see (1) and Refs. 23,42 ). In Figure 3(a) an example of a TOF spectrum after conversion to the energy transfer scale is shown. The measurement was taken in the high symmetry direction ΓM with an incident beam energy E i = 17.97 meV and at an incident angle of ϑ i = 34.25 • . The TOF spectrum consists of several peaks which are located on the creation (∆E < 0) as well as the annihilation (∆E > 0) side. The peak at zero energy transfer corresponds to elastically scattered helium atoms 27,28 . The scan curve, shown in the two centre panels of Figure 3 for phonon annihilation (blue) and creation (red) events, associates each phonon event with a specific momentum transfer ∆K. The scan curve has been backfolded into the irreducible part of the Brillouin zone and is plotted on top of the calculated dispersion. The different symbols on the scan curves, marking the main inelastic features, have been associated to phonons of different character and polarisation. The large peaks in the TOF spectra marked with the red circles in Figure 3, correspond to the Rayleigh wave (RW) as seen in the DFPT calculations. Note that in the present TOF spectra the RW exhibits typically the largest intensity of all inelastic events (cfr. the intensities in Figure 3(a)). There is a fair correspondence between the present HAS data and those previously reported by Zhu et al. 17 . Curiously Zhu et al. stated that the RW is not observed, whereas it appears in their plot, though with only a few data points, in reasonable agreement with the present one in the ΓK direction; it also occurs in the ΓM direction, once it is recognised that there is an avoided crossing, so that the RW at M is not the lowest mode. There is however an important difference with respect to Zhu et al. 17 : Present HAS data do not show any evidence of a Kohn anomaly in the ≈ 8 meV branch for wavevectors around 0.2Å −1 and associated with the nesting at the Fermi level across the Dirac cone (or more likely across the parabolic dispersion of surface quantumwell states 19,20 ). Figure 4 shows the entire experimental surface phonon dispersion (symbols) superimposed onto the DFPT calculations (grey lines). The different symbols have been associated to phonons of different character and polarisation based on the proximity to particular modes of the DFPT calculations. In total, we are able to distinguish at least 8 different branches. The polarisation analysis of the calculated surface phonon modes can be found in Figure 5 where the intensity of each polarisation projected onto the corresponding layer is given by the colour code. The left column shows the longitudinal polarisations for the first (L1), second (L2), and third (L3) layer. The right column shows the shear vertical polarisation for the first three layers (SV1-SV3), while the shear horizontal polarisation can be found in the supporting information. The theoretical dispersion curves are seen to agree quite well with the HAS data and also with the experimental Raman-active modes at Γ (green triangles in Figure 5 according to 43 ). A closer comparison of the experimental data points in Figure 4 with Figure 5 shows that mainly phonon events with the largest amplitude in the two topmost layers of the sample are observed in the experiment. In particular in the low energy region (< 10 meV), most contributions The blue (phonon annihilation, ∆E > 0) and red (phonon creation, ∆E < 0) lines in the centre panels show the scan curves superimposed onto the calculated dispersion (The colour code giving the intensity of the longitudinal L1 and L3 modes projected onto the first and third surface layer, respectively, is that of Figure 5). The symbols denote peaks in the TOF spectrum which have been assigned to phonon events (RW = Rayleigh wave). The two distinct low-energy peaks (green diamonds) observed on the creation (panel (a)) and annihilation (panel (b)) sides in the phonon gap below the RW are tentatively assigned to low-energy collective surface electronic excitations, associated with the long-period surface superstructure and revealed through HAS bound-state resonance enhancement.
come from phonons with the largest amplitude in the second layer (L2, SV2), which is a Bi layer and therefore about 2.5 times heavier than the first Se layer. The most prominent mode in the TOF spectra, the RW, corresponds predominantly to L and SV polarisations, due to its elliptical polarisation, with a particularly strong  18 suggested that the Kohn anomaly observed by Zhu et al. 7,17 may be connected to a strong e-ph interaction in the doped bulk material rather than to a surface state. They actually showed that the largest contribution to the e-ph coupling comes from an in-plane polar-type branch in the optical region between 10 and 18 meV 18 . Indeed the anomalously strong dispersion of the optical branches in that region ( Figure 5), also found in Bi and Sb tellurides 27,44 may be regarded as a manifestation of e-ph interaction. In the acoustic region and the long wavelength limit (close to Γ) the dispersion relation of the RW is linear. Its slope provides the RW group velocity in the two symmetry directions ΓM and ΓK: It appears that the RW has a velocity in the ΓM direction smaller than both transverse bulk values and is therefore a localised surface wave, whereas in the ΓK direction it has a velocity larger than that of the SH transverse sound, and is therefore a pseudo-surface wave (PSW) 46,47 . Actually in the absence of mirror symmetry for the sagittal plane in this direction, the RW is a resonance. The fact may have suggested (see Zhu et al. 7,17 ) that in Bi 2 Se 3 (111) the RW is suppressed but the present comparison with the DFPT calculation confirms that the RW is actually observed in both directions, though as a resonance along ΓK. Values for the bulk longitudinal (v L = 2.9 km/s) and transverse (v T = 1.7 km/s) group velocities of Bi 2 Se 3 have also been reported in the framework of the isotropic elastic continuum theory 8,48 . In this approximation the corresponding RW velocity, obtained by solving the cubic Rayleigh equation, 49 would be v RW = 1.56 km/s in any direction.

B. Low-energy branches
The measured HAS-TOF spectra displayed in Figure 3 show also distinct peaks yielding two branches of elementary excitations with an energy below the RW branch (green diamonds in Figure 4). On the basis of present DFPT surface phonon dispersion curves, they cannot be attributed to any possible phonon branch of the ideal surface. HAS from conducting surfaces exclusively occurs through the interaction (mostly Pauli repulsion) with the surface electron density, and therefore also electronic excitations in the THz range can be observed by HAS, with a 0.5 meV resolution and sensitivity to charge density oscillations in the 10 −6 atomic units range. Actually the observed low-energy branches are reminiscent of those recently observed with HAS in Sb(111), which have been attributed to elementary excitations (phasons/amplitons) of a multi-valley CDW 41 . The concomitant presence of a commensurate component asso-ciated with the M-point electron pockets at the Fermi level, and an incommensurate one due to the hole pockets along the ΓM direction, allows for collective excitations with a comparatively small gap at Q = 0. On the other hand no low-energy phason/ampliton modes have been detected with HAS for the perfectly commensurate multivalley CDW reported in the quasi-1D surface Bi(114) 50 , discommensuration being a requisite for depinning and a vanishing/small gap at Γ. Bi 2 Se 3 (111) has no pocket states at the Fermi level, besides the rings around Γ of the surface topological Dirac and the quantum-well states 19,20 . The satellites near the HAS diffraction peaks (see Figure 2(a)) suggest some longperiod charge-density structures and possibly low-energy collective excitations. In order to detect the associated, seemingly small inelastic features in the TOF spectra, we rely on the bound-state resonance enhancement method (Ref. 31 and Chap. 10 of Ref. 24 ), applicable to highly corrugated surfaces and successfully used to detect with HAS high-energy optical surface modes in ionic crystals 24,51 . The complete set of He-surface bound states has been measured previously 40 . Bound-state inelastic resonances occur in the HAS-TOF spectrum, with a possibly large enhancement of the inelastic intensities, at the locus of intersections of the scan curve (1) 23 with the inelastic bound-state resonance condition (2) (see Chap. 30 of Ref. 24 ). For elementary excitations with an energy ∆E = ε and wavevector ∆K = Q the equations become: (2) At the mentioned intersection of (1) and (2), an elementary excitation with (ε, Q) assists the selective adsorption of the atom of mass m, incident energy E i , wavevector k i and angle ϑ i into a bound state of energy − |E 0n |, via the exchange of a surface reciprocal lattice vector G = (G , G ⊥ ). On returning the G-vector to the surface lattice, the atom is selectively desorbed from the bound state into the final angle ϑ f . In Equation 2 the vector G has been conveniently expressed via its components parallel and orthogonal to the scattering plane, respectively. In Bi 2 Se 3 (111) the measured He-surface boundstate energies 40 are E 0n = −5.6, −3.8, −2.3, −1.2, −0.5 meV for n = 0, 1, 2, 3, 4, respectively. For a fixed scattering geometry ϑ i + ϑ f = ϑ SD (here ϑ SD = 91.5 • ), Eqs.
(1,2) provide, via the elimination of ϑ i = ϑ SD − ϑ f , the locus of intersections ε = ε Ei,n;G (Q) for any incident energy E i , bound state n and reciprocal surface vector G. The three panels of Figure 6 show some plots of ε = ε Ei,n;G (Q) in the ΓM direction for: (a) a given bound state (n = 2), a G vector (1,0) and different values of the incident energy E i ; (b) a given incident energy E i = 10 meV and different bound state energies; (c) a given incident energy E i = 10 meV and bound state n = 2 and some different G-vectors, whose functions ε = ε Ei,n;G (Q) cross the phonon gap below the RW in the first BZ. In practice the phonon gap can be fully scanned by the resonance curves ε = ε Ei,n;G (Q) by varying the incident energy, so as to detect, via resonance enhancement, weak elementary excitations. Since the low-energy data points appear to allign along two dispersion curves, independently of the incident energy, as well as of n and G, rather than being spread over the entire gap, they cannot be attributed to a resonance-enhanced many-phonon background. Furthermore, frustrated translational modes of adsorbates like CO would show no dispersion and would appear at higher vibrational energies 52 . More likely these points indicate two branches of low-energy excitations associated with the surface charge-density superstructure observed in the diffraction spectra, as anticipated above. In this respect it is worth mentioning a recent work by Shvonski et al. 35 where it is argued that a strong e-ph interaction affecting the surface 2DEG of a 3D topological crystal allows for collective polaron excitations (plasmonpolarons). Their dispersion is predicted to be that of an acoustic plasmon running below the single-particle excitation spectrum as an effect of the polaron-polaron attractive interaction. The theoretical analysis by Shvonsky et al. 35 is actually interpreting the recent observation with high-resolution electron energy loss spectroscopy (HREELS) by Jia et al. 32 of an anomalous acoustic plasmon (AAP) mode from the topologically protected states of Bi 2 Se 3 (111), with energy between 0 and 6.5 meV (and its continuation in the second zone up to ≈ 10 meV). The present HAS data do not permit to identify this AAP due to its superposition with the RW in the first BZ and in part with other phonon branches in its continuation.

C. Electron-phonon coupling
As shown in recent papers 25,26 , the temperature dependence of the Debye-Waller (DW) exponent plotted in Figure 2(b) permits to extract for a conducting surface the mass-enhancement parameter λ expressing the electron-phonon coupling strength. It is related to the DW exponent by the equations: where φ = 4.9 eV is the work function 53 , A c = 14.92Å 2 the unit cell area, I(T S ) the He-beam specular intensity, T S the surface temperature, k iz = 3.18Å −1 the normal component of the incident wavevector, and n s the number of conducting layers which contribute to the phononinduced modulation of the surface charge density 54 . The latter is estimated to be n s = 2λ T F /c 0 , where λ T F is the Thomas-Fermi screening length characterising the surface band-bending region (here ≈ 6 nm) 19 , c 0 = 9.6Å the quintuple layers (QL) thickness, and the factor 2 indicates the 2DEG multiplicity as observed with ARPES in the current Bi 2 Se 3 sample 19 .
With these values and the experimental DW derivative with respect to T S from Figure 2(b), we obtain λ = 0.51. It should be noted that, unlike in the case of low-index metal surfaces, characterised by a soft-wall repulsive potential and negligible corrugation, here the large electronic corrugation 40 implies a hard-wall potential. In this case one needs to correct k 2 iz so as to account for the acceleration impressed by the attractive part of the potential on the He atom when approaching the surface turning point (Beeby correction 28 ). Therefor k 2 iz is replaced by k ′ 2 iz = k 2 iz + 2mD/ 2 , where m is the He mass and D = 6.54 meV the He-surface potential well depth 40 . With the Beeby correction it is found λ = 0.23. The value compares quite well with values in the literature derived from other experiments, e.g., λ = 0.25 3 , and λ = 0.17 5 from ARPES measurements and λ = 0.26 55 from Landau level spectroscopy. A theoretical study by Giraud et al. 8 with phonons calculated in the isotropic continuum limit gives λ = 0.42, whereas for other ARPES measurements, where only Dirac states appear to be involved, values of λ as low as 0.076 to 0.088 have been found 4 . From the comparison it appears that the presence of a 2DEG due to quantum-well minibands (at least two in the present analysis) plays an important role in raising the e-ph coupling strength, which is quite small when exclusively due to the Dirac topological states, to values in the range of 0.2 − 0.4. The same conclusion follows from the theoretical analysis by Heid et al. 18 , who showed that raising the Fermi level from the Dirac point to above the conduction band minimum gives a corresponding increase of λ from values well below 0.1 to values in the range above 0.2, with a substantial contribution from interband coupling and in very good agreement with the present analysis. The role of n-type doping contributing to the formation of the surface quantum-well 2DEG is quite clear in the analysis of the e-ph coupling strength in Cu-doped Bi 2 Se 3 , where an analysis based on the McMillan formula 56 , indicates a value for λ as large as 0.62 4 .

III. CONCLUSIONS
In summary, we have determined the surface phonon dispersion curves of Bi 2 Se 3 along both high symmetry directions, where the largest inelastic scattering intensity is provided by the Rayleigh wave. Thus our measurements show in contrast to previous studies that the Rayleigh mode exists and is a localised surface mode in one of the high-symmetry directions (ΓM), while in the other highsymmetry direction it is actually a pseudo-surface wave (ΓK). Comparison with density functional perturbation theory calculations shows excellent agreement with the experimental data. In addition to the phonon-related losses, we observe two additional dispersion curves in the gap well below the Rayleigh mode. These two low-energy branches may correspond to collective low-energy excitations of surface electrons. The appearance of these collective electronic excitations in an unprecedentedly low energy region is probably associated with a small surface charge density and an appreciable electron-phonon coupling (λ = 0.23). However, much more detailed experiments and theoretical analysis will be needed in order to fully understand these excitations; e.g., what is the influence of the carrier concentration upon doping and what is the role of both the Dirac and the quantum-well states, with the latter providing a much larger electron-phonon interaction than the former. The analysis advocates for a more systematic study by means of elastic and inelastic HAS spectroscopy of the surface structure, the low-energy collective excitations, and the electron-phonon interaction of interesting 2D materials, where the superior space and energy resolution of HAS is hardly attainable with other current surface probes.

A. Experimental Details
The reported measurements were performed on a HAS apparatus which generates a nearly monochromatic beam (∆E/E ≈ 2%) of 4 He that is scattered off the sample surface in a fixed 91.5 • source-sample-detector geometry. The beam is produced in a supersonic expansion of He through a 10 µm nozzle followed by sampling the core of the supersonic expansion via a 310 µm skimmer. For a detailed description of the apparatus please refer to 57 . Energy dispersive measurements for inelastic scattering can be performed using TOF measurements with a pseudo-random chopper disc. After deconvolution with the pseudo random chopper sequence, the TOF signal is further transformed to an energy transfer scale which allows to determine inelastic (phonon) scattering events 57 . The scattering spectra were mainly taken with the crystal at room temperature, while a few spectra were taken with the sample cooled down to 115 K. The incident He beam energy was varied between 10 and 18 meV. The crystal structure of Bi 2 Se 3 is rhombohedral, formed of QL which are bound to each other through weak van der Waals forces 58 . The hexagonal unit cell of the Bi 2 Se 3 crystal, shown in Figure 1, consists of 3 QLs. Each QL is terminated by Se atoms, giving rise to the (111) cleavage plane that exhibits a hexagonal structure (a = 4.14Å at room temperature 59 , see Figure 1(b)). The Bi 2 Se 3 crystal was attached onto a sample holder using thermally conductive epoxy. The sample holder was then inserted into the transfer chamber using a load-lock system and cleaved in-situ 60 . The sample can be heated using a button heater on the backside of the crystal or cooled down to 115 K via a thermal connection to a liquid nitrogen reservoir. The sample temperature was measured using a chromel-alumel thermocouple.

B. Computational Details
The surface dynamical properties of Bi 2 Se 3 were studied using DFPT calculations 61 withing the Quantum-ESPRESSO package 62 .
Norm-conserving pseudopotentials and the Perdew-Burke-Ernzerhof (PBE) approximation 63 for the exchange and correlation functional were used as implemented in the Quantum-ESPRESSO package. The surface phonon dispersion was calculated using a slab consisting of 3 QLs separated from its periodic replica by 20Å of vacuum, without the inclusion of spin-orbit corrections (SOC) (see also the SI for calculations with SOC). For an accurate calculation of the surface lattice vibrations, in principle both SOC and van der Waals (vdW) corrections are necessary, both due to the presence of heavy elements in the compound and the latter to fully account for the weak bonds between the individual quintuple layers. However, as thoroughly discussed for Bi 2 Te 3 (111) 27 , it appears that for layered crystals with heavy elements SOC alone gives a general softening of the phonon spectrum, compensated by the inclusion of vdW correction, so that satisfactory results are obtained at a minor computational cost without both SOC and vdW corrections and with a better agreement with experiment 23,64 . Also for Bi 2 Se 3 better agreement with the experiment is achieved with no SOC and no vdW corrections (see section II A), likely due to a compensation of errors between the underbinding often characterising PBE functionals and SOC contributions and the overbinding due to vdW forces 27 . More precisely, the effect of SOC was found however to be weak for the low energy surface vibrational modes of typical TIs such as Bi 2 Te 3 and Sb 2 Te 3 27,44 while on the other hand it was shown that vdW corrections become important for an exact description of the low energy optical modes of Bi 2 Te 3 27 .

ASSOCIATED CONTENT
Supporting Information accompanies this paper including additional DFPT calculations, as well as additional spectra and further details about the helium atom scattering experiments.

S1. DFPT CALCULATIONS AND THE EFFECT OF SOC
The topologically protected Dirac cone forms a small Fermi circle around the Γ-point which cannot be properly described with a coarse mesh and a large smearing resulting in the impossibility to capture subtle effects, such as the proposed Kohn anomaly 17 with a standard calculation. In order to verify the effect of such states on surface phonons we repeated the phonon calculations at the q-point corresponding to the nesting vector (2k F ) by including SOC. We compared the result with that obtained by omitting the SOC. To perform these calculations we had to improve the sampling of the Brillouin zone close to the Fermi surface which FIG. S1. Surface phonon dispersion of Bi2Se3 omitting (black continuous line) and including spin-orbit coupling (red dots) at a q-vector corresponding to the nesting vector 2kF . The inclusion of SOC gives rise to a softening of the phonon modes but no evidence for a Kohn anomaly is found.
is particularly important to resolve a possibly existing anomaly. Given the peculiar shape of the Fermi surface in the Bi 2 Se 3 (111) slabs consisting of a ring around the Γ-point, we used a graded k-point mesh (equivalent to a 50 × 50 × 1 uniform mesh) near the Γ-point point and a coarser one (equivalent to a 8 × 8 × 1 mesh) near the zone boundary. The results are reported in Figure S1. A one to one comparison between phonon modes calculated with and without SOC shows that there is no evidence of a Kohn anomaly induced by the presence of the surface metallic states, involving any of the surface phonon modes. The spin-orbit coupling results merely in a overall softening of the phonon modes of at most 6%.
In addition to the shear vertical and shear horizontal phonon densities shown in the main part of the article, Figure S2 shows the shear horizontal polarisations projected onto the first, second and third layer (SH1, SH2, SH3). If the scattering plane, defined by the incoming and scattered He beam, coincides with a mirror plane of the surface, the detection of purely SH modes is in principle forbidden due to symmetry reasons 65 and we show the calculations of the SH modes here for completeness. However, phonon modes often exhibit a mixing of polarisation components and even a purely SH mode may give rise to charge density oscillations above the first atomic layer which are eventually observed in inelastic He atom scattering.

S2. HAS DIFFRACTION SCANS
Despite the diffraction scans along ΓM which are already shown in the main part of the article, Figure S3 shows also a diffraction scans along the ΓK azimuth. Note that along ΓK there is no evidence for additional peaks close to the diffraction peaks, despite two small features close to the specular reflection which can be assigned to resonances 40 . This is indicate of the hexagonal shape of the quantum well states giving rise to a multitude of connecting vectors with similar momentum transfer between the flat sides of the hexagon and thus along the ΓM azimuth.

S3. TOF MEASUREMENTS AND SURFACE PHONON DISPERSION
The phonon energies were determined by performing time-of-flight (TOF) measurements over a wide range of incident angles between the first-order diffraction peaks and at various beam energies. The phonon dispersion was then obtained by calculating the parallel momentum transfer |∆K| for each extracted phonon energy ∆E from the conservation laws of energy and parallel momentum providing the so-called scan curve for planar scattering: 23,42 where E i is the energy of the incident beam, K i is the parallel component of the incident wavevector, and ϑ i and ϑ f are the incident and final angle, respectively. Figure S4 shows some further examples of measured TOF spectra along the ΓM azimuth and the sample at room temperature. Each TOF spectrum consists of various peaks which are located on the creation (negative x-axis, ∆E < 0) as well as the annihilation (positive x-axis, ∆E > 0) side. The peak at zero energy transfer corresponds to elastically scattered helium atoms due to the diffuse elastic peak which is caused by scattering from surface imperfections such as steps 27,28 . The scan curve associates each phonon event with a specific momentum transfer ∆K based on (Equation S4), since the incoming helium atom looses or gains momentum via inelastic scattering from a phonon. For the two upper graphs the nozzle temperature was set to 50 K and for the bottom one to 80 K while the sample was held at room temperature. The vertical red dotted line corresponds to elastically scattered atoms. The dashed line describes the scan curve which connects energy transfer with momentum transfer (y-axis on the right-hand side). The symbols denote peaks in the TOF spectrum which have been assigned to phonon events (same symbols as the main part of the article). Figure S4 shows an additional set of TOF-spectra along the ΓM direction. Most measurements covering the low energy region with the acoustic phonon modes where performed with a ≈ 10 meV beam (top and centre panel in Figure S4), while for the optical phonon modes higher incident beam energies (bottom panel in Figure S4) where used, in order to cover the optical energy region on the creation side. Note that the particularly strong peak at around 5.5 eV in the centre panel is not assigned to a phonon event. This seemingly inelastic feature originates from elastic scattering and is caused by the outer tails of the velocity distribution in the incident beam. These so-called deceptons or spurions occur within the vicinity of the diffraction peaks and give rise to a large intensity in the inelastic spectra due to their elastic nature. 66 . As mentioned in the main part of the article, there is no evidence for a strong Kohn anomaly (KA) in our mea-