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 Bi2Se3.[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–Bi2Se3(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 Bi2Se3(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 ϵ, 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 4He. 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 ϑ 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 Bi2Se3 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 103 above the diffuse elastic background (Figure ), 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 Å.

Figure 3

Scattered He intensities for Bi2Se3(111) versus parallel momentum transfer. In the top (bottom) panel, the elastic scan taken along the ΓM (ΓK) azimuth with a sample temperature of 113 K (room temperature) and an incident energy E of 10.07 meV (14.49 meV) is plotted. The green squares show the area of the measured diffraction peaks, while the red stars depict the calculated values. The inset shows the R-factor (eq ) for different values of the corrugation amplitude ξpp using diffraction data at three different incident energies. The minimum of R shifts to higher values of ξpp with increasing E.

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 (Θ = 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 ϑ 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 bywhere k is the incident wavevector, and ϑ and ϑ 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 k2 is positive. If an incident He atom hits the surface, it can undergo a transition into a bound state on the surface with −|ϵ|. The process happens while the He atom is diffracted into a channel, which is kinematically disallowed (k2 < 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 ϵ of the adsorbed atom[23]where m is the He mass. From eq , it is clear that studying SARs provides access to the bound-state energies ϵ(K, G) and, more importantly, to the atom–surface interaction potential.

In this work, we have analyzed SARs upon scattering of He from Bi2Se3 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,48] Likewise, the validity of this potential for TI surfaces[49] and other layered materials such as transition metal dichalcogenides[50] has been proven. In addition, the CMP greatly simplifies the treatment of several steps within the close-coupling (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[51]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 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 V0 (i.e., without corrugation) of eq is given viaThe bound states can be described analytically bywhere n depicts a positive integer, is the Debye frequency, with m being the mass of the impinging He atom, and .

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 close-coupled 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 , the kinematic condition (eq ) can be fulfilled for specific values of ϑ. Two typical diffraction scans for the high-symmetry direction ΓM are shown in Figure . The x-axis has been transformed to parallel momentum transfer using eq . In addition to the diffraction peaks, smaller features caused by SARs can be seen.

Figure 1

Scattered He intensities for Bi2Se3(111) versus parallel momentum transfer ΔK along the ΓM azimuth with a sample temperature of 113 K. The incident beam energy E is 11.25 and 14.26 meV in the top and bottom panels, respectively. 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 . Each line is labeled with the Miller indices of the associated G-vector for the particular resonance condition.

First, we will restrict the kinematic condition given in eq to the free-atom approximation, which assumes that the surface potential is adequately described by the laterally averaged interaction potential (eq ). 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 ϵ(K, G) in eq as constants and is therefore independent of K 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 ) can be rewritten in terms of the incident angle ϑ and the incident wavevector k. 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 ϑ in the scans can be associated with a certain bound state energy ϵ and scattering channel using eq . The position of several SARs is indicated by vertical lines (based on solving eq for ϑ) in Figure and labeled according to their diffraction channel G, while the different colors correspond to different bound states ϵ. 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 is quite ambiguous since numerous combinations of ϵ and G can be thought of to fulfill eq . 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 . To monitor the drift of the resonances versus the incident beam energy, many diffraction scans with various E 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 ϵ, which fulfills eq .

Figure 2

Contour plot of the scattered He intensities (logarithmic scale) versus parallel momentum transfer and incident energy of the He atom along the ΓM azimuth and at a sample temperature of 113 K. The superimposed solid lines correspond to the resonance conditions calculated with the free-atom approximation. Each line is labeled with the Miller indices of the associated G-vector for the resonance condition. The labels ϵ0 – ϵ4 at the top abscissa denote the series of eigenvalues associated with the (10) G-vector. Scattered intensities (mainly the diffraction peaks) exceeding a certain value have been cut off to increase the visibility of the resonance effects. 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 total of 27 ϑ-scans at incident energies ranging from 9.5 to 14.5 meV was collected to construct the contour plot. The x-axis in Figure 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 . 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 ), with the color coding and associated G-vectors labeled in the same way as in Figure . 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 ). The experimentally found bound states are then used to determine the laterally averaged potential (eq ) by minimizing σ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.

Table 1

Experimentally Determined Bound State Values for the Laterally Averaged He–Bi2Se3(111) Interaction Potentiala

bound state	ϵn [meV]	Δϵn [meV]	τn [ps]
ϵ0	5.6	0.39	1.7
ϵ1	3.8	0.33	2.0
ϵ2	2.3	0.30	2.2
ϵ3	1.2	0.27	2.4
ϵ4	0.5	0.13	4.9

The corresponding internal linewidths Δϵ of the bound states (based on the experimental width of the resonances) and their lifetimes τ are also given.

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.[52] 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, and thus, no associated G-vector can be found that would explain such a curvature. One of these features, indicated by the black arrow labeled with KF in Figure , 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 ). Additional features in diffraction scans due to spin-conserving electronic interpocket transitions have been observed for the semimetal Sb(111),[53] and the observability of similar transitions in the TI Sb2Te3 has been suggested by theory.[54] The observation of additional dispersion curves in Bi2Se3[55] 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 Sb2Te3 exhibit narrow electron pockets, while the situation for Bi2Se3 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,56] 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.[55]

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 ). 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 Schrödinger equation together with the CMP. In the exponential term of the potential (eq ), the Fourier series expansion (see eq ) 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,57] (see Details on the Close-Coupling Calculations). The corresponding coupling terms for the CMP can be found in several references.[29,37,48]

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 (Θ = 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 correspondence between measured (Iexp) and calculated diffraction intensities (Isim), is found by minimizing a measure of the deviation Rwith 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 , that is, searching for the global minimum of R in dependence of ξpp. In Figure , two of these diffraction scans along the two high-symmetry directions of Bi2Se3(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 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.[58] Therefore, the above-described optimization routine was repeated by taking into account several experimental spectra recorded around three specific incident energies E, as shown in the inset of Figure . 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 = 10.3 meV), the minimum of R is found with ξpp = 0.21 Å, while for higher incident energies (E = 14.4 meV), it increases to ξpp = 0.27 Å. Optimizing eq 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 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 ) 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. Such a so-called drift spectrum shows again SARs at places where the kinematic condition (eq ) is fulfilled. As can be seen in Figure , 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 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.

Figure 4

Normalized intensity of the specular peak in dependence of the incident wavevector k. Top: Measured drift spectrum along ΓK with a sample temperature of 113 K and corrected for the He flux through the nozzle. The vertical lines show the predicted positions of the SARs using the free-atom approximation. The colors symbolize the different bound-state energies (same nomenclature as used in Figure ), with the dashed lines illustrating the threshold energy. The number next to the lines corresponds to the specific interacting G-vector. Bottom: Simulated drift spectrum using elastic close-coupled calculations with the optimized three-dimensional potential after multiplication with the corresponding Debye–Waller factor. A few vertical lines from the top panel have been extended as a guide to the eye. The turquoise dash-dotted curve illustrates the experimental broadening due to a convolution of the energy distribution in the incident beam with the SARs.

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 TN since the He flow through the nozzle is proportional to . 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 ). 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.[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 #). 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 with Figure .

Figure 9

Determination of the step heights from the oscillations present in the drift scan. The solid lines show the measured specular intensity as a function of the incident wavevector k measured at room temperature with the crystal aligned along the ΓM and ΓK azimuths, respectively. The red dash-dotted line shows intensity oscillations (eq ), which are expected in the presence of terraces separated by steps with a quintuple layer height. The dashed lines show oscillations for quintuple and subquintuple layer steps, with a different number of subquintuple layer steps for ΓM and ΓK.

We turn now to the refinement of the interaction potential based on the drift scan (top panel of Figure ), which shows the measured specular intensity as a function of the incident wavevector k. 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 . 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.[57,59] On the other hand, threshold resonances have only been observed experimentally for scattering He from a ruled grating upon grazing incidence.[60]

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 . In the bottom panel of Figure , 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,[61−63] 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.[52,62] 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.[61] Thus, we will mainly concentrate our analysis on the position of the SARs in terms of k, 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 of the SARs in the simulated data compared to the peak positions in the experimental spectrumwhere kexp are the experimental data values with uncertainties σ (σ ≈ 0.02 Å–1), and ksim are the calculated peak positions. Figure shows the effect on Rp2 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 ) 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 clearly illustrates that the final optimized values of D and κ correspond to a minimum in terms of Rp2.

Figure 5

Effect of varying the potential well depth D and stiffness κ around the optimized potential parameters. The left panel shows a measure of the agreement between the simulations and experimental data, while D is varied using eq on the left ordinate and eq on the right ordinate, respectively. The right panel shows the same reliability factors, now varying κ, while D is held constant.

Figure 8

Simulation of the drift scan using elastic close-coupled calculations for various values of ξpp with D = 6.54 meV and κ = 0.58 Å–1.

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,[64] 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 Iexp″ in the drift scan according towhere the scaling constant c = ∫ Iexpd k/ ∫ Isimd k is used to normalize the experimental and simulated curves with respect to each other.[64,65] 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 = 5.1 and 5.4 Å–1. We see in Figure that the overall trend is the same compared to using Rp as a measure, although the latter seems to be more robust in the present case.

Following this optimization process, the parameters of the final refined three-dimensional atom–surface interaction potential are found as

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 Bi2Te3(111) (6.22 meV)[49] and Bi(111) (7.9 meV),[37] while the stiffness κ of the He–Bi2Se3(111) potential is significantly smaller than for Bi2Te3(111) and for the Bi(111) single crystal, although larger compared to Sb(111).[48]

Because Bi2Se3 has a smaller vdW gap than Bi2Te3 and Sb2Te3,[40,42,66] 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 peak-to-peak corrugationfor an incident beam energy E = 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 Bi2Te3.[49]

Since the spin–orbit coupling in Bi2Se3 is stronger than in Bi2Te3, giving rise to a different electronic structure of the topological surface states,[9,67] 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 Bi2Se3 and Bi2Te3 are quite similar,[45,68] 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 Bi2Se3 and Bi2Te3 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.[57,69] 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 Δϵ of the bound state in terms of energy is then given via[69]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 τ = ℏ/Δϵ.[57]

The natural linewidths and corresponding lifetimes of all bound states are listed in Table . 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,57]

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,48,69,70] Direct experimental information is only available for the He–LiF(001) system[69] and the He–Sb(111) system.[48] The internal linewidths of the He–LiF(001) system are comparable to those found for He–Bi2Se3 in Table . 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 (69) 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 Bi2Se3 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 close-coupled 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 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–Bi2Se3(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 close-coupled 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–Bi2Se3(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 Å.