Non-magnetic impurities and in-gap bound states in topological insulators

In-gap bound states induced by non-magnetic impurities in various dimensional topological insulators are investigated based on a modified Dirac model that considers quadratic corrections to the mass term. Their existence and features greatly rely on the potential form of the impurity as well as the dimensionality of the topological insulator. It is analytically proven that the impurity potential modeled by the delta function can induce the bound states in one dimension (1D), but not in two and three. For a single non-magnetic impurity with a general isotropic potential, formal solutions are obtained and further numerical calculations are performed. In particular, the in-gap bound states induced by a non-magnetic impurity with isotropic Gaussian potentials in two-dimensional (2D) and three-dimensional (3D) topological insulators are numerically investigated. Information on how many in-gap bound states can be trapped by a non-magnetic Gaussian impurity is presented for the parameters from a series of topologically non-trivial materials.

Besides these efforts on the random Anderson disorders, the influence of ordered or even single imperfection have also attracted growing attention in the community [24]. It was shown that in three-dimensional (3D) TI Bi 0.9 Sb 0.1 [25], dislocation lines (1D imperfections) are associated with 1D fermionic excitations, which are topologically protected and not scattered by disorder. In addition, it is already known that in many systems bound states can be induced by impurities, which are zero-dimensional (0D) imperfections. A series of works have been done on the issue of a single impurity on the surface of a 3D TI starting with the gapless Dirac model [26][27][28][29][30]. However, since the boundary states are only a manifestation of the topological nature of bulk bands, one must examine the host bulk to know how these imperfections affect the system. In our recent work [31], we explicitly presented the existence of the in-gap bound states trapped by a single vacancy in 2D and 3D TIs based on a modified Dirac Hamiltonian that describes the TI bulk very well [32]. Their possible influence on the transport features of TI was also discussed. The existence of in-gap bound states in that work is due to the vacancy 3 boundary condition (BC), which demands that the electron wavefunction must vanish in the vacancy region.
In contrast to vacancies, impurities manifest themselves as additional potential terms to the Hamiltonian of the system. Although the electron wavefunction can tunnel the impurity potential and does not have to entirely vanish where the impurity locates, the electronic transport may be more or less impeded. On the other hand, in modern semiconductor techniques, intentional doping is an efficient way of improving the performance of materials by modulating their electrical properties. Therefore for TI samples, whether to understand the existing data on their electronic and transport features or to find appropriate dopants for performance improvement, it is worth performing a systematic investigation on the issue of whether and how a single impurity can induce in-gap bound states. In the simplest case, the impurity is assumed to be non-magnetic, which does not break time-reversal symmetry (TRS). Furthermore, the impurity potential can be assumed to be isotropic, which is not far away from reality but greatly simplifies our discussions. In this work, based on the modified Dirac model [32], we report that, depending on the impurity potential form as well as the dimensionality of the TIs, bound states could be induced around non-magnetic impurities with their energies lying in the bulk energy gap.
This paper is organized as follows. In section 2, we briefly review the modified Dirac model. In section 3, formal solutions of the modified Dirac model in various dimensions are presented when a single non-magnetic impurity with the general potential form appears. In particular, some analytical results are obtained when the impurity potential takes a δ-function form. In section 4, in-gap bound states trapped by a single non-magnetic impurity with isotropic Gaussian potentials in 2D films and 3D bulk TIs are numerically investigated. In section 5, the possible detection of these in-gap bound states by using a scanning tunneling microscope (STM) is discussed. Finally, a summary is given in section 6.

Modeling
Following our previous work [31,32], a modified Dirac model is employed to describe the TIs in various dimensions in a unified way, where −Bp 2 = −B( p 2 x + p 2 y + p 2 z ) is the quadratic correction to the topological mass mv 2 , and p i = −ih∂ i , i ∈ {x, y, z}, is the momentum operator. v, m and B have the dimensions of speed, mass and reciprocal mass, respectively. The Dirac matrices satisfy the anticommutation relations representation of them can be expressed as where σ i=x,y,z are the Pauli matrices, σ 0 is the 2 × 2 unit matrix and ⊗ represents the Kronecker product. The Hamiltonian (1) preserves TRS: T H 0 T −1 = H 0 , where T = −iα x α zK is the time reversal operator andK is the complex conjugate operator.
The modified Dirac Hamiltonian (1) shares the same mathematical structure as those for 2D and 3D TIs [7,18]. The topological nature carried by it can be analyzed following the Z 2 classification of Kane and Mele [5]. The main results can be found in our previous work [31,32]. In a word, the modified Dirac Hamiltonian becomes topologically nontrivial (the Z 2 index is equal to 1) if m B > 0 and topologically trivial (the Z 2 index is equal to 0) if m B < 0. A similar conclusion is also drawn by Volovik in [33]. Due to the bulk-boundary correspondence, there always exist topologically protected boundary states at the open boundaries of a TI, where the Z 2 invariants change from nontrivial to trivial [34,35]. This feature can be well described by Hamiltonian (1) with open BCs (OBCs) when m B > 0. For the vacancy BCs, such as 2D circular holes and 3D spherical vacancies [31], the corresponding solutions are called 'bound states' since they are spatially confined around the vacancies. These bound states have the same origin as those due to the Z 2 classification.

In-gap bound state solutions in the presence of a single isotropic non-magnetic impurity
In this work, we focus on the effects of non-magnetic impurities. In the presence of a nonmagnetic impurity, a potential term V (r)σ 0 ⊗ σ 0 is added to Hamiltonian (1). As the impurity potential is assumed to be isotropic, one has V (r) ≡ V (r ). In this section, we present the formal solutions of the modified Dirac model in the presence of V (r ).

Integral equation for bound state energies
Before directly solving the model, the bound state energies can be formally obtained through an integral equation. Although in most cases this integral equation cannot be solved analytically, it does provide rich information about the existence of bound states under certain impurity potentials in various dimensions.
The wavefunction (r) can be expanded by its Fourier transformation components as (r) = p u p e ip ·r/h . Thus one has where V pp = drV (r) e −i(p−p )·r/h . Once the potential is given, V pp can be calculated. After putting it back to equation (4), in principle the eigenenergies and eigenwavefunctions can be solved. However, except for several special cases, for a general potential V (r), equation (4) is hard to solve analytically.
3.1.2. δ-potential. For a delta potential V (r) = V 0 δ(r), V pp ≡ V 0 , then one has A non-trivial solution to the four-component spinor ' p u p ' requires where d is the dimensionality. For the 1D case, the modified Dirac Hamiltonian can be easily inverted. After some algebra we have where E A and E B denote the energy solutions for '+' and '−', respectively. For the 2D case, one can obtain a similar integration equation for the 2D bound state energies, where k 2 = k 2 x + k 2 y . However, the integral in equation (8) will logarithmically diverge when |k| → +∞. This means that in 2D, an impurity with δ-potential cannot trap any bound states. Similarly, in 3D, although the integration equation is more complicated, divergence also exists in the k-integration, which excludes the possibility of 3D bound states under δ-potential. However, this will not always be the real case because of two facts. Firstly, in real materials with lattice structure, the momentum k cannot be infinitely large but confined within the first Brillouin zone. Thus, the integral in the 2D and 3D cases will not diverge and finite solutions of bound state energy are possible. Secondly, the potential form of a real impurity should be distinct from the ideal δ-function. Thus the disaster of integral-divergence in 2D and 3D could be avoided.
When an impurity with δ-potential V (x) = V 0 δ(x) appears in an infinite 1D TI, the decoupled 1D modified Dirac equation reads where (x) is a two-component spinor. The continuity of the wavefunction at x = 0 requires that ( ) = (− ), where is infinitely small. In addition, the integral of equation (11) around x = 0 gives Also, the electron wavefunction should vanish when x → ±∞. All these conditions lead to two transcendental equations for bound state energy: where up to two solutions can be found, which are just the E A/B defined in equation (7). Besides, symmetric properties between E A and E B can be observed from the transcendental equation as  (10). This double degeneracy comes from TRS. The energies of these states lie inside the bulk gap and are equal to zero for the present model with particle-hole symmetry (PHS). Now we paste the two ends again with some kind of 'glue potential'; it is possible that these end states can be trapped or mixed around the connecting point and evolve into in-gap bound states. This process is sketched in figures 1(a) and (b). The shapes of the possibility density of the wavefunction of our solutions for a δ-potential support this intuitive picture for the formation of the in-gap bound states. An impurity located at x = 0, unlike the open boundary, allows tunneling between the two ends of the chain, and will affect the behavior of the wavefunction near the point of x = 0. For the 1D δ-potential in a 1D TI, as proved in section 3.1.2, the bound states induced by it are always there regardless of its strength. For comparison, a pair of bound states induced by δ-potential for m B < 0 is also possible when 0

2D
Similar to the 1D case, by changing the basis order from { 0 1,2,3,4 } to { 0 1,4,2,3 }, the modified Dirac model in 2D ( p z = 0) can be reduced to two uncoupled blocks, where h − is the time-reversal partner of h + . Here, we adopt polar coordinates (x= r cosθ , y= r sinθ) to solve the equation. Since the impurity potential is isotropic, the z-component of the angular momentumĵ z = −ih∂ θ ± (h/2)σ z still commutes with h ± , as in the circular vacancy case [31]. Thus the eigenstates here can also be labeled by its quantum numbers m j ∈ {±1/2, ±3/2, . . .}. The trial solutions of the eigenstates have the following forms: for h + and h − blocks, respectively. By putting them into the Dirac equations h ± ± = E ± , the problem is simplified to two sets of equations for the radial part,  where From these equations, it is clear that due to the appearance of the finite impurity potential, PHS is destroyed while TRS survives since the impurity is non-magnetic. Suppose that there is a solution of h + block with eigenenergy E, angular momentum m jh and radial solution [ f 1 (r ), f 2 (r )] T . We denote this eigenstate as +,E,m j , where '+' indicates that it comes from h + block. It is easy to check that the state −,E,−m j with eigenenergy E, angular momentum −m jh and radial solution [ f 1 (r ), − f 2 (r )] † is an eigenstate of h − block and also the timereversal counterpart of +,E,m j . Thus, in the original basis order { 0 1,2,3,4 }, the whole 4 × 1 eigenwavefunction with eigenenergy E has the general form up to a normalization factor, where C ± are superposition parameters. Based on this general solution, one can calculate the charge current density j c , where ± denotes the h ± blocks, respectively. In fact in 2D, the modified Dirac model (1) is the same with the model first introduced by Bernevig, Hughes and Zhang (BHZ) [7] for HgTe/CdTe QWs, as long as we adopt the following basis correspondence: (14) then describes the spin-up (-down) subspace. In this sense, equation (19) shows that the components of the total eigenwavefunction that belong to h ± blocks form a pair of states with opposite helicities.
Note that m j does not appear in the eigenvalue equation; thus there is a 2|κ|-degeneracy for the corresponding eigenenergy E. Similar to the 2D case, among the eigenstates of the 3D modified Dirac model, PHS fails, while TRS survives. TRS resides in the m j ↔ −m j pair in the 2|κ| degenerate states with energy E. Suppose that there exists a solution with eigenenergy E, angular momentum quantum number ( j, m j ), spin-orbit quantum number κ and eigenvector [φ A j,m j g 1 (r ), φ B j,m j g 2 (r )] T . We denote this eigenstate as E, j,κ,m j . Apparently, the state with angular momentum quantum number ( j, −m j ), spin-orbit quantum number κ and eigenvector [(σ x φ A j,m j )g 1 (r ), (−σ x φ B j,m j )g 2 (r )] † is the time-reversal counterpart. This state can be denoted as E, j,κ,−m j .
From the above discussions, formal solutions of 2D and 3D in-gap bound states are obtained. However, there are no further analytical results under a general isotropic potential V (r ). Numerical calculations are then needed to investigate the existence and behavior of ingap bound states under selected potential forms and material parameters.

In-gap bound states in 2D and 3D TIs induced by a single non-magnetic Gaussian impurity
As illustrated in section 3.1.2, for the δ-potential, the integral equation in 2D and 3D will diverge, so the δ-potential is not strong enough to form bound states in 2D and 3D TIs. Therefore the δ-potential is different from the vacancy or defects as discussed in our previous work, which always induce the in-gap bound states in the topological nontrivial regime [31]. In this section, based on the formal solutions obtained in sections 3.3 and 3.4, we use isotropic Gaussian potentials to mimic the influence of a real non-magnetic impurity on 2D and 3D TIs and numerically investigate the existence and features of the bound states induced by it. For the Gaussian potential, we note that the in-gap bound states can be induced for a strong strength for a TI. The effective range of the potential is also an important index of the formation of the bound states. We show that for the realistic parameters of TIs and impurity, the in-gap bound state solutions can be found.

Methodology
When an impurity appears on a 2D film or inside a 3D bulk made of TIs, the momentumk is no longer a conserved quantity and is replaced by the real space gradient operator −i ∇. Since the impurity potential V (r ) is isotropic, the 2D and 3D equations can both be reduced to 1D radial equation sets (equations (16) and (21)), characterized by good quantum numbers of the corresponding conserved quantities. Using the standard finite difference method [39], the 1D radial equations can be solved to obtain the eigenenergies and the corresponding radial eigenwavefunctions spontaneously. By checking the behavior of the eigenwavefunctions within the bulk gap, those that satisfy '| (r )| 2 r d−1 → 0 when r → +∞' are the in-gap bound states we expected, where (r ) is the eigenwavefunction and d the dimensionality.
In this section, we will use two types of Gaussian potential to mimic the influence of a non-magnetic impurity. First is the standard Gaussian distribution V G (r ), In this form, the unit of V 0 is energy multiplied by volume element in R d -space. r 0 controls the effective range of the impurity potential. Equation (22) converges to V 0 δ(r ) as r → 0. This form can help us to check whether there are bound states as the impurity potential approaches a δ-function. For parameters from 2D Bi 2 Se 3 films or 3D Bi 2 Se 3 bulk, we use a modified Gaussian distribution U G (r ), which can describe a real non-magnetic impurity in a better way: In this form, U 0 directly describes the impurity potential height with the pure energy unit; r 0 still characterizes its effective range. For example, for a single Ag atom deposited onto a 2D Bi 2 Se 3 film, U 0 ≈ 3-5 eV and r 0 ≈ 3-4 Å.

2D
From section 3.2, for a 1D TI accompanied by a δ-impurity, m B = 1/2 is a critical point that distinguishes the number of surviving bound states as V 0 becomes small enough. In higher dimensions, this feature helps us determine the choice of parameter sets. In 2D, we choose several typical parameter sets, including model parameters and those from real TI materials. On the other hand, since in the 2D modified Dirac model there always exists a twofold degeneracy originating from TRS for h ± blocks, we can simply calculate the eigenstates for the h + block.
The numbers of independent solutions for the whole 4 × 4 Hamiltonian can be obtained by just doubling that from the h + block.
4.2.1. m B = 1. We choose model parameters mv 2 = 1, Bh 2 = 1 andhv = 1 as an example of the m B > 1/2 case. The standard Gaussian impurity potential V G (r ) is adopted. Under a specific combination of (V 0 , r 0 ), the radial equation (16) is solved for each m j . After checking all m j , the total number of in-gap bound states is obtained. By sweeping V 0 and r 0 in the parameter space, the phase diagram of how many pairs of in-gap bound states can be trapped is obtained, as shown in figure 2(a). It was found that as r 0 → 0, no in-gap bound states exist. This confirms the prediction in section 3.1.2 that in 2D a δ-potential cannot trap in-gap bound states. As an example, for V 0 = 7 and r 0 = 3, there are in total six pairs of in-gap bound states, with each pair being time-reversal counterparts and having the same energy. These six energy levels are shown in figure 2(b), as well as the quantum number m j for h + block. In figure 2(c), the probability distributions of bound states are plotted. Each plot indeed describes the distribution of both degenerate eigenstates since their module squares are the same. It is clear that the states with opposite m j (for h + block) do not have opposite energies. This comes from the failure of PHS once finite V G (r ) is added. At last, as |m j | increases, the bound state energy increases and, at the same time, the eigenwavefunction broadens.

HgTe quantum well (QW).
The second example is the HgTe QW. From the experimental data [7], for a 70 nm thick HgTe QW, mv 2 = −10 meV,hv = 364.5 meV nm and Bh 2 = −686 meV nm 2 . Thus one has m B ∼ 0.05. The standard Gaussian impurity potential V G (r ) is adopted again. The phase diagram of pair number of in-gap bound states is shown in figure 2(d).
Since the dimensionless parameter m B is quite small (very close to the topologically trivial regime), the branches with m j = ±1/2 in the phase diagram are greatly distorted with each other. This makes figure 2(d) look quite different from figure 2(a). Again as r 0 → 0, no in-gap bound states exist. Here we take V 0 = 13 eV nm 2 and r 0 = 2 nm as an instance. Our calculations show that there are totally two pairs of in-gap bound states, with m j = ±1/2 (for h + block). Their energies and wavefunction probability distributions are shown in figures 2(e) and (f), respectively. Also, the two energy levels are not symmetric about zero energy since PHS does not hold any more. The pair with m j = 1/2 is much broader than that with m j = −1/2.

Bi 2 Se 3 film.
In recent years, Bi 2 Se 3 films composed of several quintuple layers have attracted more and more interest [36,37,[40][41][42][43][44]. Impurity atoms, such as Ag, can be deposited onto these films to investigate their influence on electronic and transport features. Here we choose a two-quintuple-layer Bi 2 Se 3 as our 2D TI film. The following parameters [41] are adopted: mv 2 = 0.126 eV,hv = 2.94 eV Å and Bh 2 = 21.8 eV Å 2 . This gives m B ∼ 0.32 < 1/2. To mimic a real non-magnetic impurity, the modified Gaussian potential U G (r ) is used. Following a similar process, the (U 0 , r 0 ) phase diagram of how many pairs of in-gap bound states can be trapped by a single non-magnetic impurity on a two-quintuple-layer Bi 2 Se 3 is shown in figure 2(g). As an example, the influence of an Ag atom can be modeled by the Gaussian potential U G (r ) with U 0 = 4 eV and r 0 = 3 Å. Our calculations show that in this case, there are two pairs of in-gap bound states with m j = ±1/2 (for h + block), as shown in figure 2(h). The corresponding wavefunction probability distributions are plotted in figure 2(i). Rather than the HgTe QW, these results for Bi 2 Se 3 films should be easily checked by experiments using STM, etc.

3D
In 3D, we choose a set of parameters based on those from the anisotropic model for bulk Bi 2 Se 3 [18]. In that model, one has A 1 = 2.2 eV Å, A 2 = 4.1 eV Å (correspond tohv) and B 1 = 10 eV Å 2 , B 2 = 56.6 eV Å 2 (correspond to Bh 2 ). We neglect the anisotropy and simply choose their averages as our parametershv and Bh 2 , respectively. In this case, m B ∼ 1 > 1/2. An interior non-magnetic impurity is modeled by the modified Gaussian potential U G (r ). Under a specific combination of (U 0 , r 0 ), first for a certain κ, equation (21) is solved to see whether there exist in-gap bound states. Also one should not forget the additional 2|κ| degeneracy for a specific κ. After checking all possible κ, we obtain the total number of in-gap bound states for this (U 0 , r 0 ) combination. In this way, we obtain the phase diagram of the number of in-gap bound states in (U 0 , r 0 ) parameter space, as shown in figure 3(a). In the simplest way, a cluster of Ag atoms inside a bulk Bi 2 Se 3 can be modeled by U 0 = 3 eV and r 0 = 6.5 Å [45], as indicated by the small red circle in the phase diagram, i.e. figure 3(a). Our calculations show that in this case, there are four bound-state solutions of the radial equation (belongs to κ = ±1, ±2, respectively), as well as four corresponding bound state energy levels. This is shown in figure 3(b). After considering the 2|κ| degeneracy, there are totally 12 in-gap bound states. Each one can be uniquely labeled by the three good quantum numbers ( j, κ, m j ). The module square isosurfaces of | ( j,κ,m j ) | 2 ≡ 10 −3 for the 12 bound states are plotted in figure 3(c). To be more clear, each bound state has the same color in energy level and module square isosurface. Again, PHS is broken while TRS still holds.
From equation (20), generally the angular part of a 3D in-gap bound state is a mixture of j ± 1 2 components. However, the module square of the whole wavefunction has simpler symmetry, as shown in figure 3(c). We found that the bound states for κ = ±1 (correspond to j = 1/2) have s-symmetry, while those for κ = ±2 (correspond to j = 3/2) have p-symmetry. This can be understood by directly putting the orthonormalized Laplace spherical harmonics Y m l into the wavefunctions. After some simple algebra, finally one has |ψ 3 2 ,±2,± 3 2 | 2 = 3 sin 2 θ 8π · |η 1 (r )| 2 + |η 2 (r )| 2 , where g(ζ, η) 1,2 (r ) are the solutions of equation (21) for the corresponding ( j, κ, m j ) state. Then the symmetry features of the module square of the whole wavefunction are apparent.
To end this section, we point out that the main purpose of this work is to provide a simple but inspiring way of understanding how non-magnetic impurities affect the electronic energy structure of TIs in various dimensions. The Ag impurity represents a typical example, from which the height and width of the Gaussian potential can be estimated. For Ca [46] and Pb [47] impurities, we expect the same physics. As we know, the Ca atom is lighter than the Ag atom, but Pb is heavier than the Ag atom. If an Ag impurity could induce a bound state, Pb could also induce a bound state. However, Ca may not. As long as one knows the model parameters for Ca or Pb, similar numerical calculations can be performed and the picture we proposed here should be further verified. The ring formed by these bound states (red) circulating around an impurity atom (blue) deposited on a thin TI film can be detected by STM when the Fermi surface is aligned with the bound states in the gap (b). (c) The resonance peak can be detected in the differential conductance spectrum (dI /dV ) as a function of the bias voltage (V) when the STM tip is pinned above the bound-state ring in real space. (d) The resonance peak splits for magnetic impurity.

Detecting in-gap bound states by using a scanning tunneling microscope
To detect the in-gap bound states obtained in the previous sections, we propose to observe by STM the vicinity of a positive-valence atom deposited on a gapped ultrathin film of 3D TI ( figure 4(a)), which is a direct realization of the 2D modified Dirac model with impurity barrier potential. STM can associate the local density of states (LDOS) of sample surfaces with the bias voltage dependence of tunneling current [48], allowing the spatial wavefunction at a certain energy and the energy spectrum at a certain spatial point to be probed.
When fixing the bias voltage to align the Fermi energy with the bound state energy in the gap ( figure 4(b)), the scanning of the sample surface should reveal one or more high LDOS rings around the impurity atom associated with the bound states. Note that the perfect ring predicted by the isotropic model should be reshaped with respect to the actual lattice symmetry in real materials, such as into hexagons. Alternatively, when the STM tip is pinned at the boundstate ring, the differential conductance (dI /dV ) as a function of the bias voltage (V ) should demonstrate resonance peaks at the bound state energies in the gap ( figure 4(c)). It should be noted that each energy peak should be twofold spin-degenerate for nonmagnetic impurity because of TRS. Therefore, one expects the energy peak to split if a magnetic impurity is adsorbed instead ( figure 4(d)).
By far, the STM studies of the impurity on the surface of 3D topological insulators reported no signature of the topologically protected bound states. Instead, the LDOS map is dominated by profound standing wave patterns, which originate from the quantum interference of conduction electrons from surface bands at the Fermi surface [45,49]. The reason for the absence of the bound states is because the gapless Dirac cone of the surface states of a 3D TI does not support the presence of the bound states. In addition, we emphasized that the ring pattern from the bound states should be distinguished from those from the quantum interference because the latter requires that the Fermi surface crosses a band and has a length scale of the Fermi wavelength [45].

Summary
In this work, the in-gap bound states induced by non-magnetic impurities in various dimensions are systematically investigated based on a modified Dirac model. Distinct from the vacancy case, the impurity-induced in-gap bound states can have nonzero components at which the impurity locates. However, generally they are not topologically protected. The format of the impurity potential and the dimensionality of the TI are both important for their formation. δ-potential has been proven to be able to trap in-gap bound states in 1D, while fails in 2D and 3D. For general forms of the impurity potential, formal solutions were obtained and numerical calculations were performed to explore the existence and behavior of in-gap bound states. In this paper, we focused on isotropic Gaussian potentials and obtained the phase diagrams of how many bound states can be trapped in 2D films or 3D bulks for parameter sets from a series of topological nontrivial materials. For some real experimental configurations, such as Ag atoms deposited on a Bi 2 Se 3 film or inside a bulk Bi 2 Se 3 , in-gap bound states were found to be possible. At last, we proposed a possible method of detecting these in-gap bound states induced by impurities on 2D thin film or the surface of 3D bulk by STM measurements. These results should be of interest to people working in this field, and should attract further attention to the bound state physics of TIs.