Observation of a topological insulator Dirac cone reshaped by non-magnetic impurity resonance

The massless Dirac electrons found at topological insulator surfaces are thought to be influenced very little by weak, non-magnetic disorder. However, a resonance effect of strongly perturbing non-magnetic impurities has been theoretically predicted to change the dispersion and physical nature of low-energy quasiparticles, resulting in unique particle-like states that lack microscopic translational symmetry. Here we report the direct observation of impurities reshaping the surface Dirac cone of the model three-dimensional topological insulator bismuth selenide. A pronounced kink-like dispersion feature is observed in disorder-enriched samples, and found to be closely associated with the anomaly caused by impurity resonance in the surface state density of states, as observed by dichroic angle-resolved photoemission spectroscopy. The experimental observation of these features, which closely resemble theoretical predictions, has significant implications for the properties of topological Dirac cones in applied scenarios that commonly feature point-defect disorder at surfaces or interfaces. The electronic properties of topological insulators are robust against perturbations, including the presence of non-magnetic impurities. However, surface impurities can give rise to resonant states near the Dirac point, and if their density becomes sufficiently high it is predicted that they can substantially modify the dispersion of the Dirac cone and develop a collective behaviour that results in the formation of particle-like states that lack microscopic translational symmetry. L. Andrew Wray at Purdue University and at the New York University Shanghai, and colleagues, used angle-resolved photoemission spectroscopy to experimentally observe the reshaping of the surface Dirac cone in a defect-rich sample of the topological insulator Bi2Se3. These results indicate that surface impurities can provide a useful handle to control the properties of topological insulators.


INTRODUCTION
Three-dimensional topological insulators are bulk semiconductors with spin-helical Dirac cone surface states that span the bulk band gap. 1,2 Since their discovery around 2007, 3-6 the topological Dirac cone has appeared at the heart of a wide range of proposals for novel emergent quasiparticles and next-generation electronics, such as for the realization of exotic dyon-based, axion-based, or Majorana fermion-based physics. 1,2,7-10 Moreover, the spin-helical Dirac surface states are very robust against dilute non-magnetic impurities due to intrinsic immunities to backscattering and Anderson localization, [11][12][13][14][15][16] and numerous studies have shown the surface Dirac cone remains qualitatively intact in the presence of weak non-magnetic disorder. [17][18][19] However, local probe studies have found that the effect of disorder on the real space electronic structure can be remarkably strong. Scanning tunneling microscopy (STM) experiments and complementary theoretical works have identified that even the simplest form of the non-magnetic impurity, a crystallographic point defect, will tend to give rise to resonance states very close to the Dirac point in two-dimensional (2D) Dirac fermion systems such as topological insulators and graphene. [20][21][22][23][24][25][26][27] In systems like bismuth selenide (Bi 2 Se 3 ), these surface defects are associated with distinctive patterns in STM topography images, and can redistribute local density of states (DOS) by building up a new peak tens of millielectronvolts above the surface Dirac point.
The relative significance of impurity resonances for real space electronic structure as opposed to momentum space is easy to understand if the impurities are low in density, because the associated resonance states will adhere locally to the sparse defects. However, if there is a sufficiently high density of nonmagnetic impurities, the dispersion of the topological surface Dirac cone is expected to change significantly. Recent models have shown that when the density of impurities reaches an experimentally achievable high level, the impurity resonance states will behave collectively much like a flat band that hybridizes coherently with the upper Dirac cone, breaking the Dirac band into upper and lower branches. 28,29 Measurements and theory suggest that electrons near the resonance adhere spatially to the disordered impurity lattice, which lacks translational symmetry. This scenario is at odds with the standard definition of a quasiparticle as a near-eigenstate of momentum. The momentum operator (K) is the generator of translations, and the spatial translation operator is defined as T x ð Þ ¼ expðÀix Á KÞ, implying that (near-) eigenstates of one operator should be (near-) eigenstates of the other. When an electron system is so disordered that translational symmetry in single-particle wavefunctions is absent on the length scale of one de Broglie wavelength, the result is termed a "bad metal," and it is assumed that the quasiparticle picture no longer applies. However, the resistance of topological surface electrons to Anderson localization and backscattering is thought to result in a unique quasiparticle-like character for electrons at the impurity resonance. Though these electrons profoundly lack translational symmetry, simulations suggest that their width cross-section in momentum space is narrower than one inverse wavelength, 28,29 and meets the nominal definition for a "good" quasiparticle. This is quite surprising, and means that inducing disorder at a topological insulator surface may enable the first experimental realization of an itinerant quasiparticle that propagates with a well-defined momentum, but occupies a spatial basis that lacks even nearneighbor translational symmetry. A quantitative correspondence with STM measurements of the spatially resolved DOS distribution has been cited as experimental evidence for this novel quasiparticle character. 28 However, the larger scale morphology of such highly disordered samples has been problematic for highresolution angle-resolved photoemission (ARPES), and a topological Dirac cone reshaped by non-magnetic disorder has not been previously reported.
In this study, we present a high-resolution linear-dichroic ARPES (LD-ARPES) investigation of defect-enriched Bi 2 Se 3 , and report the experimental observation of a topological surface Dirac cone reshaped by non-magnetic impurity resonance. The LD-ARPES spectra reveal an anomalous kink-like feature in the Dirac cone dispersion 40 meV above the Dirac point, which matches the predicted signature of coherent hybridization with an impurity resonance. Integrating over momentum shows that the dispersion anomaly is associated with a DOS peak, and that both features are progressively attenuated when successive intervals of low temperature (LT) annealing are applied to reduce disorder. All of these experimental results are highly consistent with theoretical predictions and with previous STM findings.
Modeling impurity resonances in a topological surface Dirac cone Bi 2 Se 3 is widely seen as a model topological insulator system, with a single Dirac cone surface state connecting across a bulk band gap of roughly 0.3 eV. An idealized model of the ARPES spectral function of the linearly dispersive Bi 2 Se 3 surface state is shown in Fig. 1a. When a surface is simulated with randomly distributed non-magnetic impurities represented by scalar delta-function potentials, the Dirac bands are broadened and their dispersion is changed (Fig. 1b). The bands develop an apparent kink at E4 0 meV above the Dirac point, corresponding with a blurred region in the simulated spectrum (dashed lines in Fig. 1b). Weighting the spectra with the participation ratio (PR) of each eigenstate, a method to reveal spatially inhomogeneous states, 30 highlights the impurity resonant states as they are more concentrated around the impurity sites (Fig. 1c). In the weighted spectrum, the blurred feature then can be distinguished as being composed of two bands (white dashed lines) that are broken by the impurity resonance. The two disconnected dispersions are easier to distinguish at higher defect densities, 28 and the corresponding features broaden as they approach the resonance energy of E~40 meV. Integrating the simulation over the 2D momentum space extracts the DOS, which shows a hump at the resonance energy (Fig. 1d) where a local DOS peak has been observed for point defects by STM. [20][21][22]

RESULTS
Crystallographic point defects occur in typical bulk-grown Bi 2 Se 3 samples with relatively low densities of~<0.05% per 5-atom formula unit. In this study, defect-enriched Bi 2 Se 3 bulk samples are synthesized by abbreviating the final annealing stage of sample growth (see Methods), following the procedure in ref. 31 . X-ray diffraction (XRD) data (Fig. 2a) and STM topography (Fig. 2b) indicate that the resulting Bi 2 Se 3 samples realize a large density of point defects while maintaining good single-phase crystallinity. The resulting defect species are labeled on the STM topography map, based on previously identified correspondences. 32 The primary type of impurity is interstitial Se atoms residing on the outer surface or between the first and second Bi 2 Se 3 quintuple layers, with a combined density of ρ~0.08% per surface unit cell (i.e., per 0.15 nm 2 surface area). A small portion of anti-site Bi Se and Se vacancy defects are also observed with densities of ρ~0.03% and ρ~0.01%, respectively. This defect density has a good correspondence with the bulk Hall carrier density, but the precise numbers seen by STM fluctuate considerably from region to region (see Methods). The theoretical precondition for impurities significantly reconstructing the electronic structure is that there must be a defect population with local resonances at approximately the same energy E R relative to the Dirac point, and a density that exceeds a cutoff proportional to (E R ) 2 . 28 In the sparsedefect limit, E R is proportional to the negative inverse of the effective interaction strength between defects and the surface state (U eff , defined in Methods), 27 suggesting that defects can be neglected if the ratio ρ/(U eff ) 2 is significantly smaller than a critical threshold. Moreover, the surface state skin depth in Bi 2 Se 3 and ARPES measurement depth are both limited to <~1 nm, 33 meaning that the DOS contribution from deeper-lying resonances will not be strongly observed. Based on these considerations, the 2D density adopted for simulations is ρ = 0.06%, representing the assumption that roughly half of the impurities observed in the top nanometer of the crystal will effectively share a common resonance energy, and collectively exceed the critical threshold.
High-resolution LD-ARPES was used to perform a targeted study of defect-derived changes in the surface electronic structure. Measurements were performed at the Advanced Light Source MERLIN beamline (BL4.0.3), which also provides a small 50 μm beam profile to minimize feature broadening from macroscopic inhomogeneity. Two linearly polarized photoemission experimental geometries were selected as shown in the Fig. 2c, d. The σpolarization condition places the electric field parallel to the Fig. 1 Impurity resonance in momentum space. a A simulated Bi 2 Se 3 surface Dirac cone without impurities. Red lines trace the linear band dispersion. b A simulated Dirac cone reshaped by scalar impurities (U eff = 1 eV, density ρ = 0.06% per 2D surface unit cell). A kink-like dispersive feature occurs at the impurity resonance energy, and is traced with dashed lines. c The simulated Dirac cone is weighted by participation ratio to highlight impurity-resonant states. The kink-like feature is resolved more clearly as containing the split dispersions predicted in ref. 28 (white dashed lines). d The corresponding simulated DOS curves of (blue) a pristine Dirac cone, (red) a Dirac cone with impurities, and (yellow) a participation ratio (PR) weighted simulation with impurities Observation of a topological insulator Dirac cone L Miao et al.
sample surface, and the π measurement condition polarization is inside the scattering plane, with a primarily out-of-plane (crystalline c-axis) orientation. Incident photon polarization was switched between π and σ geometries via control of the synchrotron beamline elliptically polarized undulator, to keep exactly the same beam spot on the sample. Previous ARPES measurements on the topological Dirac cone mainly used the π-pol experimental geometry, which gives strong emission from the primary P zorbital component of the Dirac cone. The σ-pol condition is rarely used due to its low efficiency. However, the weakness of emission from the electronic structure of a pristine Dirac cone is also a positive feature of this polarization condition, as the derivative term in the electron-photon interaction couples in-plane polarization directly to nanoscale in-plane inhomogeneity, which is greatest within defect resonance states and defect-rich regions of the sample surface. The asymmetric properties of these polarization matrix elements are reviewed in the Online Supplementary Material, and qualitatively associate the π and σ geometries with selective sensitivity to defect-poor and defect-rich surface regions, respectively.
High-resolution ARPES measurements of the Dirac cone at a freshly cleaved sample surface are shown for the π-pol and σ-pol geometries in Fig. 3a. The surface Dirac bands cross at the Dirac point and the ARPES matrix element effect gives the bands different intensity profiles under different polarizations. The lower Dirac cone is more difficult to trace due to the close proximity of the bulk valence band, and both the lower Dirac cone and the valence band show much higher intensity in σ-pol measurements. Momentum distribution curves (MDCs) of the Dirac bands are shown in Fig. 3e, and show a large discrepancy between σ and π polarizations. At energies greater than E~> 30 meV relative to the Dirac point, the σ-pol ARPES-imaged Dirac bands always appear to be centered at smaller momenta, indicating a dispersion different from the states seen under π-pol. The MDCs of the simulated impurity-rich surface show a spectrum very similar to the Dirac cone imaged under σ-pol. In the simulation (Fig. 3d), the Dirac bands have smaller momenta above the resonance energy (E~> 30 meV) within the kinked Dirac cone. This polarization dependence between the dispersion of ARPES-imaged Dirac bands is not possible for a perfectly homogenous sample, but the σ-pol dispersion closely resembles the emergent "kink-like" feature associated with higher impurity densities. To reveal the quasiparticle dispersions more clearly, the MDCs of simulated and dichroic ARPES-imaged Dirac cones (Fig. 3d, e) are fitted with Voigt functions to track their dispersions (Fig. 3h, i). For simplicity, the "kink-like" feature in the Dirac cone with resonance states was also treated as being composed of just two peaks for the purposes of the fitting procedure (see Supplementary Fig. 1 for the curve- common. An expanded inset shows distinctive defect profiles used for characterization. c The π-polarization ARPES measurement geometry, with the electric field of incident photons mostly normal to the sample surface (projecting 82% onto the z-axis). d The σ-polarization ARPES measurement geometry gives an electric field that projects 100% onto the in-plane y-axis (the ARPES analyzer slit axis) Observation of a topological insulator Dirac cone L Miao et al.
by-curve fitting of MDCs). The σ-pol and π-pol dispersions have a relatively constant momentum offset at energies high above the Dirac point. At lower energies approaching the Dirac point, the σpol group velocity appears to become very large, as is common on the low-energy side of dispersion kinks, and the two dispersions merge rapidly beneath E~40-50 meV. These anomalous features of the σ-pol dispersion can be closely reproduced by the simulation of a Dirac cone with dense impurities (Fig. 3h).
LT annealing provides a relatively safe way to mobilize the interstitial Se atoms in the van de Waals layer without creating new impurities of Se vacancy or Bi Se anti-sites. The temperature T 120°C is chosen to remain well below the thermal activation energy for creating new point defects, and beneath the T~> 150°C low-energy cutoff for eliminating larger morphological defects. 34,35 To reduce disorder, the same sample was treated with intervals of LT annealing between synchrotron beamtimes, and re-cleaved for each new ARPES measurement. In a second experiment following 1 h of LT annealing, ARPES data (Fig. 3b) show significantly sharper bands, which is consistent with a reduced defect density. From the dichroic MDC comparison (Fig.  3f), the discrepancy of dispersion between the π-pol and σ-pol imaged Dirac cone still exists, but is visibly smaller. The corresponding traced dispersion (Fig. 3j) shows the same kinklike feature around E~40 meV in σ-pol imaged band structure, but with a smaller amplitude and an onset slightly closer to the Dirac point. The momentum offset between σ-pol and π-pol dispersions is reduced to roughly 60% of the value before annealing, and is less prominent outside of a narrow energy window from E5 0-70 meV above the Dirac point. A final experiment followed 2 more hours of LT annealing (3 h in total), and found almost no difference between the π-pol and σ-pol ARPES spectra (Fig. 3g). Unlike the first 1 h LT anneal, very little change is noted in the sharpness of the bands following this final 2 h LT anneal. The traced bands (Fig. 3h) show the same dispersion under both polarizations, and there is no such anomalous kink-like feature. Though changes in band dispersion are of particular interest for physics and applications, a more basic property of impurity resonance known from STM and theory is the build-up of a DOS peak at the resonance energy. [20][21][22][23][24][25][26][27][28][29] This DOS feature cannot be identified from the ARPES band dispersions, as Luttinger's theorem does not apply to disordered systems. 36,37 To evaluate DOS, we instead sum the ARPES spectral function over the 2D momentum space, making use of the continuous rotational symmetry near the Dirac point to define DOS(E) = P k 2π k j j Ã IðE; kÞ. This symmetrization procedure is applied separately to the σ-pol and π-pol dichroic ARPES data to obtain differently weighted approximations for the surface state DOS (Fig. 4a). In all cases, the π-pol imaged DOS curve has a linear trend in the upper Dirac cone, matching expectations for a massless 2D Observation of a topological insulator Dirac cone L Miao et al.
Dirac fermion. However, the σ-pol data reveal an extra hump around E~40 meV in the non-annealed base sample and after 1 h of annealing. To exclude the possibility that this feature is from an extraneous matrix element effect, a corresponding simulation for a pristine Dirac cone is shown in Fig. 4b based on the energyresolved ARPES matrix elements identified in ref. 38 . The symmetrized DOS curves for both measurement conditions show a good linear character, and neither has an anomalous peak feature, suggesting that the extra hump in the symmetrized DOS curves represents a real anomaly within the Dirac surface state. A cleaner view of the dichroic feature is obtained by subtracting the π-pol curve from the σ-pol DOS (Fig. 4c), which reveals the anomalous hump as a peak at E~40 meV in the dichroic DOS curve of the non-annealed sample. This peak continuously loses intensity with LT annealing, and is no longer identifiable after 3 h of annealing. The PR-weighted impurity resonance peak simulated in Fig. 1d is reproduced at the bottom of Fig. 4c, and is energetically consistent with the dichroic DOS feature. This correspondence can only be qualitatively meaningful as it is filtered through ARPES matrix elements; however, it indicates that the kinked band dispersion observed in σ-pol ARPES spectra is appropriately aligned with the impurity resonance in exactly the way that was theoretically predicted.

DISCUSSION
These results show that Bi 2 Se 3 samples with a high point defect density can exhibit significant changes in surface state dispersion and DOS, both of which are reversible via an aging process that includes LT annealing. Numerical modeling closely reproduces the experimental features, and suggests that the dispersion kink and higher energy momentum shift in σ-pol ARPES band structure originate from a small discontinuity in the Dirac bands due to hybridization with impurity resonance states. While this may seem like a disruption of the topological band connectivity, it does not involve time-reversal symmetry breaking, and should not necessarily be viewed as a band gap. Band topology constraints are absent in the case of strong disorder, and the DOS distribution is peaked, rather than gapped, at the resonance energy.
With respect to earlier studies, it is noteworthy that the π-pol data do not show a similar kink-like dispersion feature, suggesting that the measurement conditions used for most ARPES experiments will give little weight to impurity-rich surface regions (see Supplementary Note 1). Moreover, the in-plane σ-pol orientation chosen for this study was selected to achieve particularly favorable matrix elements for resolving impurity physics (see Supplementary Note 2).
Unlike the conventional manifestation of a non-dispersive impurity band, the impact of disorder is seen to play out over a wide energy scale of several hundred millielectronvolts. Together with the theoretical protection against Anderson localization, 28 this suggests that the Dirac cone electrons have a unique coherent relationship with the disordered impurity lattice, which preserves reasonably sharp quasiparticle-like character in spite of allowing large changes to the overall electronic structure. The DOS peak found at the dispersion kink is consistent with real space STM investigations of impurity resonance, and has been predicted to greatly enhance the susceptibility to magnetic order in mean-field modeling for an ordered impurity lattice. 27 Further theoretical work will be needed to more fully understand the many-body susceptibilities of a dispersive quasiparticle-like mode that is profoundly inhomogeneous on a <100 nm scale, and for which the number of states in a "band" cannot be treated as constant throughout momentum space (i.e., for which Luttinger's theorem does not hold). These results and predictions suggest that the non-magnetic impurities found ubiquitously at the surfaces and buried interfaces of real-world samples and devices may provide a far-reaching mechanism for shaping the physical properties of a TI surface.

High-resolution linear dichroic ARPES measurements
All ARPES measurements were performed at the BL4.0.3 MERLIN ARPES endstation at the Advanced Light Source, with a Scienta R8000 analyzer and base pressure better than 5 × 10 −11 Torr. The sample was maintained at T~20 K, the time between s-polarization and p-polarization measurements is roughly 4 h, and surface band structure features were observed to be stable within~10 meV on the time scale of the each LD-ARPES experiment (~<24 h). The chemical potential for all samples was E-E D3 00 meV above the Dirac point, which matches the expectation of E-E D = 325 meV above the Dirac point for the Hall effect carrier density of 2.6 × 10 19 cm −3 . This surface potential is calculated using the band bending model in ref. 39 , which treats the surface state and bulk charge carriers on equal footing in the modified Thomas-Fermi approximation.
Surface doping of Bi 2 Se 3 by adatoms and photon exposure is a natural concern in quantitative ARPES experiments, and multiple doping mechanisms can come into play. 19,33,[39][40][41][42][43][44][45][46][47] However, the observed stability of the measured surface in this particular case is consistent with strong bulk screening associated with the high bulk carrier density, and with the observation that aging effects from residual gas and photon exposure tend to saturate at lower surface potentials of E-E D~< 290 meV. 40 The overall energy resolution was between 15 and 20 meV, and photon flux was well below the regime on which photo-gating effects have been observed. 41 Measurements on the base sample were performed at hv = 30 eV, and later post-annealing measurements were performed at hv = 34 eV to obtain a higher photon flux, due to the rapid loss of photon throughput on the high-energy grating beneath 50 eV.
Growth and STM characterization of defect-enriched Bi 2 Se 3 samples High-quality defect-enriched samples were grown following the procedure described in ref. 31 . The resulting crystals show no sign of impurity phases.  Fig. 4 The impurity resonance DOS peak. a ARPES DOS curves estimated from (red) σ-pol and (blue) π-pol are shown for different levels of LT post-annealing (PA). The curves are offset in increments of 2, and normalized at E = 120 meV above the Dirac point (dashed line), an energy that is higher than (shaded region) the expected impurity resonance and low enough to avoid the bulk conduction band. b Simulated ARPES DOS curves for a pristine Dirac cone, based on empirical photoemission matrix elements from ref. 38 . c Dichroic DOS curves the subtracting the π-pol ARPES DOS from σpol ARPES DOS. The dichroic DOS curves are compared with (bottom) the PR-weighted impurity resonance simulation from Fig.  1d. The anomalous peak-like feature that vanishes with annealing is indicated with filled triangles Observation of a topological insulator Dirac cone L Miao et al.
The Hall effect carrier density was 2.6 × 10 19 cm −3 , corresponding to an expected defect density of ∼0.19%, which corresponds reasonably well with the roughly~0.12% defect count in the top nanometer of the crystal seen by STM topography maps described in the main text. ARPES and STM measurements were performed close to the sample center, where the defect density is expected to be lower than on the periphery. The Hall carrier density and surface chemical potential are not greatly changed following the aging/annealing process A slight reduction in the surface chemical potential suggests that the bulk carrier density is reduced by 10-20% in later measurements, relative to density beneath the beam spot in the initial measurement. Additionally, it is expected that the postannealed defect distribution may be more homogeneous across microcrystalline domains.
STM data were obtained at LT (T < 50 K) with tips calibrated on an Au (111) surface. The topographic map in Fig. 2b was measured using a bias of −0.7 V and a tunneling current of 100 pA. Rapid topographical maps of multiple~100 nm regions were sampled, and were consistent with previous analyses that suggest a stochastically random placement of impurities within local regions (see the Supplemental material in ref. 28 ). However, the standard deviation in defect density for different surface regions separated by several microns was much higher than the variation expected from the Poisson distribution. This variation reveals large fractional differences (up to a factor of~2) in the local impurity densities that are averaged over in the ARPES spectral function.

Modeling the TI surface
The surface state is modeled as a spin-helical 2D Dirac cone on a hexagonal lattice, perturbed by scalar delta-function-like impurities, as described in ref. 28 . The impurity potential is effectively reduced by the fact that the delta-function-like potential is not well resolved in the model, which imposes a high-energy cutoff on the kinetic basis for state diagonalization. An effective value for the potential is calculated as the trace of the potential Hamiltonian for a single impurity (U eff = Tr(H U ) = −1 eV), whereas the uncorrected potential would have a value of U = −35 eV, as defined in ref. 28 . The impurities are randomly distributed with a density of ρ = 0.06%. Spectral features were convoluted by a Lorentzian function with 30 meV peak width at half-maximum, except where otherwise noted.
Intensities in Fig. 1c are weighted by the PR P of each single-particle eigenstate. The PR gives higher intensity for states that are less evenly distributed throughout space, and is used to highlight emission from defect resonance states. It is defined as: where the sum is overall sites in the system, and n α,i is the local DOS on site i of |α〉, which is an eigenstate of the full Hamiltonian.

Data availability
The data and source code that support the findings of this study are available from the corresponding authors on reasonable request.