Connection topology of step edge state bands at the surface of a three dimensional topological insulator

Topological insulators in the Bi$_2$Se$_3$ family manifest helical Dirac surface states that span the topologically ordered bulk band gap. Recent scanning tunneling microscopy measurements have discovered additional states in the bulk band gap of Bi$_2$Se$_3$ and Bi$_2$Te$_3$, localized at one dimensional step edges. Here numerical simulations of a topological insulator surface are used to explore the phenomenology of edge state formation at the single-quintuple-layer step defects found ubiquitously on these materials. The modeled one dimensional edge states are found to exhibit a stable topological connection to the two dimensional surface state Dirac point.


Introduction
Three dimensional topological insulators (TI) are materials with Z 2 topological order that manifest conducting two-dimensional (2D) Dirac cone surface states protected by time-reversal symmetry [1,2]. The surface state electrons resist scattering from weak non-magnetic perturbations [1][2][3][4][5][6][7][8][9], however recent studies have shown that strongly perturbing point-and step-like surface defects can introduce new in-gap states and modify the band structure near the Dirac point [10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26]. Here, we present a numerical and analytic analysis of the singlequintuple layer step defect of Bi 2 Se 3 -family TIs, to explore the nature of the associated edge state. These investigations establish that edge states with a stable topological connection to the 2D Dirac point of the TI surface state can exist, in scenarios that do not require fine tuning of the Hamiltonian. As with the surface states of Weyl semimetals, the connection is defined by a linear dispersion of the lower-dimensional (here 1D) state, converging on the Dirac point of higher dimensional (2D) band structure. The protected nature of this connection is considered with respect to broken symmetries and disorder, and it is shown that the occurrence of such states can be ubiquitous across a wide range of parameters for describing surface steps in a 3D TI tight binding model. Bi 2 Se 3 is widely seen as a model system for studying TI surface physics, as it has one of the largest band gaps presently known in a TI system (∼300 meV), which is spanned by a relatively ideal single Dirac cone surface state [27][28][29]. The crystal structure of Bi 2 Se 3 is shown in figure 1(a), with weak Van der Waals bonding between stacked quintuple atomic layers.
Step defects one quintuple layer in height are very common at the surface of thin film samples in this material as shown in figure 1(a) [12,30]. Figure 1(b) shows the arrangement of Se atoms at a cleaved surface near a single step, seen from above. Viewed on a micron scale, the steps tend to run parallel to the in-plane nearest-neighbor axis, resulting in triangular plateaus as shown in figure 1(c) [12,30].
Protected 1D edge states are typically associated with 2D topological order, however several classes have also been proposed to occur at the surfaces of 3D material systems under specific circumstances. Hinge states can occur at the intersection of two non-parallel faces of a so-called 'higher order topological insulators' [31,32], and are expected in special cases for traditional 3D TIs [33]. In the 3D TI case, if the intersecting faces have Dirac points at the same energy, the edge can host an in-gap mode that intersects the 2D Dirac point, associated with a near-realization of the Jackiw-Rebbi Hamiltonian [34]. However, structurally simple intersections of this type are challenging to create and study experimentally. Another topologically protected edge state scenario has been identified at certain classes of step edges on a topological crystalline insulator surface, and can be understood from extrapolation to a scenario in which particle-hole symmetry is unbroken [35].
Other 1D in-gap bound states that converge on the Dirac point of a massless 2D Dirac Hamiltonian have been identified for specific models, such as the bound states underneath a gate electrode [36] or a 1D Gaussian potential [37]. These scenarios are intriguing because although they do not require fine tuning of the Hamiltonian, they nonetheless appear to be almost coincidental, and have not been identified explicitly with a distinct topological invariant of the system. The principal result of this paper will be to show that this class of topologically connected edge states is insensitive to the specific form of the TI Hamiltonian or the 1D feature to which edge states are bound, and is expected to be a generic feature of line-like defects on the surface of 3D TIs (or, equivalently, planar defects in a 3D Weyl kinetic Hamiltonian). This paper is organized as follows: in section 2, we establish a broad set of conditions allowing the existence of edge states with non-trivial band connectivity, bound to a 1D scalar potential. In section 3, we present numerical simulations on a 2D lattice model motivated by real TI surfaces, and demonstrate the robustness of the connectivity of 1D edge states against symmetry breaking. In section 4, we show examples of this topological connectivity within a 3D tight binding model resembling Bi 2 Se 3 . (c) Large scale schematic of the triangular step plateaus found on Bi 2 Se 3 . The image is not drawn to scale, as real terraces tend to be over 100 nm in size.

Allowed existence of an edge state with non-trivial band connectivity
Bound states at quintuple layer step edges have been observed for the related TIs Bi 2 Se 3 and Bi 2 Te 3 [16,38], and a recent STM analysis has been interpreted to suggest that they may be a form of 1D electron gas brought on by an effective scalar potential at the step [38]. The in-gap surface states of these materials have large wavelengths (λ>5 nm), meaning that the real-space structure of the step edge potential is not well resolved on this basis. Because of this, and for the sake of analytical tractability, we will begin by considering the scenario of a 1D delta function potential added to a 2D system defined by a massless Dirac Hamiltonian: where v D is the Dirac velocity,k is the momentum operator and σ is the Pauli vector. The scalar potential is parametrized by a constant U with units of energy. An additional constant w with units of length can be thought of as the real-space width of the step edge potential. Momentum along the 1D potential axis is a conserved quantity (k P ), and the 1D bound state dispersion in this scenario is already known from [36] to converge on the 2D Dirac point with a constant group velocity. Changing the delta function potential strength modifies the velocity of the bound state, but not the momentum space connectivity. The connection to the 2D Dirac point can thus be called a topological property, in that it is robust for a continuous range of 1D perturbation strengths. However, the edge state connectivity is not generically protected from Coulomb perturbations, as a surface potential in the right form could in principle cause the group velocity of the edge state to exceed the Dirac velocity, pushing the edge state dispersion into the 2D Dirac cone, making it no longer a well defined edge state (dashed lines in figure 2(a)).
The projection of the bound state onto kinetic basis states in the 2D Dirac cone is localized, and is negligible at energies outside of a > L  | | E v k window around the 2D Dirac point, for sufficiently large values of a constant L v (see discussion in appendix A). We observe that this has the consequence that a non-delta function potential will generally give a very similar linear dispersion converging on the Dirac point (at k P ∼0). So long as the Fourier transform of the potential (V(q)) is nonzero at q=0 and effectively flat on a momentum scale of  d L  q v k v D , the potential will be indistinguishable from a delta function within the low energy basis that composes the step edge state. This scenario is played out in a practical context in [37], which establishes the linear dispersion for a Gaussian 1D potential. Qualitatively speaking, as the amplitude of k P is increased, nonlinearities in the surface state dispersion will emerge when the curvature of V(q) becomes significant on a momentum scale proportional to k P , and/or when the kinetic state basis at energies within < L  | | | | E k v deviates from a massless 2D Dirac Hamiltonian (see turning point indicated by arrows in figure 2(b)). Corrections to the edge state energy from these factors act at lowest order in proportion to k P (see appendix A), meaning that they can influence the edge state velocity, but not the point of convergence.
Deviation from a massless 2D Dirac Hamiltonian is inevitable at large momentum in a real material, and must be accounted for to understand the band structure connectivity of the other end of the 1D edge state in momentum space. If the 2D Dirac cone Hilbert space is curtailed by a high energy cutoff such that kinetic basis states outside a momentum window < L | | k are disregarded, the edge state will follow an arching dispersion like that drawn in figure 2(b) (see also an analogous simulation in figure 3(a)). This scenario is close to what is expected at a real step edge, since the 2D Dirac cone surface states of TIs are only well defined over a small energy range inside the bulk band gap. Once the band momentum becomes sufficiently large, the edge state of a delta function potential curves in a direction opposite to the sign of the edge potential (U, defined below). It is required to merge with the 2D Dirac cone at or before the momentum cutoff, due to the lack of states for the 1D potential to couple between as momentum along the edge approaches the cutoff. This connection bridged by the edge state between the 2D surface Dirac point and the state continuum immediately above or below the 2D Dirac point is a further property that can be used to classify an edge state, and will be discussed later in the context of a more realistic model with full 3D topological order (section 4). The diagram in figure 2(c) shows a summary of the topological and non-topological dispersions that can be expected for a 1D edge state bound to a scalar potential parametrized by U. At a critical value of U=U C , the edge state will converge towards the Dirac point with a velocity v S equal to the 2D Dirac velocity (v S =v D ). Increasing U causes the dispersion to intersect with the 2D Dirac continuum (white region), so that the edge state does not include a well defined dispersion that intersects with the Dirac point. For a continuous range of lower potentials U<U C , the edge state converges on the Dirac point (red region, labeled Topological) until velocity of convergence becomes equal to the negative Dirac velocity (v S =−v D , not shown in figure 2(c)). The allowed existence of a topological connection that is protected over a finite range of constant prefactors for essentially any spatial form of the 1D potential is robust, and is not conditional on symmetries that do not destroy the Dirac point, such as positive reflection symmetry across the step, or particle-hole symmetry in the kinetic Hamiltonian. Numerical simulations in the next sections will show that Dirac point-intersecting 1D bound state dispersions are a common, and possibly ubiquitous, feature in typical models of step edges at TI surfaces.  (A7)). The potential used to describe the step edge is U=3.8 eV, as fitted to STM on Bi 2 Te 3 in our previous work [25].
shown as a function of the step potential (U) in the discrete lattice simulation and the analytic 2D continuum-limit theory. (d-f) LDOS showing step edge state dispersion for systems with U = 3.8 eV and broken reflection, particle-hole, or time-reversal symmetry, respectively.
3. Edge states of a 1D scalar potential in a 2D lattice model To relate the above picture more closely to the step edge states seen in real TI materials, we consider a 2D hexagonal real-space lattice that resembles the Bi 2 (Se/Te) 3 surface, following the modeling implementation in [25]. The step edge is described as a scalar potential that repeats along a chain of surface sites extending along the crystallographic a-axis (nearest-neighbor direction). The modeled system has translational symmetry along this axis, and is simulated with a large (effectively infinite) number of sites along the y-axis. The single-electron Dirac Hamiltonian for this discretized scenario can be written as where the kinetic Hamiltonian  0 is unchanged from equation (1), and the potential term  U is described in terms of the one-site number operator n α acting on sites α intersected by the step. The Dirac velocity is taken to be 3 eV Å, the hexagonal lattice constant is 4.2 Å, and energy cutoff for the kinetic basis is Exact diagionalization is used to obtain the eigenstates and energies of the system.
The local density of states (LDOS) on top of the modeled step edge was shown in our previous work [25], and closely matches the momentum-integrated LDOS profile of a step edge state seen by STM on Bi 2 Te 3 , at a step potential U=3.8 eV. The dependence of step edge LDOS on momentum parallel to the step (k P ) is shown in Increasing the step potential causes the edge state to converge towards the upper Dirac cone dispersion (figure 3(b)-(c)). State dispersions closely resemble the continuum-limit (CL) expectation for a delta function potential, derived in appendix A. In this idealized case, the edge state velocity approaches the Dirac velocity as  ¥ U , meaning that there is no finite value for U C , and the system has topological connectivity for all values of U. The CL model includes just one surface state, and large deviations from the model are expected to occur when a new surface state appears, as will be shown in the 3D tight binding model simulations in section 4.

Protection of the edge state
The allowed existence of edge states with a protected connection to the Dirac point as established in section 2 is not conditional on Hamiltonian symmetries such as reflection symmetry or particle-hole symmetry, so long as the 2D Dirac point itself is not destroyed, and the Fourier transform of the 1D potential does not vanish at q=0. For example, reflection symmetry of the 1D potential is broken in figure 3(d) by parallel lines with U=0.1 eV and U=−0.1 eV one lattice site above and below the U=3.8 eV 'step', respectively. Similarly, the simulation in figure 3(e) breaks particle-hole symmetry of the kinetic Hamiltonian by reducing the negative energy cutoff to −0.3 eV. In each case, the Dirac point connectivity of the edge state is unchanged, but the velocity with which it intersects the Dirac point is slightly altered.
However when the 2D Dirac point itself is gapped by the addition of a Zeeman Hamiltonian term that breaks time-reversal symmetry, the protected connection is also necessarily broken as shown in figure 3(f). In this case, the full Hamiltonian is where B z defines the strength of a magnetic exchange field oriented normal to the TI surface. The edge state disperses through the gap with a new extremum at zero momentum that would create an additional square root anomaly in sufficiently high-resolution LDOS measurements. The precise dispersion can be extrapolated by noting that the Zeeman and k P terms in the model Hamiltonian combine to create a Pauli vector that is orthogonal to the only other Pauli vector in the Hamiltonian (which comes from k ⊥ ). As such, the edge state energy at momentum k P in the presence of a perturbing magnetic field will be identical to the unperturbed energy at a different momentum ¢ Disorder is another phenomenon that can render topological band features ill-defined, by causing momentum to no longer be a good quantum number. However, because the edge states appear in spin-chiral time-reversal pairs, disorder that does not violate time-reversal symmetry will only mix the edge states with bulk states, and not with the time-reversed partner. The effect of time-reversal invariant disorder is therefore expected to be proportional to the bulk density of states, which vanishes at the Dirac point. While the states may still be destroyed due to strong disorder, the connection to the Dirac point has qualitatively more protection than dispersions at other energies, suggesting that it will be a robust phenomenon.

Physical steps in a 3D tight binding model
To explore the phenomenology of the step edge state connectivity in a less idealized scenario, we briefly examine single-and double-quintuple layer step edges in a minimal 3D tight binding model. Bulk and surface state kinetics in the model are similar to vacuum-cleaved Bi 2 Se 3 , and it incorporates the same TI topological invariants ([1; 000]) derived from a band inversion at the 3D Γ-point. In this approach popularized by [2], a single Bi 2 Se 3 formula unit is reduced to two p z orbitals displaced along the surface normal axis. Parameters of the model are listed in appendix B, and the projection of 2D momentum space onto the 1D momentum axis parallel to the step is shown in figure 4(d). The step aligns with the 2D Γ-K axis, and creates a new periodicity from Ḡ to Ḡ which is traced in red in figure 4(d), with Kramers degeneracy required at theM point.
In this case, incorporating the physical one quintuple layer step drawn in figure 1(a) involves modifying only the kinetic Hamiltonian, but nonetheless results in the appearance of an edge state that connects between the Dirac point and the lower state continuum (see figure 4(a)). Unlike other effectively 1D topological in-gap states, the edge state shown in figure 4(a) does not depend on extrapolation to a scenario with particle-hole symmetry, and actually vanishes when particle-hole symmetry if restored to the Hamiltonian (by setting h 0 =0; see appendix B). Increasing the height of the step to 2 quintuple layers results in a larger edge state group velocity near the Dirac point ( figure 4(b)), but does not change the Dirac point or lower continuum connectivity. We note however that a sufficiently large step edge will effectively introduce a new 2D state continuum that overlaps with and obscures the 2D Dirac point.
Turning on a positive scalar potential (U>0) acting at nearest-neighbor sites bordering the step causes the edge state group velocity to increase at the intersection point with the Dirac point ( figure 4(e)), matching the behavior predicted in the figure 2(c) phase diagram. The edge state merges with the upper Dirac continuum above a critical potential of U>U C ∼0.3 eV, and ceases to connect to the Dirac point. However, rather than the system entering an extended 'non-topological' phase region as posited in figure 2(c), a new Dirac-connected surface state emerges as an antibound state of the lower Dirac cone. This new state is indicated by a green arrow in figure 4(c), and is the Kramers partner of the original step edge state at theM point. Tracing the dispersion, we see that these two edge states effectively connect between the 2D Dirac point and the upper band structure continuum, whereas the original state connected to the lower continuum. This topological distinction is indicated by shading in figure 4(e).
Comparing with the analytic CL solution for a 2D Dirac surface with a delta function potential (CL model, see appendix A), we find that the analytic model gives a close match for the low-momentum dispersion as the step edge states emerge from the lower Dirac cone, but diverges as they approach the upper Dirac cone at V S /V D >0.5 in figure 4(e). For this plot, the analytic curves have been shifted to align with the tight binding model on the U-axis, and the input U parameter has been rescaled upward by a factor of 3.9 for the left curve and 2.4 for the right curve. Taken together, this large upward rescaling of the U-axis, as well as the non-infinite critical potential for changing the Dirac point connection topology ¹ ¥ ( ) U C , and the appearance of successive Dirac-connected states as a function of U, reveal that including 3D structure and coupling with bulk band symmetries can have an important role in defining the edge state properties. This finding is consistent with other recent literature suggesting that impurity state LDOS may include a significant bulk band component at a TI surface [14,26,39,40]. However, the rescaled 2D CL model can do an excellent job for edge states with group velocities not far removed from their band of origin (the lower or upper Dirac cone), and accurately describes the low-momentum dispersion for over half of the parameter space explored in figure 4(e).

Summary
In conclusion, we have established that step edge states in single-particle models of a 3D TI surface can manifest topological phase regions featuring protected connections to the 2D surface Dirac points, with no reliance on fine tuning in the bulk or surface Hamiltonian. This analysis builds on the previous observation of Dirac point connectivity in specific 2D models, and provides a guiding principle for understanding the likely form of in-gap states observed by STM at step edges. Realistically parametrized simulations of a Bi 2 (Se/Te) 3 single-quintuple layer step edge are performed using a 2D Dirac cone Hamiltonian and a 3D tight binding model, and establish that topologically connected step edge states are physically plausible in this material family.

Appendix A. Exact CL model
In this section, we solve the exact continuum model for a 2D Dirac cone surface with a delta function potential, and a high energy cutoff. Imposing a high energy cutoff on the kinetic basis is shown to have no effect on bound state dispersion near the 2D Dirac point. Complementary derivations for very similar scenarios can be found in [36] and [41], and the analytic real-space wavefunction derivation in [41] can be manipulated to identify the localized E −2 decay trend of the DOS projection onto high energy kinetic eigenstates.

A.1. The gapless Hamiltonian
Consider a (2+1)−D Dirac fermion in the presence of a singular delta potential parallel to the y-axis. In units where ÿ and v F are unity, the effective Hamiltonian is Here, W has units of energy times length, and is equivalent to the U step potential parameter defined in equation (2) multiplied by the width of the potential barrier. Since translation symmetry exists along the y-axis, the eigenstates can still be labeled by the eigenvalues ofp y . For a given eigenvalue p y , we are thus left with a onedimensional problem of a Dirac fermion in the presence of a delta function potential at the origin.

A.2. The bound state energy: poles of the T-matrix
The energy of the bound state,  ( ) p b y , can be obtained by finding the real poles of the T-matrix of the problem, defined via is the retarded Green's function for the potential-free problem. For the Delta function potential, the T-matrix has a simple form . A4 x x x y y 0 1 The momentum cutoff, Λ, has been introduced above. The analytic structure is simple when E lies in the spectral gap, < + ( ) E p p  The bound state energy, given by the real poles of t(E) within the spectral gap, is found via the condition The solution to this is given by This needs to be solved numerically to obtain the bound state dispersion at a finite Λ. In order to maintain approximate rotational symmetry, Λ should be replaced by Lp y 2 2 . Some results are shown in figure A1.
Taking the first correction due to finiteness of Λ, we have Following the procedures outlined above for the gapless case, we obtain the bound state energy to be This needs to be solved numerically to obtain the bound state dispersion at a finite Λ. In order to maintain approximate rotational symmetry, Λshould be replaced by L -p m y 2 2 2 . Some results are shown in figure A2.
A.5. The bound state dispersion: analytic results for the gapped case Taking L  ¥ yields a scaled relativistic dispersion: Taking the first correction due to finiteness of Λ, we have We adopt the minimal tight binding model framework for topological order, in which a single Bi 2 (Se/Te) 3 formula unit within one quintuple layer of the crystal is reduced to two p z orbitals displaced along the surface normal axis [2,42]. These two-orbital unit cells are arranged in an AA stacked hexagonal lattice, with an in-plane lattice constant of a=4.2 nm and a z-axis lattice constant of c, the precise value of which is not physically relevant.
Taking into account of the spin degree of freedom, the four states are   ñ  | ( ) P , z , where±and   ( ) indicate the parity and spin of the state, respectively. To simulate a surface and step edge with translational symmetry along the x-axis, the modeling basis includes N z =10 z-axis layers, which was sufficient to decouple the top and bottom surfaces. Along the +ˆ( ) ) x y 2 3 2 nearest-neighbor axis inside the plane of the surface, N y =150 sites were used for full Brillouin zone spectral function maps, and a larger N y =1000 system was used to eliminate finite size effects in all other panels. Coupling around the in-plane axis resulted in a small Dirac point gap that scaled approximately as -N y 1 , and had a value of 3 meV for N y =1000. To avoid finite size effects in the analysis, the smallest nonzero momentum value considered in the manuscript has an amplitude of 0.005 Å −1 , giving an intrinsic kinetic gap in the ideal 2D Dirac cone that is a factor of 6.7 ( = D ( ) 6.7 3 meV T ) larger than the finite size gap. The displayed plots at this momentum amplitude were qualitatively insensitive to ∼50% fractional changes in system size, which is understandable as for a gapped 2-state system, an off-diagonal matrix element with this relative amplitude would account for just a 1.1% change in the gap energy, and 0.55% mixing of PDOS.
A z-axis slip in the linkage of repeating boundary conditions is used to create the surface edge, as shown in figure B1. A Coulomb potential U term was applied on the site closest to the step, to account for energetic factors that cannot be described by the simple implementation of z-axis slip. This is described by a Hamiltonian term = å ( ) H U n U i i , where n i is a number operator, and the sum over i indexes all four orbital and spin states in the unit cell with an open in-plane boundary (i.e. the 2 real-space orbitals connected by the A 1 dashed line in figure B1(b)). Though the top and bottom surfaces were effectively decoupled, an identical H U term was applied to the step edge on the bottom of the slab for the sake of symmetry. When represented in real-space, the hopping Hamiltonian term around the boundary includes a phase factor proportional to k P , due to the fact that the plane of the slab is not orthogonal to the step axis (x-axis).
In this model, there are two types of hopping-intra-layer and inter-layer hopping terms of which the vectors connecting the corresponding unit cells are represented by  a 1,2,3 and  a 4 : where A i and B i are the parameters for the intra-layer and inter-layer hopping terms, respectively. In our simulation, the parameters are given as following: With the above hopping parameters, the resulting model is a strong 3D TI with topological invariants [1; 000] and band inversion at the 3D Γ point as for Bi 2 Se 3 -family TIs. The bulk band gap is 0.3 eV, similar to that of Bi 2 Se 3 , and the Dirac cone has a velocity of ∼2 eV Å −1 .