Disorder enabled band structure engineering of a topological insulator surface

Three-dimensional topological insulators are bulk insulators with Z2 topological electronic order that gives rise to conducting light-like surface states. These surface electrons are exceptionally resistant to localization by non-magnetic disorder, and have been adopted as the basis for a wide range of proposals to achieve new quasiparticle species and device functionality. Recent studies have yielded a surprise by showing that in spite of resisting localization, topological insulator surface electrons can be reshaped by defects into distinctive resonance states. Here we use numerical simulations and scanning tunnelling microscopy data to show that these resonance states have significance well beyond the localized regime usually associated with impurity bands. At native densities in the model Bi2X3 (X=Bi, Te) compounds, defect resonance states are predicted to generate a new quantum basis for an emergent electron gas that supports diffusive electrical transport.

rapidly than the band slope for momenta k 0.05Å, however it is difficult to trace the band slope with sufficient accuracy in this region.
Histograms of the modeled participation ratio are compared with DOS and band structure calculations in Supplementary Figure 4, with dashed lines bracketing the electronic structure most influenced by defects. Summing all states between the dashed lines gives a number equal to twice the number of defects, showing that approximately 2 resonance states with large participation ratios can be identified per defect. This is the behavior that would be expected if the resonance states were energetically-isolated bound states, as all bound states are doubly degenerate in the absence of time reversal symmetry breaking.

SUPPLEMENTARY NOTE 2: ROTATIONAL ANISOTROPY OF DEFECTS
The STM image of a defect resonance state in Fig. 1a has 3-fold rotational symmetry, with 3 bright spots and 3 dark spots. At a fixed radius, the ratio of intensity maxima to minima is approximately R = Imax I min 1.5, and systematic error in the background may add up to 0.2 to the ratio. Defining normalized basis states for for s-and f-wave components of the wavefunction as Ψ s (θ) = 1/ √ 2πr and Ψ f (θ) = sin(6θ)/ √ πr, the ratio of maximum to minimum intensities can be calculated as obtained with an f-wave partial density of states of A f 2 = 0.01, giving an upper bound of roughly 1% for the f-wave admixture.
Hexagonal warping of the Fermi surface is not included in this study, as it scales with the third power of momentum [7], and the present focus is on small momenta close to the Dirac point. With realistic parameters (λ 250eVÅ −3 for Bi 2 (Se/Te) 3 [7]) the effect of hexagonal warping is a very slight increase to the f-wave component of resonance states, to a degree that is not easily visible on the color scale in Fig. 1(a) of the main text.
Momentum-resolved spectra are highly isotropic, and have been rotationally averaged in the main text to provide a more continuous map of momentum dependence in the ARPES spectral function. When plotted along a single momentum axis, the features are effectively identical, but are broken into broadly spaced discrete momentum states that intersect the chosen momentum axis. An example of what this looks like for the ARPES simulations in The energy cutoff of W=400 meV is chosen to roughly match the range over which a Dirac-like dispersion is known to be realized in Bi 2 Se 3−x Te x . This cutoff means that the basis cannot nicely resolve the atomic structure of defects. In the main text, point defects are created as 3-site clusters with a perturbation strength of U/3 on each site. However, changing this to a larger 6-site triangular cluster (U/6 on each site) or a single atomic site perturbation has no qualitative impact on simulted spectra, as seen in Supplementary Figure   6. The value of U required to maintain the resonance state energy relative to the Dirac point is slightly smaller for larger clusters, as described in the figure caption. However, setting the resonance state energy to match the experimental scenario for defects discussed in the main text results in very similar spectra, which are just past the crossover point of manifesting an intensity local maxiumum at the D 1 Dirac point.

SUPPLEMENTARY NOTE 3: ADDITIONAL MODELING DETAILS
Uniformly distributed random numbers used to generate the defect configurations were obtained using a Mersenne twister generator [8]. We note that the distribution of defects observed on real Bi 2 Se 3 surfaces is also highly uniform and uncorrelated on a ∼100nm scale, as seen from the Poisson fit in Supplementary Figure 1b. Full diagonalization of the Hamiltonian was performed using LAPACK drivers [9], with numerical accuracy at least 5 orders of magnitude better than the nearest significant (or visually resolvable) digit of extracted quantities. Individual simulations of system sizes of 300×300 sites or greater yielded visually smooth and reproducible ARPES spectra. Smaller systems are only considered in Fig. 4 of the main text, where configurational averaging over multiple defect distributions is used to achieve smooth trend curves for the participation ratio and twist velocity.
The simulation in Fig. 3f is seeded from defect locations observed by STM in Supplementary Figure 1a, rather than a random distribution. A 10nm border was cropped from the edge of the simulation output, to reduce inaccuracy associated with the repeating boundary conditions of the model. Unfortunately, the mapped region has a surface area more than 50% smaller than the 300 × 300 system size qualitatively associated with good convergence.
The small system size contributes to roughness in the LDOS curves and is problematic for extracting the LDOS curve minima, which are intrinsically unstable against noise. To obtain a reliable spatially-resolved estimate of the LDOS minimum, LDOS curves were convoluted by a 20meV peak width at half maximum (pwhm) Gaussian, and the local minima were averaged over an 8nm radius. The 16nm diameter of averaging is similar to the wavelength of electrons roughly 200meV from the Dirac point, and was chosen to maximize accuracy while remaining near the resolution limit of real physical effects near the Dirac point energy.