Generation of spin-polarized hot electrons at topological insulators surfaces by scattering from collective charge excitations

Topological insulators (TIs) are materials which exhibit topologically protected electronic surface states, acting as mass-less Dirac fermions. Beside their fascinating fundamental physics, TIs are also promising candidates for future spintronic devices. In this regard, generation of spin-polarized currents in TIs is the first and most important step towards their application in spin-based devices. Here we demonstrate that when electrons are scattered from the surface of bismuth selenide, a prototype TI, not only the elastic channel but also the inelastic channel is strongly spin dependent. In particular collective charge excitations (plasmons) excited at such surfaces show a large spin-dependent electron scattering. Electrons scattered by these excitations exhibit a high spin asymmetry, as high as 40%. The observed effect opens up new possibilities to generate spin-polarized currents at the surface of TIs or utilize the collective charge excitations to analyze the electrons’ spin. The results are also important to understand the spin polarization of the photo-excited electrons excited at TIs surfaces. Moreover, our finding will inspire new ideas for using these plasmonic excitations in the field of spin-plasmonics. A key feature of topological insulators, such as Bi2Se3 is their protected electronic surface states where the direction of spin is locked to their momentum. Here, using spin-polarized electron energy loss spectroscopy, the authors demonstrate that collective charge excitations can initiate spin-dependent electron scattering at the surface of Bi2Se3, and this occurs for both the elastic and inelastic scattering channels.

G enerally when electrons are scattered from an atom, the scattering cross-section is spin dependent 1 . The same is true when electrons are scattered from an ordered array of atoms, e.g. the surface of a crystal. This leads to different intensities of the scattered spin-up and spin-down electrons [1][2][3][4] . Here the direction of the spin is defined as the projection of the spin operator along the scattering plane's normal vector. In the case of magnetic atoms, the main mechanism leading to this spindependent scattering cross-section is the exchange process 5 . However, in the case of nonmagnetic atoms the underlying mechanism is the so-called spin-orbit scattering. In fact, the difference in the scattering cross-section is a result of the interaction of the electrons' spin with the orbital angular momentum, as it scatters off the atoms. Hence, a large spin-dependent scattering cross-section is expected for the surfaces exhibiting a large spin-orbit coupling (SOC), depending on the symmetry of the crystal and the direction of the scattering plane. The effect has already been used to create a spin-polarized electron beam or to analyze the spin polarization of an electron beam with an unknown spin polarization [6][7][8][9][10] .
Topological insulators (TIs) represent a class of materials with a large SOC. The most prominent feature of TIs is that their electronic surface states are protected by topology [11][12][13][14] . As a consequence, the spin of the surface state electrons is locked to their momentum, making the material robust against nonmagnetic defects because the electrons must undergo a spin reversal when they backscatter [15][16][17][18] . It has been predicted that the collective excitations associated with such electronic states shall exhibit an unconventional spin character, reflecting the chiral spin texture of the surface electronic states [19][20][21][22][23][24] . All these compelling properties of TIs have made them promising candidates for the next generation of spin-based devices, useful for spintronics and topological quantum technologies [25][26][27] . Bi 2 Se 3 is a prototype three-dimensional (3D) TI and exhibits a rather simple band structure and a cone-like Dirac surface state with a Dirac point located below the Fermi level 17 . 3D TIs can host magnetism 28 and superconductivity 29 . This fact makes their integration in spin-based technologies more feasible.
Here we demonstrate that, in the case of TIs, not only the elastically scattered electrons but also the inelastically scattered ones show a large spin-dependent scattering. In particular, electrons scattered by the collective charge excitations exhibit a large spin asymmetry. This implies that in addition to the elastic spinpolarized scattering events one can take advantage of inelastic processes, such as those leading to the creation of collective charge excitations, and thereby produce spin-polarized hot electrons at the surface. The observed effect in combination with the fascinating topological properties of TIs may find applications in spin-based devices. Collective charge excitations can be excited using photons, providing a unique way of photonic-controlled spintronics. Furthermore, the effect is essential to understand the spin polarization of photoexited electrons. Since photoemission experiments are the key to address spin-polarized surface states, it is of prime importance to consider this effect, while deducting the spin polarization of the topological electronic states from the photoemission data.

Results and discussion
Spin-dependent elastic and inelastic scattering. We performed high-resolution spin-polarized electron energy-loss spectroscopy (SPEELS) on Bi 2 Se 3 (0001). The experiments were carried out under ultrahigh vacuum condition on freshly prepared surfaces (see "Methods"). Figure 1a presents the atomic structure of the Bi 2 Se 3 (0001) surface. The corresponding measured and simulated low-energy electron diffraction (LEED) patterns of this surface together with a schematic drawing of the surface Brillouin zone are shown in Fig. 1b. The sharp LEED pattern indicates an atomically ordered and flat surface. It is well known that in situ cleaved samples exhibit a well-ordered surface and are terminated as Se-Bi-Se, which is the expected structure, assuming that the surfaces cleave naturally along the van der Waals gap 30,31 . The spin-polarized elastic and inelastic scattering experiments are performed along the main symmetry Γ-M direction, as depicted in Fig. 1c. Figure 1d shows the intensity of the elastically scattered electron beam as a function of the incident beam energy. The data are recorded for the two opposite orientations of the spin of the incoming electron beam (I ↓ and I ↑ ). The results indicate that the scattering intensity is spin dependent. The asymmetry A ¼ is provided in the lower panel of Fig. 1d. The zero asymmetry is defined as the case in which no difference is observed for the intensity of the scattered electrons when the incoming beam polarization is switched from ↓ to ↑. In order to be able to follow the spin asymmetry over a wide range of incident beam energy, the beam was optimized at two different incident energies (4 and 9 eV) and the results are shown for both cases. The results indicate a spin asymmetry of about 30% for the electrons with an incident energy of about 3-3.5 eV. Likewise for the higher energies between 10.5 and 13.5 eV, one observes also a large spin asymmetry of about 40%. The observed asymmetry of the elastically scattered electrons is due the fact that they are scattered by the relatively large spin-orbit potential of the surface. Although the surface is likely Se terminated, the large SOC, mainly originating from the Bi atoms, acts on the spin of the incoming electrons. This results in a large spin asymmetry, as it is observed in the experiment. Data presented in Fig. 1d can be understood based on a simple scattering model, e.g. the one introduced in refs. [1][2][3][4] . Consider a fully spin-polarized incoming electron beam in the form of plane waves that is scattered from a regular array of atoms. The scattered electrons wavefunctions shall in principle include a spatial part and a spin part, which, in a general case, should include all the possible spin components. In such a condition, the scattering cross-section may be calculated using a classical approach. It is rather straightforward to imagine that the scattering intensity calculated in this way is determined by the coherent superposition of the scattered amplitudes. Consequently, the intensity provides information regarding the scattering geometry, incident angles (the polar θ and azimuthal ϕ scattering angle), incident energy E i , scattering potential as well as the crystal structure. As the scattered wavefunction includes spin components, one can easily imagine that the partial scattering cross-sections will not be identical when the spin-orbit potential is taken into account. This leads to the fact that the scattering cross-sections are different for different spin polarizations of the incoming beam. Hence, a nonzero spin asymmetry is expected. As mentioned above, the intensity is determined by the atomic arrangement and the asymmetry is mainly due to the intrinsic properties of the involved atoms and the strength of the SOC. The variations in the partial intensities observed in the upper panel of Fig. 1d represent mainly the atomic arrangement of the surface atoms and the scattering potential. The spin asymmetry reflects mainly the effect of the SOC of the surface. The spin asymmetry depends on the initial state and hence depends on the energy of the incident electrons as well as the scattering angles (θ and ϕ). In the present case, the conditions are such that the asymmetry exhibits a maximum for incident energies in the range of 10.5 < E i < 12 eV. In principle, a full multiple scattering theory shall reproduce the data presented in Fig. 1d in great details 3 .
As a side remark, a spin-flip scattering process can easily be observed in a 'complete' experiment, i.e. an experiment in which a spin polarized beam and a spin detector is used. In such a case, one would be able to measure all the possible scattering rates, e.g. ↑ to ↑ (non-flip), ↑ to ↓ (flip), ↓ to ↑ (flip) and ↓ to ↓ (non-flip) 5 .
Here ↑ and ↓ refer to the up and down spin states of the incident and scattered electrons, respectively. However, since the scattering cross-section depends on the spin of the incoming beam, the spin-flip scattering process is manifested in the intensities of the scattered beam. In such a case, even if no spin-resolved detection of the scattered beam is performed, one would be able to observe different intensities of the scattered electrons when the spin state of the incoming beam is switched from ↑ to ↓ 32,33 . Hence, one would observe a spin asymmetry associated with the spin-flip scattering processes occurring during the scattering process.
We note that data similar to those reported in Fig. 1d may be obtained by spin-polarized low-energy electron diffraction (SPLEED) experiments. However, in SPEELD one probes only the elastically scattered electrons without any momentum resolution. The experiments are usually performed at higher incident beam energies (typically >10 eV up to 200 eV, where the LEED fine structures are expected to be observed). The main idea here is to unravel the spin character of collective surface excitations. This is only possible if one can probe both the elastically and inelastically scattered electrons with a very high spin, energy and momentum resolution.
In order to understand the effect of SOC on the inelastic excitations, we investigated the spin resolved spectra recorded at various incident energies and for different momentum transfers of the incident electron beam. Figure 1e shows the SPEEL spectra recorded at an incident electron energy of 3.5 eV and at the specular geometry. In each spectra, one observes the elastic peak located at the zero energy loss. In addition, one observes also an intense peak at an electron energy loss of about 90 meV. This peak is associated with the collective charge excitations at this surface (see the discussion below). The most important result is that electrons scattered by these collective excitations show also a high spin asymmetry as shown in the lower panel of Fig. 1e. Note that the fine details of the asymmetry spectra are the consequence of several effects. For a detailed description of the features observed in the asymmetry spectra, we refer the reader to Supplementary Note 1.
Origin of the observed loss features. In order to unambiguously identify the origin of the feature at the energy loss of about 90 meV, additional experiments were performed with a higher energy resolution at different energies and momentum transfers. Moreover, the loss spectra were numerically calculated using the semi-classical dielectric response theory, as described in "Methods". In the calculations, one Drude term was considered to describe the electronic contribution and two Lorentz terms to describe the ionic contribution to the dielectric function. It turned out that the most important term to consider is the Drude term caused by the electronic contributions. The other terms have only a minor effect on the final loss spectra. This assumption is fully supported by the optical measurements showing the large spectral weight of the electronic contribution to the optical conductivity data [34][35][36] . The most important input parameter for our calculations is the carrier density n, which has been independently measured by our Hall effect measurements. All the other parameters are taken from the available experimental data 17,34,35,37,38 . A detailed list of the parameters used for our calculations is provided in Supplementary Table 1).
The results of both experiments and calculations are summarized in Fig. 2. Figure 2a presents the experimental spectrum measured at the Γ-point and at an incident energy of E i = 4.0 eV. The numerically calculated spectrum for this Fig. 1 Surface structure, spin-dependent elastic and inelastic scattering on Bi 2 Se 3 (0001). a The atomic structure of the Bi 2 Se 3 (0001) surface and b its corresponding low-energy electron diffraction (LEED) pattern, recorded at an electron energy of 78 eV. The experimental (left) and simulated (middle) LEED patterns are shown together with a schematic drawing of the surface Brillouin zone (right). c The scattering geometry. θ denotes the scattering angle. The red and blue arrows represent the spin of the incoming beam parallel and antiparallel to the scattering plane's normal vector n. The scattering plane is shown by the shaded light-blue area. E i (E f ) and k i (k f ) denote the energy and the momentum of the incident (scattered) electron beam, respectively. d Intensity of the elastically scattered electron beam as a function of the incident beam energy. Data shown by filled and open symbols were recorded when the beam was optimized at 4 and 9 eV, respectively. e Spin-polarized electron energy-loss spectra recorded at a beam energy of 3.5 eV and at a momentum transfer of q = 0.05 Å −1 . In d, e, the incident angle, i.e. the angle between the incident beam and the surface normal, was set to θ = 40°. The scattering plane was parallel to the high symmetry direction Γ-M of the surface Brillouin zone. ↓ and ↑ represent the spectra when the spin polarization of the incoming electron beam was parallel (↓) and antiparallel (↑) to the scattering plane's normal vector n. The upper panels show the spin polarized and the total (I ↓ + I ↑ ) spectra. The lower panels show the asymmetry A ¼ spectra. The error bars represent the statistical uncertainties. The strong noise (or large error bars) in the asymmetry spectra for the loss energies >200 meV is due to the fact that the total intensity decreases substantially. The total intensity appears in the denominator of spin asymmetry. Hence, the larger uncertainties in the total intensity lead to larger errors in the asymmetry. geometry is presented in Fig. 2b. Our calculations reproduce the measured loss spectra rather well. The region showing the plasmonic excitations is magnified in Fig. 2c, d. In addition to the main loss feature at about 90 meV, our calculations predict the higher-order plasmonic excitations at about 180 meV ( Fig. 2d). At this energy, one also observes a small feature in the experimental spectra (Fig. 2c).
It is important to mention that the calculations of the energyloss spectra based on the approach described above are only valid for the dipolar scattering regime, i.e. close to the elastic scattering (the Γ-point). In order to describe the measurements at the M-point, we construct the loss spectra considering the following assumptions.
(i) We consider the loss feature, which satisfies the condition Re εðω; qÞ Â Ã ¼ À1. This part describes practically the loss peak observed also at the Γ-point.
(ii) Since in the present case the momentum transfer is rather large, one expects that the excited plasmonic excitations propagate along the interface between the depletion layer and the bulk. The condition under which a loss feature shall be observed is given by Re εðω; qÞ Â Ã ¼ Àεð1Þ. Such a condition would lead to a plasmonic loss to be observed at ω sp j q≠0 ¼ ω p = ffiffiffiffiffiffiffiffiffiffiffiffi ffi 2εð1Þ p 39 . (iii) Since the contribution of the quasielastic scattering at the M-point is substantially reduced, one would also observe the optical phonons of the system 40 . This fact can be considered by adding a phonon line at the given phonon frequency.
We, therefore, construct the loss spectrum as an overall spectrum of the above-mentioned excitations. The results are summarized in Fig. 2f. Figure 2e shows the experimental loss spectrum measured at the M-point. One observes three main loss features at 22, 71 and 96 meV. The first peak is due to the excitation of an optical phonon as has also been observed in other experiments 37,[41][42][43][44][45] . The last two peaks are the results of the plasmonic excitations satisfying the two conditions mentioned above. The calculated loss spectrum is shown in Fig. 2f. The black curve represents the case in which the elastic contribution is not subtracted. The blue curve shows the case in which the elastic contribution is subtracted. We note that typically the momentum resolution of the electron scattering experiments, like ours, is on the order of 0.03 Å −1 . Hence even at the specular geometry one shall observe contributions from both plasmonic modes. This may also explain the slight lower energy and the larger broadening of the peak measured at the Γ-point [Fig. 2a, c].
The final conclusion of the comparison provided above is that the loss feature observed in the experimental spectra is predominately determined by the collective excitations of all charge carriers, i.e. the surface plasmonic excitations of the bulk electronic states projected onto the surface as well as the Dirac states. The signature of plasmonic excitations has also been observed in earlier non-spinpolarized high-resolution electron energy-loss spectra 22,43,44 and also in optical experiments 21,23,24,34,35,38,46 . It has been discussed that the observed excitation may even be a mixed state of the socalled Dirac plasmons of the topologically protected Dirac electronic states as well as the ordinary surface plasmons of the 3D charge carriers of the projected bulk states 22 . We do not aim to disentangle these two. The most important feature that we aim to discuss here is The black curve in f is the calculation including the elastic peak. In the blue curve, the elastic contribution is removed. Note that in the spectra presented in a, c, e the error bars are also shown. Since they are smaller than the symbol size, they are not visible. The error bars represent the statistical uncertainties of the measured count rates.
the spin dependence of these excitations at TI surfaces and the large observed spin asymmetry, while they are excited by spin-polarized electrons.
Another important result of this comparison is that, since the energy (frequency) of these excitations is directly linked to the carrier concentration n, this provides a unique possibility and easy way to vary the plasmon frequency over a wide range of frequency. This is important from the application point of view. In principle, by just changing the carrier concentration one can tune the excitation frequencies going towards the frequencies where the optoelectronic (or opto-spintronic) devices are usually operated.
Origin of the observed high-spin asymmetry. Irrespective of the fine details of the plasmonic modes observed in the experiment and described above, the most prominent feature of these excitations is that they exhibit a rather large spin asymmetry (see, for example, the lower panel of Fig. 1d). Note that spectra recorded at higher incident energy (10.5-13 eV) show also a large spin asymmetry in the region of the plasmonic peaks. The spectra look very much the same as those shown in Fig. 1e. The value of the spin asymmetry at the peak associated with the plasmonic excitations is very similar to that of the elastic peak at zero energy loss. In order to investigate the contribution of the collective charge excitations in the observed spin asymmetry, one may tune the energy of the incoming electron beam such that the zero-loss peak shows no spin asymmetry. Looking at the data presented in Fig. 1c, this is possible by tuning the energy of the incident beam to E i = 4.0 eV, where the observed spin asymmetry of the elastic scattering is zero. In Fig. 3, the map of total intensity and spin asymmetry in the vicinity of the Γ-point recorded at E i = 4.0 eV is presented. The spectra were recorded in the momentum space by changing the scattering geometry, i.e. a continuous rotation of the sample about its main axis in steps of 0.5°. In each step, the spectra were recorded for the two different spin directions of the incoming beam. For each individual energy-loss value, two intensities were recorded, subsequently. In the total intensity map presented in Fig. 3a, one observes the (quasi-)elastic peak, located at the energy-loss of zero. In addition to that, one observes a large intensity at an energy loss of about 90 meV as a result of collective charge excitations. Looking at the asymmetry map presented in Fig. 3b, one realizes a high-spin asymmetry in the energy range, where the collective charge excitations are located. At the zeroloss peak, one observes a very low-spin asymmetry. Details of the asymmetry for loss energies <20 meV are mainly governed by the presence of the phonon peaks of the system. At 30-35 meV, the onset of plasmonic excitations appears and consequently the asymmetry increases. At high loss energies, losses due to multiple plasmonic excitations as well as single particle expiations are possible. These excitations lead to the remaining large spin asymmetry for higher energy losses. As the total intensity decreases with the energy loss, this leads to a slight linear increase of the asymmetry (for a detailed description of the asymmetry spectra, see Supplementary Note 1). The large asymmetry of the inelastically scattered electrons by collective charge excitations indicates that the electrons travelling in front of the surface can undergo a spin reversal, while scattered by such excitations. Such a process can lead to a high asymmetry at the plasmonic losses. We would like to emphasize that the spin-flip scattering observed here is of spin-orbit nature, meaning that the electrons that contribute to the plasmonic excitations are then scattered by the spin-orbit potential of the surface. It is important to mention that, if these electrons would have been scattered elastically or quasi-elastically, they would not have shown a large spin asymmetry. Electrons that contribute to the excitation of the collective charge excitations and travel in front of the surface have the possibility to undergo a spin-flip spin-orbit scattering and hence produce a considerably large spin asymmetry.
In our model used to describe the plasmonic peaks, the spin dependence of the electronic surface states is not taken into account (calculations presented in Fig. 2). In fact, a full theoretical consideration of spin-dependent inelastic electron scattering together with an accurate description of the spin-dependent electronic surface states would explain the observed phenomena in great details. However, such an advance calculation requires: (i) a detailed consideration of multiple scattering processes with an unambiguous consideration of the involved scattering matrix elements (not only elastic but also inelastic) and (ii) a detailed consideration of all electronic states at the surface (including their spin character). Note that the surface charge carriers do not only include the Dirac states but also the projected bulk states at the surface. Although all states are spin polarized (since they are subject to a large SOC), in practice, it is not easy to distinguish between the Dirac states and the other electronic states. The evidence for the importance of the consideration of all the electronic surface states comes from our calculations of the loss spectra. The input of our calculations is the carrier concentration measured by the Hall effect measurements (without considering their spin character). The agreement between the calculated and measured spectra indicates that indeed all these states are important and have to be considered in order to adequately explain these plasmonic features.
Experiments performed at higher wavevectors reveal that the intensity of the peak associated with the collective charge excitations reduces with wavevector. This is expected from the semi-classical dipolar-scattering theory. The dipolar-scattering contribution strongly reduces while going from the specular geometry towards higher momentum transfers. The most interesting observation is that, meanwhile, the spin asymmetry increases. For example, at q = 0.07 Å −1 the asymmetry is as high as 30%. This large asymmetry is observed in the loss region where the collective charge excitations exist. In order to illustrate this, the asymmetry spectra recorded at 0, 0.04, 0.07 and 0.1 Å −1 are shown in Fig. 3c (for individual spectra, see Supplementary Note 1 and Supplementary Fig. 1). We attribute this effect to the fact that at higher momentum transfers the electrons contributing to the collective charge excitations exhibit a longer travelling distance (time) in front of the surface and hence the probability of contributing to a spin-orbit spin-flip process becomes higher. The reason is that when electrons transfer a certain momentum to the system the parallel component of the scattering momentum becomes smaller. The maximum asymmetry is observed at slightly lower energy losses than the main peak. This is due to the fact that the main loss peak is, in fact, a superposition of two modes, with different intensities (see Supplementary Note 1 and Supplementary Fig. 1).
The description of the spin asymmetry of the plasmonic peak may seem to be similar to that of the elastic scattering (both are of spin-orbit nature). However, they are not exactly the same. This is evident by the fact that the observed spin asymmetry of the inelastic scattering does not show the same E i dependence as the one of the elastic scattering. The plasmonic excitations exhibits a large spin asymmetry even when the incident beam energy is tuned such that the elastic spin asymmetry is negligible. This observation may be understood by the two following scenarios. (i) The dynamical electric field generated by the response of the collective charge excitations to the incident electrons has a spin-orbit-induced component, which can differently act on the electrons having opposite spins. This means that the presence of such a dynamical field alters the scattering cross-section. The cross-section is differently altered for the two spin orientations of the incoming electron beam. (ii) From a classical point of view, the energy-loss process associated with the excitations of the plasmonic excitations slows down the electron so that it can better be "captured" by the spin-orbit potential of the surface. We note that the picture provided above is oversimplified and a more in-depth theoretical confirmation is indeed very helpful here.

Conclusion
In conclusion, we showed that electrons scattered by collective charge excitations at the Bi 2 Se 3 (0001) surface exhibit a large spin asymmetry. We anticipate that the observed phenomenon is general and can occur at the surface of other TIs and twodimensional (2D) quantum materials with a large SOC. The importance of the effect is threefold. First, it can be used to generate spin-polarized hot electrons at TI surfaces. Consider an unpolarized electron beam that is scattered by the collective charge excitation at the surface of TIs. Since electrons with opposite spin exhibit different scattering cross-sections, one can produce a spin-polarized beam in this manner. In a similar manner, one may use the effect to analyse the spin of an unknown electron beam. Generation of spin-polarized hot electrons together with the fascinating topological properties of TIs may find its applications in spin-based devices. As such spin-dependent plasmons may be excited using photons, this would open up new opportunities for photonic-controlled spintronics and lead to the creation of an exciting research direction. One of the important properties of these plasmonic excitations is that their frequency is directly related to the carrier density. As demonstrated by our theoretical calculations of the loss spectra, this quantity can easily be tuned by tuning the doping. This provides a unique possibility and easy way to vary the plasmon frequency over a wide range of frequencies going towards the frequencies where the optoelectronic (or opto-spintronic) devices are operated. Moreover, the position of the Dirac point of the material can also be tuned by tuning the material's composition, as has been shown in ref. 47 . This fact indicates that the observed effect is not limited to the studied compound and can be extended to a large family of TIs. The effect can be tailored to a vast variety of 2D and quantum materials. Moreover, the advance fabrication techniques provide a great opportunity to tailor the effect to low-dimensional topological solids in the form of thin films, multilayers and nanostructures. With this, one would be able to also couple different plasmonic modes originating from different interfaces. Second, the effect must be considered while interpreting the spin dependence of photoexcited electrons 17,48-50 . Generally when an electron is photoexcited from a state below the Fermi level to a state above the vacuum level, it can, in principle, couple to these plasmonic excitations and thereby lose energy. Meanwhile, the electron's spin state may be changed. The process is highly probable due to the fact that the coupling of electrons to such excitations is rather strong. This strong coupling is evident from the intense peaks observed in the loss spectra [see, for example, Fig. 3a]. The effect can, therefore, lead to a fictitious spin polarization, which has nothing to do with the initial spin state. All these need to be considered for an unambiguous interpretation of the spin-dependent processes in all the experiments based on the photoexcited electrons, e.g. spin-resolved photoemission or dichroism experiments. Third, the observed effect may inspire ideas in the emerging field of spin plasmonics. Since these plasmonic excitations enable optical control of spin states, they can open new perspectives for photonic-controlled spintronics. As the excitation cross-section of these plasmonic excitations is strongly spin dependent, one would be able to couple the electrons' spin to these plasmonic modes. Note that the intensity of the plasmonic excitation is 5% of the elastic peak. This is a huge value and clearly indicates that these excitations can be excited very efficiently by electrons. The experimental spectra show a spin asymmetry of about 40%, which is also a very high value. One may either utilize the spin-polarized electrons to excite a certain plasmonic mode or, vice versa, generate a spin current via an excited plasmonic mode. The plasmon frequency can be tuned by tuning the carrier concentration. Combining the observed effect with the topological character of TIs would lead to new functionalities. This would even be more attractive by combining TIs and magnetic materials 51 .

Experiments
Sample preparation. All the spectroscopic experiments were performed under ultrahigh vacuum (UHV) conditions. Single crystalline Bi 2 Se 3 samples were loaded into the UHV chamber and were cleaved under UHV conditions at room temperature. This cleaving method results in an atomically clean Se-terminated Bi 2 Se 3 (0001) surface with relatively large traces of one quintuple layer height 30,31 . The high quality of the surface is evident by our LEED patterns showing very sharp (1 × 1) spots. The cleanliness of the surface was checked by Auger electron spectroscopy indicating a clean and contamination-free surface. Experiments were carried out on a freshly prepared surface, immediately after cleaving.
Spin-polarized electron scattering experiments. A spin-polarized monochromatic electron beam with an energy resolution between 5 and 11 meV was used 32,33 . The spin-polarized electron beam is generated by photoemission from a strained GaAsP heterostructure 52 . The spin polarization of the incoming electron beam is either parallel or antiparallel to the scattering's plane normal vector. The degree of polarization was estimated by performing spin-polarized elastic scattering from a W(110) single crystal, resulting in a value of about 70%. The total scattering angle, i.e. the angle between the incident and the scattered beam, was set to 80°. The scattering intensity was recorded simultaneously for two different spin polarizations of the incoming electron beam. More precisely, first the scattering geometry was fixed and then the intensity of the scattered electrons was measured as a function of the energy-loss. When recording the intensity of the scattered electron beam for each value of energy loss, two values for the intensity were measured: one for spin-up electrons and the other one for spin-down electrons. Flipping the beam polarization was realized by reversing the helicity of the laser beam used to emit the spin-polarized electrons from the photocathode. This means that the two spectra (one for spin-up and the other for spin-down electrons) were recorded at the same time within a single run of the measurement. In order to ensure a reasonably high intensity, the incoming beam was optimized at two different energies (4 and 9 eV). The scattered beam was energy analysed and its intensity was measured without further spin analysis.
In order to measure the inelastic excitations, the spectra were recorded in offspecular geometry at a certain wavevector transfer q ¼ k i sin θ À k f sinðθ 0 À θÞ, where k i (k f ) is the magnitude of the wavevector of the incident (scattered) electrons, and θ (θ 0 ) is the angle between the incident beam and sample normal (the scattered beam). Different wavevectors were achieved by changing the scattering angles. In the experiments, the total scattering angle θ 0 was kept at 80°a nd θ was changed by rotating the sample. The spectra were recorded along the main symmetry Γ-M direction of the surface Brillouin zone. The momentum resolution of the experiment is given by Δq ¼ ffiffiffiffiffiffiffiffiffiffi 2mE i p =_ðcos θ þ cosðθ 0 À θÞÞΔθ. E i represents the energy of the incident beam and Δθ depends on the spectrometers (Δθ = 2°in our case). The momentum resolution of the spectrometer is about 0.03 Å −1 .
Hall effect measurements. The 3D carrier density was measured by probing the Hall coefficient of the sample. The standard Hall experiments were performed on the same sample with the lateral dimensions of 5 × 3 mm 2 and the thickness of 0.5 mm. We measure a Hall coefficient of R H = −18.24 × 10 −8 Ωm/T at 20 K and R H = −16.7 × 10 −8 Ωm/T at 300 K, determined from the linear behaviour of the Hall resistivity as a function of the magnetic field. The carrier density n is inversely proportional to the Hall coefficient. Our analysis results in carrier concentrations of n = 3.42 × 10 25 m −3 and n = 3.7 × 10 25 m −3 at T = 20 K and T = 300 K, respectively.
Theory. In order to understand the origin of the peaks observed in the energy-loss spectra, we have developed a numerical scheme to calculate the loss spectra. The scheme is based on the semi-classical dielectric response theory. Electrons scattered from a surface interact with the dynamical electric field at the surface leading to the excitation of the collective charge excitations. Note that the so-called dipole mechanism is the main mechanism for the specular beam and is mediated by the Coulomb interaction. As the Coulomb interaction is rather long range and the time scales involved are rather long (about 10 fs), the microscopic details of the interaction potential are not needed in order to describe this scattering mechanism. Using the fluctuation-dissipation theorem, the correlation function for the electric field fluctuations is related to the dielectric properties of the medium. Therefore, the classical single-loss probability of a scattered electron at 0 K and close to the specular geometry (q → 0) can be calculated by where a 0 is the Bohr radius and ε(q, ω) is the dielectric function of the medium and in its general form is given by Here the first term represents the contribution of the charge carriers to the dielectric function. ε(∞) is the high-frequency dielectric constant and ω p is the plasma frequency of the free carriers and is related to the carrier density n, the effective mass m* and the vacuum permittivity ε 0 by ω p ¼ ffiffiffiffiffiffiffi ne 2 ε 0 m Ã q . The quantity γ p is the linewidth of the plasmon peak and is determined by the plasmon relaxation time. The second term represents the Lorentz oscillators, describing the ionic contribution to the dielectric function via m transversal optical phonons. Q j , ω TO,j and γ TO,j represent the oscillator strength, the phonon frequency and linewidth, respectively. Equation (1) describes only the so-called single-loss probability for an electron with the wavevector k i to be scattered inelastically from a semi-infinite surface and lose the energy ω at T = 0 K. The multiple scatterings, the zero-loss peak, and temperature effects are included using the formalism developed by Lucas and Šunjić [53][54][55] .
The most important input parameter for our calculations is n, which has been independently measured by our Hall effect measurements. m* = 0.19m e , ε(∞) = 29 and the phonon frequencies and linewidths are taken from the available experimental data 17,37 . We consider one Drude and two Lorentz terms in the description of the dielectric function. The parameters used for our numerical calculation of the loss spectra are listed in Supplementary Table 1).

Data availability
The data sets generated and/or analysed during the current study are available from the corresponding author on reasonable request.

Code availability
The codes associated with this manuscript are available from the corresponding author on reasonable request.