Helium–Surface Interaction and Electronic Corrugation of Bi2Se3(111)

We present a study of the atom–surface interaction potential for the He–Bi2Se3(111) system. Using selective adsorption resonances, we are able to obtain the complete experimental band structure of atoms in the corrugated surface potential of the topological insulator Bi2Se3. He atom scattering spectra show several selective adsorption resonance features that are analyzed, starting with the free-atom approximation and a laterally averaged atom–surface interaction potential. Based on quantum mechanical calculations of the He–surface scattering intensities and resonance processes, we are then considering the three-dimensional atom–surface interaction potential, which is further refined to reproduce the experimental data. Following this analysis, the He–Bi2Se3(111) interaction potential is best represented by a corrugated Morse potential with a well depth of D = (6.54 ± 0.05) meV, a stiffness of κ = (0.58 ± 0.02) Å–1, and a surface electronic corrugation of (5.8 ± 0.2)% of the lattice constant. The experimental data may also be used as a challenging benchmark system to analyze the suitability of several van der Waals approaches: the He–Bi2Se3(111) interaction captures the fundamentals of weak adsorption systems where the binding is governed by long-range electronic correlations.


■ INTRODUCTION
The material class of topological insulators (TIs) has lately received broad attention 1−6 due to their protected metallic surface states and the insulating bulk electronic structure. 7,8 An archetypal TI and one of the most studied examples is the here presented Bi 2 Se 3 . 9,10 Topological surfaces show modifications of their electronic structure upon adsorption of atoms and molecules. 11−14 However, the interaction of TI surfaces with their environment, that is, atom−surface interaction potentials are barely investigated by the experiment despite the fact that topology can have implications far beyond electronic transport properties and topological materials, provides a perfect platform for studying phenomena such as heterogeneous catalysis or sensing applications. 15,16 Understanding the atom−surface interaction on topological insulators is interesting from a fundamental point of view and may help to obtain a deeper understanding of the interaction with gases and molecules in the physisorption regime. 17,18 For example, as shown recently, the long-range part of the potential can be responsible for band bending effects upon adsorption. 19 Here, we present a study of the atom−surface interaction potential for the He−Bi 2 Se 3 (111) system. The shape of the interaction potential between the surface and atoms can be extracted from atom−surface scattering experiments, and we follow this approach using helium atom scattering from the surface of Bi 2 Se 3 (111). Helium atom scattering (HAS) is strictly surface-sensitive and allows the investigation of surface structure and dynamics of conducting as well as insulating materials. 20 The technique permits measurements of the atom−surface interaction potential to a very high accuracy via selective adsorption resonances (SARs). 21,22 Studying SARs provides access to the bound-state energies ϵ n , which are supported by the atom−surface interaction potential, and more importantly, to the atom−surface interaction potential itself. Previous experimental studies of SARs have mainly been performed on salts with the NaCl structure, 20−26 while only recently, semimetals and semiconductors have been studied using this approach. 27 −29 In addition to experimental measurements of SARs, ab initio approaches have been employed to determine a numerical atom−surface interaction potential. 30−32 On the one hand, the extremely small adsorption energies of He atoms on surfaces (in the few meV region) and the delocalized nature and mobility of electrons on conducting surfaces make such systems particularly challenging even for state-of-the-art van der Waals (vdW)-corrected density functional theory (DFT) approaches. 32 On the other hand, since measurements of the atom−surface potential give insight into the atom−surface interaction dynamics within the vdW regime, experimental results can be used to test the ability of dispersion-corrected DFT approaches to simulate nonlocal interactions. 30,33,34 The atom−surface interaction potential is also a necessary ingredient for quantum mechanical calculations of elastic scattering intensities, 35−37 allowing for a comparison with experimentally observed He diffraction peak intensities. In this context, the influence of vdW forces in atom−surface scattering calculations of noble gases has recently been studied, 35,38 and the experimental diffraction intensities may even be used as a benchmark to test the performance of different vdW-corrected DFT approaches. 39 However, the comparison with diffraction intensities merely considers a small number of diffraction channels, which are accessible in the experiment, and thus, a comparison of quantum mechanical scattering calculations with SARs provides an even more rigorous test in terms of the sought atom−surface interaction potential. Following this approach, we use experimental SARs together with quantum mechanical He−surface scattering calculations to determine an accurate three-dimensional atom−surface interaction potential. Herein, experimental data from complex surfaces such as TIs may be especially valuable, in particular since it was shown only recently that, for an accurate theoretical description of the layered structure of TIs, the inclusion of vdW corrections is essential. 40

■ EXPERIMENTAL DETAILS
The experimental data in this work was obtained at the HAS apparatus in Graz, which is able to generate a nearly monochromatic beam of 4 He. The scattering geometry is defined by a fixed source-sample-detector angle with 91.5°(for a description in greater detail, see ref 41). Briefly, the He beam is generated via a supersonic expansion from a base pressure of 50 bar to 10 −6 mbar through a cooled nozzle of 10 μm. The central part is selected with a skimmer (310 μm orifice), creating a He beam with an energy spread of ΔE/E ≈ 2%. By varying the nozzle temperature, the beam energy can be tuned between 9 and 20 meV. The beam then hits the sample in the main chamber under ultrahigh-vacuum (UHV) conditions (p ≤ 2 × 10 −10 mbar) and is further detected using a quadrupole mass analyzer. For varying the incident angle ϑ i in the fixed source-sample-detector geometry, the sample can be rotated.
For a detailed description of the sample growth procedure, see the study of Bianchi et al. 42 After in situ cleavage of the sample in a load-lock chamber, 43 the cleanness and purity of the sample can be further studied using low-energy electron diffraction (LEED) and Auger electron spectroscopy (AES). The rhombohedral crystal structure of Bi 2 Se 3 is built of quintuple layers (QL), which are bound to each other through weak van der Waals forces. 42 The unit cell consists of three QLs and shows Se termination upon cleavage. The surface along the (111) cleavage plane has a lattice constant of a = 4.14 Å at room temperature. 44 The sample was fixed on a sample holder using thermally conductive epoxy.
The intensity of the specular reflection throughout the measurements typically reached values with a signal-to-noise ratio of 10 3 above the diffuse elastic background (Figure 3), while the full width at half maximum (FWHM) was typically about 0.015 Å −1 . Hence, the angular broadening of the specular peak is mainly limited by the angular broadening of the apparatus, giving rise to an estimate (lower limit) for the quality of the crystal 45 with domain sizes larger than 1000 Å.
The sample temperature was varied between cryogenic temperatures (113 K, via a thermal connection to a liquid nitrogen reservoir) and room temperature (300 K). For the first characterization, the scattered specular intensity was measured in dependence of the sample temperature to determine the surface Debye temperature (Θ D = 122 K 46 ). All scattering calculations presented in this work have been corrected by the corresponding Debye−Waller attenuation based on the experimentally determined Debye temperature.
After the first characterization of the crystal, various elastic diffraction scans at 113 K and room temperature were collected in both high-symmetry directions and at different incident energies. Furthermore, measurements of the specular reflection in dependence of the beam energy (measured at 113 K) allow us to obtain further details of the SARs and the atom−surface interaction potential (see Refinement of the Interaction Potential).

■ RESULTS AND DISCUSSION
After an introduction to the atom−surface interaction potential and the kinematic analysis, we begin our analysis of the SAR features following the free-atom model. Based on the approximate surface-averaged potential, we are then determining the corrugation amplitude of the potential from diffraction measurements to acquire a three-dimensional potential. Finally, we are going to compare the experimental SARs with quantum mechanical scattering calculations to further refine the three-dimensional atom−surface interaction potential.
Atom−Surface Interaction Potential. As stated by Bragg's law, if an atom is scattered by a periodic surface, the change of the wavevector component parallel to the surface, K, must be equal to a surface reciprocal lattice vector, G.
Here, we present mainly measurements where the polar (incident) angle ϑ i is varied around the corresponding axis, while the scattered beam intensity is detected. For elastic scattering, the momentum transfer parallel to the surface, ΔK, is then given by where k i is the incident wavevector, and ϑ i and ϑ f are the incident and final angles with respect to the surface normal, respectively. Selective adsorption phenomena, which may appear upon scattering of atoms from a periodic surface, occur due to the attractive part of the atom−surface interaction potential. For an elastic process, the kinematically allowed G-vectors are those for which the wavevector component perpendicular to the surface k f, z 2 is positive. If an incident He atom hits the surface, it can undergo a transition into a bound state on the surface with −|ϵ n |. The process happens while the He atom is diffracted into a channel, which is kinematically disallowed (k f, z 2 < 0). Such SARs can only happen if the difference between the energy of the incident atom and the kinetic energy of the atom moving parallel to the surface matches the binding energy ϵ n of the adsorbed atom 23 where m is the He mass. From eq 2, it is clear that studying SARs provides access to the bound-state energies ϵ n (K i , G) and, more importantly, to the atom−surface interaction potential.
The Journal of Physical Chemistry C

Article
In this work, we have analyzed SARs upon scattering of He from Bi 2 Se 3 using a corrugated Morse potential (CMP). Despite the deviation of the CMP from the expected z −3 asymptotic behavior, it has been shown that the overall shape of the CMP represents the measured bound states well enough. A comparative study of various potentials with different asymptotic behavior has shown a similar outcome when subsequently used in close-coupling calculations. 37,47 Likewise, the validity of this potential for TI surfaces 48 and other layered materials such as transition metal dichalcogenides 49 has been proven. In addition, the CMP greatly simplifies the treatment of several steps within the closecoupling (CC) algorithm, allowing for an analytical solution in those cases, which leads to a reduced computational cost.
The three-dimensional CMP, written in dependence of the lateral position R on the surface and the distance z with respect to the surface, is 50 Ä with κ being the stiffness parameter, D being the depth of the potential well, and ν 0,0 being the surface average of the exponent of the corrugation function. The electronic surface corrugation is described by ξ(R), where R describes a periodically modulated surface corresponding to a constant total electron density. For ξ(R), a two-parameter Fourier ansatz was used, which is described by a summation of cosine terms. The Fourier series expansion is based on the sixfold symmetry of the topmost layer of the surface (see eq 14 in the Appendix), which is the only relevant layer when considering the energies of the impinging He atoms used in this study. 29,37 The corrugation magnitude is then typically given in terms of the peak-to-peak value ξ pp of ξ(R).
The laterally averaged surface potential V 0 (i.e., without corrugation) of eq 3 is given via The bound states can be described analytically by To determine the three-dimensional potential, we start first with the laterally averaged atom−surface interaction potential, trying to identify SARs in various diffraction scans. After determining a laterally averaged atom−surface interaction potential based on the SAR positions and the free-atom model, we determine the corrugation amplitude ξ 0 of the potential by comparison of the diffraction intensities based on closecoupled calculations with the experimentally determined diffraction intensities. We are then further refining the potential using close-coupled calculations of the resonance positions, optimizing the agreement with the experimental measurements of SARs. The whole procedure is described in more detail below.
Determination of SARs in the Elastic Scans. In a first step of determining the atom−surface interaction potential, the laterally averaged atom−surface interaction potential is used to identify the SARs in various elastic scans. When performing a diffraction scan according to eq 1, the kinematic condition (eq 2) can be fulfilled for specific values of ϑ i . Two typical diffraction scans for the high-symmetry direction ΓM are shown in Figure 1. The x-axis has been transformed to parallel momentum transfer using eq 1. In addition to the diffraction peaks, smaller features caused by SARs can be seen.
First, we will restrict the kinematic condition given in eq 2 to the free-atom approximation, which assumes that the surface potential is adequately described by the laterally averaged interaction potential (eq 4). Strictly speaking, it holds only in the case where the surface corrugation approaches zero, which is not possible in reality, since corrugation is necessary to provide the G-vector for the resonance processes. We will then use and refine the full three-dimensional atom−surface interaction potential at a later point within the formalism of close-coupled calculations.
Nevertheless, the free-atom approximation is useful for a first identification of SARs as it treats the binding energies ϵ n (K i , G) in eq 2 as constants and is therefore independent of K i and G. The introduction of a corrugation may give rise to changes in the resonance positions, which tend to become more important with increasing surface electronic corrugation. Since the peak-to-peak corrugation of similar materials is in the region of 5 − 9% of the lattice constant, 37,45 one needs to keep in mind that considerable shifts of the resonance positions may occur.
Using the free-atom model as a starting point, the position of the SARs (eq 2) can be rewritten in terms of the incident angle The intensity scale has been expanded, and in addition to the diffraction peaks that are cut off due to their high intensity, small peaks and dips corresponding to selective adsorption processes can be seen. The vertical lines in various colors illustrate the kinematic conditions for five bound-state energies ϵ 0 −ϵ 4 according to eq 6. Each line is labeled with the Miller indices of the associated G-vector for the particular resonance condition.
The Journal of Physical Chemistry C Article ϑ i and the incident wavevector k i . The latter is given by the beam energy, which is determined by the nozzle temperature. The scattering vector G can be separated into two parts G = (G ∥ , G ⊥ ) normal and parallel to the incidence plane The peaks or dips at a particular incident angle ϑ i in the scans can be associated with a certain bound state energy ϵ n and scattering channel using eq 6. The position of several SARs is indicated by vertical lines (based on solving eq 6 for ϑ i ) in Figure 1 and labeled according to their diffraction channel G, while the different colors correspond to different bound states ϵ n . For better visibility of the SARs, the sample was cooled down to 113 K to minimize the inelastic background as well as the linewidth of the resonances. The scan in the upper panel was measured with an incident beam energy of 11.25 meV and the lower panel at 14.26 meV. The vertical lines are according to the first five bound-state energies ϵ 0 −ϵ 4 of the later determined laterally averaged interaction potential defined by their color, while the numbers next to the lines denote the corresponding reciprocal lattice vectors.
It becomes immediately evident that the assignment of resonances in a single diffraction scan such as Figure 1 is quite ambiguous since numerous combinations of ϵ n and G can be thought of to fulfill eq 6. For a better assignment of the bound states and G-vectors, a collection of numerous scans was put together, as shown in the contour plot of Figure 2. To monitor the drift of the resonances versus the incident beam energy, many diffraction scans with various E i were performed. Putting this collection of measurements together allows us to identify the bound-state resonances much easier. In particular, the curvature of a certain resonance feature allows us to determine the corresponding G-vector, which greatly simplifies the search for the associated bound state ϵ n , which fulfills eq 6.
A total of 27 ϑ-scans at incident energies ranging from 9.5 to 14.5 meV was collected to construct the contour plot. The xaxis in Figure 2 corresponds to the parallel momentum transfer and the z-axis to the scattered intensity. The y-axis is formed by plotting the scans with various incident energies; that is, a cut at y = constant would result in a graph such as Figure 1. The plot was constructed by connecting the individual scans on a two-dimensional grid using a linear interpolation, while the intensity (z-axis) is plotted on a logarithmic scale. Scattered intensities (mainly the diffraction peaks) that exceed a certain value have been cut off to increase the visibility of the resonance effects. The superimposed solid lines correspond to identified bound states according to the free-atom approximation (eq 2), with the color coding and associated G-vectors labeled in the same way as in Figure 1. A number of lines of high and low intensities, which we identify as selective adsorption resonances features, can be seen to run across the data set.
Following the described approach in analyzing the position of these SARs, a set of five distinct eigenvalues of the laterally averaged potential is obtained ( Table 1). The experimentally found bound states are then used to determine the laterally averaged potential (eq 5) by minimizing σ ϵ N 1 n N n n 0 with N being the number of bound states included. In doing so, a laterally averaged potential with the parameters D = (6.6 ± 0.2) meV and κ = (0.58 ± 0.07) Å −1 is determined. The obtained potential supports a total of seven bound states, with ϵ 5 and ϵ 6 being quite close to zero, that is, to the threshold condition.
Note that not all lines of high and low intensities can be explained using SARs based on the free-atom approximation. We will later see that using a three-dimensional potential will give rise to shifts compared to the free-atom approximation, while considering inelastic channels can give rise to changes from maxima to minima and vice versa. 51 There are also some features present that cannot be explained by SARs even when considering a full three-dimensional potential. The features we are referring to show only a weak dependence of ΔK with respect to the incident energy E i , and thus, no associated Gvector can be found that would explain such a curvature. One of these features, indicated by the black arrow labeled with KF The black arrow labeled with KF indicates a region of increased intensity, which is likely to correspond to kinematic focusing (KF) rather than a resonance effect. a The corresponding internal linewidths Δϵ n of the bound states (based on the experimental width of the resonances) and their lifetimes τ n are also given.
The Journal of Physical Chemistry C Article in Figure 2, is likely to correspond to a kinematic focusing (KF) effect (see ref 29). Moreover, there appear to be features with increased intensity next to the first-order diffraction peaks (see, e.g., the top panel in Figure 1). Additional features in diffraction scans due to spin-conserving electronic interpocket transitions have been observed for the semimetal Sb(111), 52 and the observability of similar transitions in the TI Sb 2 Te 3 has been suggested by theory. 53 The observation of additional dispersion curves in Bi 2 Se 3 54 in analogy to Sb(111) suggests a similar assignment of the additional peaks in the diffraction spectra. On the other hand, the surface electronic structures of both Sb(111) and Sb 2 Te 3 exhibit narrow electron pockets, while the situation for Bi 2 Se 3 is different, with a single Dirac cone close to Γ̅ that typically evolves into multiple states due to the formation of quantum well states. For the latter case, it is still under debate whether both storage in the UHV chamber and exposure to intense ultraviolet light are necessary ingredients. 42,55 The distance of the above-mentioned features in terms of ΔK with respect to the diffraction peaks would make an electronic transition induced by the helium atom possible (via scattering from −k to +k). However, the question remains whether this would be a spin-conserving transition, for example, between the Rashba spin-split quantum well states, or whether, for example, a phonon would be required in the process to allow for a spin-flip. In any case, since we can exclude the fact that these features are related to SARs, their origin is not relevant for a determination of the atom−surface interaction potential. Hence, the assignment and discussion about these features are beyond the scope of this work and will be treated separately. 54 Calculation of the Scattering Intensities. Once a first estimate of the laterally averaged atom−surface interaction potential has been established, we are going to consider the three-dimensional CMP (eq 3). A comparison of quantum mechanical calculations of the scattered intensities with the experimentally found ones yields a value for the corrugation ξ pp . The process of elastically scattering a He atom from a surface can be described by the time-independent Schrodinger equation together with the CMP. In the exponential term of the potential (eq 3), the Fourier series expansion (see eq 14) is used, which yields a set of coupled equations for the outgoing waves. These waves are solved using the close-coupling algorithm for a finite set of closed channels 28,56 (see Details on the Close-Coupling Calculations). The corresponding coupling terms for the CMP can be found in several references. 29,37,47 To determine the corrugation of the sample surface, the elastic peak intensities are calculated with the close-coupling algorithm and compared with the experimentally measured values of the peak areas. 45 The calculated values of the diffraction peak intensities are corrected with the Debye− Waller attenuation (Θ D = 122 K 46 ), while the measured intensities are determined from the peak area using a fitted Voigt profile. 37 The reason for using the measured peak areas rather than the peak heights is the need to account for the broadening due to the energy spread of the He beam and additional broadening of the diffraction peaks caused by the angular resolution of the apparatus as well as domain size effects of the crystal surface. While the values for D and κ of the potential were held constant, the corrugation ξ pp was varied between 0.01 and 0.6 Å with a step width of 0.001 Å. The corrugation amplitude ξ pp , which describes the best corre-spondence between measured (I G exp ) and calculated diffraction intensities (I G sim ), is found by minimizing a measure of the deviation R 8) with N being the number of experimentally measured diffraction peaks. 20 For the analysis, we used a total of 35 elastic scans at various incident energies, trying to minimize eq 8, that is, searching for the global minimum of R in dependence of ξ pp . In Figure 3, two of these diffraction scans along the two high-symmetry directions of Bi 2 Se 3 (111) are shown together with the measured and calculated peak areas (symbols in the figure).
The equipotential surface describing the electronic corrugation ξ(R) corresponds to the classical turning points of the potential, that is, the locus of all points for which V(R, z) = E i holds. Hence, one may expect that, with increasing beam energy, the turning point shifts to distances closer to the ion cores and the scattered He atoms experience a larger corrugation amplitude ξ pp . The corrugation then shows a dependence on the incident energy of the molecular beam, typically following a monotonic increase with E i . 57 Therefore, the above-described optimization routine was repeated by taking into account several experimental spectra recorded around three specific incident energies E i , as shown in the inset of Figure 3. Indeed, when only diffraction scans taken at the lowest incident energy are considered in the analysis explained above, the best fit value of ξ pp decreases. For medium incident energies (E i = 10.3 meV), the minimum of R is found with ξ pp = 0.21 Å, while for higher incident energies (E i = 14.4 meV), it increases to ξ pp = 0.27 Å. Optimizing eq 8 based on all recorded diffraction spectra yields ξ pp = 0.25 Å, which is thus an average corrugation over the whole range of beam energies considered in this study. We will later observe (see Refinement The Journal of Physical Chemistry C Article of the Interaction Potential) that, for the SARs, small changes in terms of ξ pp are much less important compared to changes of the well depth D and the stiffness κ.
As mentioned above, the free-atom approximation used in the first approach (see Determination of SARs in the Elastic Scans) of determining a potential neglects the corrugation of the surface, meaning that the coupling term vanishes and the band structure would solely consist of parabolic bands. In considering surfaces with larger corrugations, higher-order Fourier components of the surface potential have to be considered. The three-dimensional corrugated surface potential (eq 3) may then give rise to substantial deviations from the free-atom parabolic bands. Hence, we will use again quantum mechanical calculations to accurately describe the resonance positions and, in doing so, refine the full three-dimensional potential in the following. We will also rerun the optimization routine for the electronic corrugation ξ pp , which was, at first, determined using D and κ as obtained from the averaged potential, with the refined three-dimensional potential.
Refinement of the Interaction Potential. In addition to elastic ϑ-scans, the intensity of the specular peak can be measured as a function of the incident wavevector k i . Such a so-called drift spectrum shows again SARs at places where the kinematic condition (eq 2) is fulfilled. As can be seen in Figure  2, the SAR conditions "move through" the specular condition and cause the intensity to change. However, in contrast to the ϑ-scans where most SARs appear as weak features in between the diffraction peaks, a measurement, while staying at the specular reflection, allows for the highest signal-to-noise ratio in the experiment. The top panel in Figure 4 shows such a drift scan along the ΓK azimuth where a multitude of broad and narrow peaks and dips is visible. Hence, such a measurement gives access to the detailed shape of the potential, and by comparison with calculations, we will further refine the potential, looking into shifts with respect to the free-atom approximation.
At the same time, additional effects may complicate the analysis of SARs in a drift scan. First, the incident beam intensity decreases with increasing nozzle temperature T N since the He flow through the nozzle is proportional to T 1/ N . The effect can however be easily accounted for by correcting the scattered He intensity in the measurement with the corresponding factor (as done in the top panel of Figure  4). Second, surface imperfections such as terraces and steps can give rise to variations of the scattered intensity: The interference of outgoing waves, which are scattered from different terraces, can cause periodic oscillations of the detected signal as a function of k i . 20,29 From the major peaks and dips of these oscillations, the terrace height(s) of the investigated sample can be calculated. 20,29 These peaks are typically much broader than the peaks and dips caused by SARs, and we will attempt to analyze those in the Appendix.
Moreover, the energy distribution of the incident beam will give rise to a broadening of the natural linewidths of SARs, which can be incorporated by numerically convoluting the elastic intensity with the appropriate distribution in incident energy (see the turquoise dash-dotted curve in the bottom panel of Figure 4 # ). Finally, SARs tend to become less pronounced and broader with increasing sample temperature due to a linewidth broadening and the increasing importance of inelastic effects, as can be seen when comparing Figure 4 with Figure 9.
We turn now to the refinement of the interaction potential based on the drift scan (top panel of Figure 4), which shows the measured specular intensity as a function of the incident wavevector k i . The sample was aligned along the ΓK azimuth and held at a sample temperature of 113 K while changing the nozzle temperature from 44 to 100 K. The colored lines display the SAR positions according to the free-atom approximation using the optimized surface potential. The numbers next to the lines denote the corresponding reciprocal lattice vectors following the same nomenclature as for Figure 1. In cases where several G-vectors lead to the same solution of the kinematic equation, the Miller indices of only one G-vector are given. In addition, the vertical dashed lines correspond to the threshold energies of the surface potential. In principle, threshold resonances can also give rise to intensity variations, which have been predicted to be experimentally detectable for scattering of atoms from highly corrugated surfaces. 56, 58 On the other hand, threshold resonances have only been observed experimentally for scattering He from a ruled grating upon grazing incidence. 59 The measured drift spectrum can be simulated using calculations based on the elastic close-coupling formalism. The simulated spectra can then be compared to the SAR positions in the experimental data, and in doing so, the surface potential can be further refined. For these calculations, the corrugation values from above were used, while the values of D and κ were varied in the neighborhood of the first estimated values. After the potential parameters D and κ, which describe the closest agreement with the measurements, have been found, the corrugation ξ pp is further refined by minimizing eq  Figure 4, the result of the simulation using the close-coupling formalism is plotted, which has further been multiplied with the corresponding Debye−Waller factor (see Close-Coupling Formalism for further details about the calculations).
Note that, due to the above-described additional effects, it is not possible to resemble the actual shape of the whole measured drift spectrum using the elastic close-coupling simulation since, at finite temperature, the linewidth and shape of the SARs will be influenced by inelastic channels, while at the same time, the oscillations caused by terraces are superimposed onto the experimental spectrum. From a theoretical point of view, these effects have been mainly considered in the limit of low corrugated surfaces, 60−62 but it is well known that inelastic events are expected to account for the attenuation of the line shapes as observed in the experiment and can even turn maxima into minima. 51,61 In the elastic theory of resonant scattering, it has been shown that the occurrence of minima, maxima, and mixed extrema can be explained and predicted by establishing some general rules. However, in many cases, their applicability is limited since these were derived for weak coupling conditions and hard model potentials. 60 Thus, we will mainly concentrate our analysis on the position of the SARs in terms of k i , as described below. Despite the complications caused by linewidth broadening and the additional oscillations caused by the terraces, most peaks and dips of the measured data can be identified in the simulated drift scan. The position of the SARs based on the free-atom approximation, denoted by the colored lines, are shifted with respect to the peaks of the quantum mechanical calculations using the full three-dimensional potential. Generally, resonances with higher corresponding G-vectors tend to show a larger shift compared to those with lower indices (see Close-Coupled Calculations of the Drift Scan for a set of simulated curves showing how the SARs change upon variation of D, κ, and ξ pp ). However, it is hard to assign each peak in the calculation to the corresponding resonance condition of the free-atom approximation.
Instead of defining a global χ 2 parameter for the goodness of fit that would not distinguish between different aspects, we will concentrate on optimizing the position of a number of specific features. Due to the above-described broadening effects and superimposed oscillations from the terraces, it is impossible to define a simple parameter that adequately describes the overall agreement between the experiment and simulation. An approach that considers a global χ 2 parameter will be highly sensitive to unimportant aspects such as the energy broadening of the incident beam or the described oscillations due to the terraces that modulate the intensities. Hence, we concentrate on optimizing the position of a number of specific "target" features, which can be identified directly in the experimental measurement and are then checked qualitatively to give an improvement over the entire data set. Target features should be representative of the entire potential, that is, being associated with different G-vectors and bound states. Furthermore, for practical purposes, we choose features that can be unambiguously identified, clearly above the background and preferably with no other resonances close enough to cause confusion with a different SAR channel or an intersection with another feature. Finally, SARs leading to peaks are generally preferred to those leading to dips. 21 We concentrate on testing the positions k i of the SARs in the simulated data compared to the peak positions in the experimental spectrum where k i exp are the experimental data values with uncertainties σ i (σ i ≈ 0.02 Å −1 ), and k i sim are the calculated peak positions. Figure 5 shows the effect on R p 2 if the potential well depth D and stiffness κ are varied around the optimized parameters. The whole set of simulated drift spectra for this purpose can be found in Close-Coupled Calculations of the Drift Scan. Note that the effect of varying ξ pp (Figure 8) on the resonance positions is rather small, with the main effect being the changes in the shape and amplitude of the resonances. Thus, ξ pp is better determined in comparison with the diffraction peak intensities, and we concentrate on a refinement of D and κ based on the position of the SARs in the drift scan. Figure 5 clearly illustrates that the final optimized values of D and κ correspond to a minimum in terms of R p 2 . Another possibility for a measure of the agreement between the simulation and experiment is to try and adopt an R-factor in analogy to what has been used in LEED in experiments. According to Pendry, 63 a good R-factor for comparison with calculated LEED experiments should be sensitive to the peak position. It should not be sensitive to the absolute intensities but take into account the relative intensities of features that are close in energy. Based on those criteria, an R-factor for SARs can be defined by considering a measure of the curvature, that is, the second derivative of the intensity I exp ″ in the drift scan according to where the scaling constant c = ∫ I exp d k/ ∫ I sim d k is used to normalize the experimental and simulated curves with respect to each other. 63,64 Unfortunately, the number of additional effects in the experimental spectrum makes such an approach quite difficult. We can restrict the approach to certain regions in the calculated and simulated spectra, for example, considering the region between k i = 5.1 and 5.4 Å −1 . We see in Figure 5 that the overall trend is the same compared to using R p as a measure, although the latter seems to be more robust in the present case. Compared to the results from Determination of SARs in the Elastic Scans, the well depth D decreased, while the stiffness κ showed only a subtle change. By using this optimization process, we can define a measure for the agreement between the close-coupled calculations and the experiment. With respect to the first approach based on the free-atom approximation, the uncertainties of all values are significantly reduced.
In comparison to previously studied systems, the value of the well depth D is between those found for Bi 2 Te 3 (111) (6.22 meV) 48 and Bi(111) (7.9 meV), 37 while the stiffness κ of the He−Bi 2 Se 3 (111) potential is significantly smaller than for Bi 2 Te 3 (111) and for the Bi(111) single crystal, although larger compared to Sb(111). 47 Because Bi 2 Se 3 has a smaller vdW gap than Bi 2 Te 3 and Sb 2 Te 3 , 40,42,65 a different interaction of the outermost layer in the weak vdW regime might be expected, and indeed, this is confirmed by the determined potential in terms of the different stiffness κ. Due to the polarizability of both Bi and Se, one would expect a "soft" potential with a significant long-range attractive part. On the other hand, the determined potential parameters indicate that the potential is actually "stiffer" than the typical potential of a simple flat metallic surface. Detailed simulations will be required to resolve the importance of the above-mentioned effects, and the presented data may provide a benchmark to test challenging vdW approaches in DFT.
Finally, revisiting the calculation of the diffraction peak intensities with the refined potential parameters yields a peakto-peak corrugation for an incident beam energy E i = 11.7 meV. The average overall beam energies considered in this study correspond to a surface electronic peak-to-peak corrugation of (5.8 ± 0.2)% of the surface lattice constant. The value is similar to the semimetal Bi(111) (5%) 37 while being smaller than the reported 9.6% for Bi 2 Te 3 . 48 Since the spin−orbit coupling in Bi 2 Se 3 is stronger than in Bi 2 Te 3 , giving rise to a different electronic structure of the topological surface states, 9,66 this may also cause a different electronic surface corrugation, although it is of course difficult to make any connections between the localized electronic bands in terms of k-space and the "global" surface electronic corrugation. At the same time, the work functions of Bi 2 Se 3 and Bi 2 Te 3 are quite similar, 45,67 and hence one might not expect large differences in terms of the surface electronic corrugation. Nevertheless, we can certainly conclude that the surface electronic corrugation of both Bi 2 Se 3 and Bi 2 Te 3 is of the same order of magnitude and the charge smoothing due to the Smoluchowski effect is definitely less pronounced for both TIs compared to the observations of flat metal surfaces. 20 Linewidth and Lifetime of SARs. The angular broadening of SARs in experimental measurements is related to the lifetime of the corresponding bound state, that is, the time that the He atom spends in the bound state before it leaves the surface. The experimentally determined (external) width has to be corrected for resolution aspects of the apparatus to determine the internal or natural linewidth of the bound state. 56, 68 The external full width at half maximum (FWHM) was determined by fitting the SAR features in several angular diffraction scans with a Gaussian function. The external width Δϑ total is corrected for the resolution aspects based on the convolution of two Gaussian distributions, yielding an internal angular width (Δϑ res ) 2 = (Δϑ total ) 2 − (Δϑ app ) 2 . Δϑ app accounts for the angular resolution of the apparatus as well as the energy spread in the beam and was obtained from the FWHM of the associated diffraction channel. The natural linewidth Δϵ n of the bound state in terms of energy is then given via 68 where G ∥ is the parallel component of the G-vector associated with the resonance. Finally, the lifetime τ of the bound states can be established from the uncertainty principle using τ n = ℏ/ Δϵ n . 56 The natural linewidths and corresponding lifetimes of all bound states are listed in Table 1. The lifetimes increase with increasing quantum number n of the bound state. The trend is expected since bound states with a higher n are further away from the surface and thus experience less of the surface corrugation, decreasing also the probability of scattering events that would cause the He atom to leave the surface.
In general, various factors will affect and limit the lifetime of a bound state. For elastic scattering, the lifetime is limited by the probability of being scattered out of the bound-state channel, which relates to the form and amplitude of the lateral corrugation in the atom−surface potential. Considering inelastic processes, the natural lifetime will be further reduced by factors such as defect or phonon scattering, with the latter becoming more important with increasing temperature. 27,56 Due to the required high experimental resolution, information about the linewidth and lifetime of SARs is limited to a very small number of systems. 27,47,68,69 Direct experimental information is only available for the He− LiF(001) system 68 and the He−Sb(111) system. 47 The internal linewidths of the He−LiF(001) system are comparable to those found for He−Bi 2 Se 3 in Table 1. Based on elastic scattering events and the similarity of the atom−surface interaction potentials (similar depth although larger corrugation in the case of He−LiF(001)), one would expect similar linewidths in both cases, and indeed, this is confirmed by the measurements. On the one hand, with increasing surface temperature (room temperature in ref 68 compared to 113 K in our study), inelastic events may become more important; on the other hand, the higher Debye temperature of LiF compared to Bi 2 Se 3 suggests that phonon scattering will again be similar when comparing both studies.
The influence of different atom−surface interaction potentials on the linewidth has also been subject to previous studies. As noted by Tuddenham et al., 27 based on closecoupled calculations, a Morse potential with the same corrugation as a corresponding reference potential gives features whose linewidth is similar to those seen in the experiment. Hence, we hope that the experimental determination of the linewidths presented in this study will initiate further work in this direction and, for example, in comparison with inelastic close-coupled calculations, eventually allow to The Journal of Physical Chemistry C Article rule out whether elastic or inelastic scattering channels are mainly responsible for the lifetime of bound states.

■ SUMMARY AND CONCLUSIONS
In summary, we have determined an atom−surface interaction potential for the He−Bi 2 Se 3 (111) system by analyzing selective adsorption resonances in helium atom scattering spectra. For a first approximation, we start with the free-atom approximation and a laterally averaged atom−surface interaction potential, which is then further improved and refined based on closecoupled calculations to obtain an accurate three-dimensional atom−surface interaction potential. The free-atom approximation cannot provide information about the shape of the resonances, and by comparison with close-coupled calculations, we are able to obtain the complete experimental band structure of atoms in the corrugated surface potential. Following a systematic analysis, the He−Bi 2 Se 3 (111) potential is best represented by a corrugated Morse potential, which exhibits a well depth of D = (6.54 ± 0.05) meV and a stiffness of κ = (0.58 ± 0.02) Å −1 . The surface electronic corrugation varies slightly depending on the incident beam energy with an average of (5.8 ± 0.2)% of the lattice constant.
Inelastic processes and phonon-mediated resonances have been proven to play important roles, and a precise atom− surface interaction potential as determined in this study is a necessary ingredient to investigate the effects such as the temperature dependence and linewidth of selective adsorption resonances. From the angular width of selective adsorption resonances in the scattering spectra, we are able to obtain the natural linewidth of the resonances and an estimate for the lifetime of the bound states. Moreover, since a meV He beam is scattered in the low-density region dominated by the tails of Fermi-level surface states, studying selective adsorption resonances provides access to the interaction of TI surfaces within the weak adsorption regime. Hence, we hope that the present data will encourage future ab initio studies to test the ability of vdW corrections on the current system.
As a side note, in the Appendix, we use intensity oscillations due to the interference of the He beam being scattered from different terraces to analyze the step heights of the cleaved crystal surface. The analysis confirms the existence of steps with a quintuple layer height with an indication that subquintuple layer steps may exist also. Despite the existence of terraces, the angular broadening in the diffraction spectra speaks for the high quality of the cleaved sample with domain sizes larger than 1000 Å.

■ APPENDIX Details on the Close-Coupling Calculations
Close-Coupling Formalism. Considering the periodicity of the surface lattice, both the surface potential V(r) and wave function ψ(r) can be written as a Fourier series Inserting the latter into the time-independent Schrodinger equation yields 13) for the z-direction, where k G, z 2 is the z-component of the particles kinetic energy after the surface interaction, and V 0 is the laterally averaged interaction potential.
Equation 13 is a set of coupled equations, with V G − G′ being the coupling terms, which can be written in a matrix form, since the close-coupling method treats each ψ G as a scattering channel. Finally, the matrix equation is solved by using a Numerov algorithm. 29,56 The two-parameter Fourier ansatz used for the corrugation where x and y depict the coordinates, and ξ 0 determines the corrugation amplitude. The value of the corrugation is determined by the peak-to-peak corrugation ξ pp of eq 14.
For all close-coupled calculations presented in this work, all open channels were considered in the simulations. The number of closed channels was adjusted until the error was determined to be <10 −4 . In the drift simulations, a value of less than 100 closed channels was sufficient, while for the calculation of the diffraction intensities, a minimum of 150 channels was required for the result to converge. For the drift spectrum calculated in Figure 3, 100 closed channels were used, and the integration boundaries in terms of z were set to [−9, 20] Å with an energy resolution of dE = 0.001 meV.
As mentioned in the main part of the paper, the value of the corrugation ξ pp exhibits a feeble dependence on the incident energy E i . Therefore, the best fit value of ξ pp was determined for several sets of diffraction spectra, with each set taken at a different incident energy. In a first approximation, we assume that the corrugation increases linearly within the considered energy region, according to  Figures 6 and 7 are used in Refinement of the Interaction Potential for the refinement of the potential and the calculation of the reliability factors. For these simulations, the corrugation magnitude was varied according to the incident energy, following eq 15.
In Figure 6, the well depth D is varied within a range of 6.3 to 6.8 meV, which clearly shows different shifts of the individual features in the drift spectrum depending on D. With increasing well depth, the peaks wander toward lower incident energies E i . In addition, it becomes clear that certain features shift at various "speeds" upon varying the well depth. For instance, the position of the sharp feature at E i ≈ 5.6 meV barely changes along the calculated range of D. Therefore, this feature might be caused by a threshold condition. In general peaks originating from lower index, G-vectors show a slower shift, as already mentioned in the main part of the paper.

Article
In the next step, the well depth was held constant, while the stiffness κ was varied in a range from 0.54 to 0.62 Å −1 . The peaks in Figure 7 tend to shift to higher incident energies with increasing stiffness. The most prominent feature changes its position rapidly from 4.7 to 5.1 meV over the whole calculated range of κ.
Finally, the corrugation ξ pp was varied from 0.14 to 0.32 Å. In Figure 8, it becomes clear that the influence of the corrugation within the considered range of ξ pp is much more subtle in terms of the position of the peaks. However, ξ pp significantly changes the shape and intensity of the individual features. For example, the peak at E i ≈ 4.9 meV changes its shape from a minimum to a pronounced maximum when increasing ξ pp . In addition, some features might not even be visible at lower values of the corrugation, changing to sharp maxima at higher values of ξ pp .
Detailed Analysis of the In-Phase/Anti-Phase Conditions due to Steps If the specular intensity is dominated by interference of waves scattered from different terraces, the drift spectrum can also be used to determine the terrace height(s) of the investigated sample surface. 20,29,70 When changing the incident wavevector k i , the specular intensity shows oscillations due to constructive (in-phase) and destructive (antiphase) interference of the outgoing waves. The phase difference φ upon scattering from two terraces separated by a step height h is given by 29,70 hk hk h k 2 cos 2 where Δk z is the momentum transfer perpendicular to the surface. In Figure 9, the specular intensity is plotted as a function of the incident wavevector k i for both high-symmetry directions ΓM and ΓK with the sample held at room temperature. The measured intensity has again been corrected for the decreasing He flux with increasing nozzle temperature T N . 29 The variation of the intensity can then be described by a simple model that overlaps the plane waves that are reflected of the different terraces from the sample. The intensity as a function of Δk z is described as     17) for N different step heights h j with a probability a j . I 0 is the intensity of an ideal surface without any steps, and the exponential prefactor accounts for the Debye−Waller attenuation, with W being the Debye−Waller factor. 29,70 The latter considers that the Debye−Waller factor changes with increasing perpendicular incident momentum k iz 29 during a drift scan.
The individual QLs in Bi 2 Se 3 as well as other layered TIs are bound to each other through weak vdW forces, which gives easy access to the (111) surface by cleavage. 42,48 Hence, steps that occur upon cleaving are expected to correspond to the QL distance along the c-axis of the conventional hexagonal unit cell (9.6 Å at room temperature 44 ). Indeed, the main oscillation in Figure 9 suggests a step height of 9.9 Å (red dash-dotted curve according to eq 9) in good agreement with the separation between the Se termination layers of a QL. 71 The small deviation from the theoretical QL distance may be due to electron spill-out, that is, the fact that the electronic step height observed in HAS does not necessarily correspond to the crystallographic one as measured, for example, with LEED. 72 On the other hand, it appears in Figure 9 that an additional oscillation, giving rise to a further modulation of the intensities from the QL steps and according to a smaller step height, is present also. Sub-QL steps have previously been observed during thin-film growth of TIs, 73,74 while different terminations have also been found for cleaved single crystals depending on the sample treatment. 75,76 Inclusion of a second step height (h 2 = 3.6 Å, dashed curves in Figure 9) seems to improve the agreement with the experimental data, which suggests that terraces with the internal Se layer exist in some cases. 76 Here, considering a different number of steps along ΓK and ΓM gives a better agreement with the experimental data, which would imply that the orientation of the terraces is not completely random and depends on the crystal orientation. The fact that Bi 2 Se 3 has a smaller vdW gap than Bi 2 Te 3 and Sb 2 Te 3 might explain the existence of sub-QL steps of the cleaved samples.
Note, however, that the number of parameters in eq 17 and the overlapping with SARs makes an exact determination of the step height difficult. While some of the peaks appear in the ΓK direction as well as in the ΓM direction, the peak at around k i = 4.9 Å −1 can only be found in the ΓK direction. It should also be considered that SARs with lower interacting G-vector "move" slower across the specular peak, which may lead to peaks that look like oscillations caused by the steps. From the above analysis, it is clear that oscillations due the QL spacing are present, while there is also a hint toward sub-QL steps. At the same time, we would like to stress that, during such a drift scan, HAS is particularly sensitive to single defects, and hence the appearance of sub-QL steps is likely to be negligible. In fact, as mentioned in Experimental Details, the small width of the specular and diffraction peaks in the elastic scans suggests a very high quality of the cleaved crystals.
Finally, as already mentioned above, the measurements at room temperature ( Figure 4) show quite nicely upon comparison with the low temperature measurements ( Figure  9) that the resonance effects are more pronounced for the latter. At higher temperatures, the natural linewidth of the SARs increases, which makes it harder to distinguish them from noise. At lower temperature, mainly elastic scattering channels need to be considered, while at elevated temperatures, the inelastic channels become important also, making the inclusion of these channels in the quantum mechanical calculations necessary.