Persistence of symmetry-protected Dirac points at the surface of the topological crystalline insulator SnTe upon impurity doping

Abstract

We investigate the effect of a non-magnetic donor impurity located at the surface of the SnTe topological crystalline insulator. In particular, the changes on the surface states due to a Sb impurity atom are analyzed by means of ab initio simulations of pristine and impurity-doped SnTe. Both semi-infinite and slab geometries are considered within the first-principles approach. Furthermore, minimal and Green's function continuum models are proposed with the same goal. We find that the Dirac cones are shifted down in energy upon doping; this shift strongly depends on the position of the impurity with respect to the surface. In addition, we observe that the width of the impurity band presents an even–odd behavior by varying the position of the impurity. This behavior is related to the position of the nodes of the wave function with respect to the surface, and hence it is a manifestation of confinement effects. We compare slab and semi-infinite geometries within the ab initio approach, demonstrating that the surface states remain gapless and their spin textures are unaltered in the doped semi-infinite system. In the slab geometry, a gap opens due to hybridization of the states localized at opposite surfaces. Finally, by means of a continuum model, we extrapolate our results to arbitrary positions of the impurity, clearly showing a non-monotonic behavior of the Dirac cone.

---

1. Introduction

The discovery of topological insulators (TIs) in the past decade has fueled great interest in Dirac matter. These materials are insulating in the bulk, but host protected electronic states in their boundaries; namely, surface or edge modes in the three-dimensional or two-dimensional case, respectively.^1,2 Such gapless boundary states are protected by time-reversal-symmetry; they have a Dirac spectrum and show spin-momentum locking, a promising property for their application in spintronic devices and quantum information processing.³ TIs are successfully characterized by the [Doublestruck Z]₂ topological invariant, which is related to the number of gapless Kramers pairs in their boundaries. In fact, they show an odd number of Dirac cones at their surfaces. The discovery of topological insulators has led to the recognition that symmetry-protected topological phases of matter are more abundant and ubiquitous than expected.^4,5

One of the families within the broad class of symmetry-protected matter are topological crystalline insulators (TCIs). Unlike [Doublestruck Z]₂ insulators, TCIs are protected by crystal symmetries and present an even number of Dirac cones in their surfaces, so that their [Doublestruck Z]₂ index is zero. These symmetries include rotations, reflections, or even glide planes. The first TCIs discovered were IV–VI semiconductors, such as SnTe and Pb_1−xSn_xTe ternary alloys, which are protected by mirror symmetries.^6–9 When a mirror symmetry is responsible for the topological protection, the relevant topological invariant is the mirror Chern number n_M. There are two types of surface states in TCIs, depending on the location of the Dirac points, with different properties. Type-I surface states have their Dirac points at time-reversal invariant momenta, whereas in type-II TCIs the Dirac cones are displaced from such time-reversal invariant k-points.¹⁰

SnTe has a rocksalt crystalline structure with small bandgaps at the L points of the Brillouin zone (BZ), where the conduction and valence bands are inverted.^6,11 The mirror plane that protects these topological states is (110) [see Fig. 1(c)], so that an even number of robust Dirac cones appear at the (001), (110) and (111) crystal surfaces, which are symmetric with respect to the (110) mirror reflection. For the (111) termination, surface states are of type-I, whereas (001) and (110) surface states are identified as type-II. For the [001] surface termination, which is of type II, two equivalent L points of the BZ are projected in the ΓL₁L₂ plane to the same X̄ point, as depicted in Fig. 1(b). Hence, two parent Dirac cones coexist at X̄, and their interaction gives rise to two pairs of massless Dirac fermions deviated from X̄, while the high-symmetry point is found to be gapped¹⁰ [see Fig. 1(d)]. SnTe undergoes a crystalline phase transition at low temperatures (below 100 K) to a rhombohedral lattice. Such distortion has a small effect in the bulk bands, but destroys the surface Dirac cones, as any other perturbation that breaks the mirror (110) symmetry would do.¹² Interestingly, within this rhombohedral phase, SnTe has been recently shown to be a higher-order TI, showing hinge states along specific directions.¹³

Fig. 1 (a) SnTe 8-atomic-layer slab unit cell. Sb impurities are highlighted in blue: in the Sb1 case, the whole atoms are colored; in the Sb2, Sb3 and Sb4 cases, the impurity atoms are circled. Higher slabs can be constructed adding new atomic layers in the dotted region. A 16 Å-wide vacuum space was considered in the z-direction in order to avoid interaction among the slabs during the first-principles calculations. (b) Face-centered cubic Brillouin zone and projected (001) plane. (c) Top view of the c(4 × 4) (light blue) and 1 × 1 (dark blue) unit cells. The mirror symmetry plane is represented by a dotted line. (d) SnTe (001) surface band structure of the semi-infinite system.

On the other hand, PbTe, which also has a rocksalt structure, is not a TCI and its surface states are gapped. Notwithstanding, a topological phase transition can be induced either by pressure or alloying. In particular, as anticipated above, Pb_1−xSn_xTe is a topologically non-trivial material (n_M = −2) for x ⪆ 1/3.^7,9 Actually, a band inversion as a function of the alloy composition was reported long time ago.¹⁴

Thin films of these TCIs also show nontrivial topological behavior. However, hybridization between the surface states opens a gap that decreases for increasing layer thickness.^15,16 In fact, a non-monotonic damped oscillation has been found, that it is ultimately related to the alternation between nonsymmorphic and symmorphic symmetry with the number of atomic layers.^15,17 Specifically, if a [001]-oriented SnTe slab has an even number of layers, it possesses nonsymmorphic symmetry. This reduces the hybridization between surface states, whereas for an odd number of layers the symmetry is symmorphic, and more sizeable band gaps appear.

Additionally, symmetry-protected surface states in TCIs can be tuned by introducing impurities. The effects of doping on the electronic structure of bulk SnTe have been addressed in ref. 18 and 19. As for the (001) surface, In doping,²⁰ hydrogen adsorption²¹ and diluted impurities²² have been considered. Doped SnTe, specially by In impurities, is an ongoing subject of discussion regarding topological superconductivity.²³

Impurities in TCI thin films may open a gap or modify the surface states, allowing for the control of topological phase transitions in these materials and constituting an additional tool to modulate their properties.^24,25 Indeed, tuning the electronic properties by doping or by applied external fields opens a way for the application of these materials in nanoscale devices.¹⁶ Such modification by doping may be also relevant for thermoelectric applications, another arena in which topological materials have been proven to be of interest.^18,19

In this work, the effect of non-magnetic impurities in the topological crystalline insulator SnTe is studied. Specifically, we consider Sn substitution by a Sb atom, which acts as a donor impurity. In fact, Sb-doped SnTe has been found to be a superconductor.²⁶ Here, we focus on the changes in the TCI dispersion relations induced by the presence of the Sb impurity. As mentioned above, doping may alter or modify the surface states by hybridization or just by opening a gap, therefore changing the electronic properties of the system. This is the main issue we address herein. We explore both semi-infinite and slab geometries, with a twofold goal. Firstly, ultrathin layers of TCIs are currently grown and incorporated in quasi-2D materials and devices as building blocks or integrating parts, so the thin slab geometry has an interest in itself. Secondly, in most first-principles calculations, a finite slab geometry is employed to model surfaces, so it is equally important to evaluate the differences between finite-size and semi-infinite calculations aimed at modeling the same system. To this end, DFT-based ab initio simulations are used to characterize SnTe and the effect of the non-magnetic donor impurity near the crystal surface. We analyze the hybridization between surface and impurity states as well as the effect of the impurity position on the electronic structure and crystal symmetries. Furthermore, both perturbation theory and continuum models are proposed. The former is used to predict the behavior of the Dirac cone upon doping, while the latter describes the impurity states.

Our main results are summarized as follows: (i) the Dirac cones appearing at the TCI surfaces are shifted down upon the n-doping due to Sb substitutional impurities. (ii) The shift of the Dirac cones is strongly dependent on the position of the impurity with respect to the surface. (iii) The width of the impurity band presents an even–odd behavior with the layer position of the dopant, both in the semi-infinite and the slab geometries; we relate it to the location of the nodes of the wave function with respect to the surface. (iv) In the slab system, a gap opens in the Dirac points due to coupling of the states at opposite surfaces, even though they are almost 50 Å apart. We verify that the surface states remain gapless in the semi-infinite geometry, with unaltered spin textures. (v) Finally, we derive continuum models in order to extrapolate the behavior of this doped system for larger sizes, unreachable with more accurate, ab initio simulations.

Our work provides a complete picture of the interplay between impurity and surface states in TCIs, highlighting the differences between finite-sized slabs and single surfaces, and additionally including simpler continuum models which yield a better interpretation of the numerical first-principles results.

2. First-principles approach

2.1. Geometry and methods

Our calculations were performed using the pseudopotential DFT SIESTA code^27,28 within the generalized gradient approximation under the PBE parametrization (GGA-PBE) for the exchange and correlation functional.²⁹ Spin–orbit coupling was considered within the fully-relativistic pseudopotential formalism³⁰ for all calculations unless otherwise stated. Accordingly, a relativistic pseudopotential was employed, following the Troullier–Martins method.³¹

A double-ζ singly-polarized basis set was chosen for all the atoms involved, with a radial extension of the strictly localized orbitals determined from a confinement energy (energy shift) of 100 meV. The resolution of the real space grids was set to around 0.06 ų (equivalent to a mesh cutoff of 600 Ry). For the integration over the BZ we employed k-supercells of size around (17 × 17) relative to the SnTe-(1 × 1) bulk unit cell, with a temperature in the Fermi–Dirac occupation function of 100 meV.

A lattice parameter of 6.39 Å was used following ref. 6 and 18, which yields a satisfactory agreement to experimental results and a successful magnitude of the band gap.¹⁴ This allows us to model rocksalt room-temperature SnTe, avoiding the aforementioned rhombohedral phase transition.⁹ In order to tackle the robustness of the surface states against film thickness, several slab calculations were performed using different number of layers (not shown here). We found an optimal compromise between accuracy and our computational capability for a nonsymmorphic 16 atomic-layer SnTe slab. Therefore, all the slab calculations presented in this work have been performed with this thickness.

Substitutional donor impurities have been considered by replacing a Sn atom by Sb within a c(4 × 4) supercell—that is, in a defective layer one out of eight Sn atoms is substituted by a Sb atom. We have considered up to four different locations of the Sb: at the top surface layer (referred to as Sb1 hereafter), the second layer (Sb2), the third layer (Sb3) and the fourth layer (Sb4). Furthermore, in order to avoid net dipoles in the slab, a second Sb substitutional impurity was always placed at an equivalent position, but with respect to the bottom layer [see Fig. 1(a)]. In this way, both impurities are located at the same distance from the corresponding nearest surface. Thus, the simulated impurity doping in the slabs corresponds to a small value of x = 0.016.

Atomic relaxation was performed for each structure. For the Sb1 slab, the two top and bottom layers were relaxed while for the Sb2, Sb3 and Sb4 slabs, the 3, 4 and 5 layers counted from the top and bottom surfaces were included in the relaxations until forces on the atoms were below 0.02 eV Å⁻¹. The top and bottom layers of the final equilibrium structure are the most affected by the relaxation. To be specific, some atoms of the Sb1 and Sb3 slabs are displaced up to 0.17 Å compared to 0.14 Å for the Sb2 and Sb4 slabs.

Although minimal, the hybridization of opposite surface states still has an effect on the band structure. Hence, the electronic and spin structure of the Sb-doped SnTe(001) surface will be mainly presented as calculated under a semi-infinite geometry via projected density of states and magnetization maps, PDOS(k,E) and M_α(k,E) respectively, with α = x,y,z indicating the direction.^32,33 To this end, we have employed surface Green's function matching techniques in order to couple the Hamiltonian matrix elements of the top half of the slab with those of a bulk SnTe calculation. Self-consistency was preserved in the matching process, as reflected by the fact that maximum deviations were always smaller than 10 meV between the on-site energies of the Sn and Te central layers of the slab and those in the SnTe bulk. PDOS(k,E) and M_α(k,E) maps were computed with a resolution of 0.003 Å⁻¹ in reciprocal space and 5 meV in energy. Accordingly, the self-energy that determines the spectral width of any surface states (imaginary part of the energy entering the Green's function) was set to 5 meV. For zoom-ins of the Dirac points we computed additional maps with an increased resolution (0.0016 Å⁻¹ and 3 meV) around the [capital Gamma, Greek, macron]-point and over a small energy window. In the semi-infinite calculations, a single impurity is used. For topologically-protected surface states, we expect them to remain unaltered, since the placement of a single impurity does not change the (110) mirror symmetry, as outlined in Fig. 1(c). On the contrary, the 16 atomic layer slabs do break this symmetry when taking into account two impurities, as explained before. Since two-dimensional (few-layer) topological crystalline insulators are also the subject of great attention, we consider this comparison to be relevant.

2.2. Ab initio results

Top panel of Fig. 2 displays the effect of the substitutional Sb on the surface electronic structure for the four locations of the impurity as it is moved from the surface layer towards the bulk, namely Sb1 to Sb4. The band structures have been computed along the [capital Gamma, Greek, macron]–X̄ direction for 16-layer slabs (that is, under an ultrathin film geometry). Note that, due to band folding, the X̄ point in the (1 × 1) BZ back folds into the [capital Gamma, Greek, macron] point of the supercell's BZ. For the sake of comparison, we also include in the left column the k-folded SnTe pristine surface computed under the same c(4 × 4) supercell as that used for the defective structures. There are several features immediately apparent in the plots. First of all, in all defective surfaces the Sb develops a split impurity band that disperses across the gap. Interestingly, the band is narrower when the impurity is located at an odd-numbered layer (particularly the Sb1 case but also Sb3), due to a higher localization of the states at the surface. On the contrary, when the impurity is located at an even-numbered layer (Sb2 and Sb4 cases), the band becomes more dispersive and runs closer to the edge of the bulk conduction band. Second, the presence of the impurity induces a clear n-type doping as the Fermi level E_F shifts towards the minimum of the conduction band in all defective structures, while it is pinned just above the valence band (VB) maximum for the pristine surface. This doping effect will be further addressed in section 3 using a minimal continuum model. Despite the large thickness of the slabs employed, all Dirac points present a gap due to the long decay length of the topological surface states, so that there is still a sizeable interaction between the top and bottom surfaces. The size of this gap is minimal for the odd cases (Sb1 and Sb3), being around 5 meV, while it attains substantially larger values for the even cases. Note that it is around 50 meV for Sb2 and 30 meV for Sb4. For the pristine slab, for which the gap is due to finite-size effects, it is 55 meV. Moreover, the breaking of the (110) mirror symmetry in the two-impurity geometry used in the slabs affects the symmetry-protected surface states and contributes to the gap opening, as depicted in panels (b)–(e). Although such even–odd behavior in the Dirac cone gap value with the impurity position is reminiscent of the nonsymmorphic–symmorphic effect with respect to the number of layers previously reported,^15,17 we note that in our case its origin should be different, because all slabs are composed of the same number of layers.

Fig. 2 Electronic structure of pristine SnTe (left column) and Sb-doped SnTe, with the donor impurity located in the first, second, third and fourth atomic layer. Top panel shows the band structure of the 16-layer slabs, as well as the bulk band projections of pristine SnTe represented by the gray shaded area. Bottom panel displays the projected density of states PDOS(k,E) on the first atomic layers for the semi-infinite system.

The bottom panel of Fig. 2 displays the equivalent electronic structures to the top panel but in the form of PDOS(k,E) maps, since they have been computed under a semi-infinite geometry (see previous section). The Sb-doped surfaces, overall, show a similar behavior to that found for the 16-layer slabs, namely, the appearance of Sb impurity bands across the gap showing an even–odd effect in their dispersion, a clear n-type doping and the presence of Dirac cones slightly away from [capital Gamma, Greek, macron]. However, these Dirac cones remarkably persist in all defective cases of the semi-infinite structure. The splitting of the impurity bands induced by the spin–orbit coupling remains small until it approaches the [capital Gamma, Greek, macron] point, where the two branches start to deviate from each other by more than 100 meV.

Focusing now on the electronic structure around the [capital Gamma, Greek, macron] point, the most relevant feature is the survival of a gapless Dirac point in all structures, which is consistent with the fact that the (110) mirror plane still holds in the doped structures and, hence, the Dirac points remain symmetry-protected. However, there are qualitative differences with respect to the pristine surface which are best seen in Fig. 3, where zoom-ins with an enhanced resolution are presented. First of all, and in accordance with the n-type doping mentioned above, the Dirac points consistently shift by around 150–200 meV towards lower energies. Precisely at the [capital Gamma, Greek, macron] point the interpretation of the electronic structure becomes more complex due to the profusion of extra bands arising from both the Sb impurity and the band folding.

Fig. 3 Zoomed (a) PDOS(k,E) as well as (b) M_x and (c) M_y maps of the semi-infinite structure. In panel (a), the Sb impurity contribution is highlighted in blue. PDOS(k_x,k_y,E) maps at an energy 20 meV (d) below and (e) above the Dirac point. The spin textures are depicted with arrows that outline the spin direction. See main text for an explanation of labels a, b and f in the upper panel.

All relevant bands within the gap are better resolved in Fig. 3(a), where we present enlarged views of the PDOS(k,E) maps. The projections have been taken over the surface SnTe layer, superimposing the Sb contributions in blue.

For the clean surface case, the branches forming the Dirac cone can be clearly identified, as well as the Sb-split impurity bands, denoted as a and b. There is a further band inside the gap, labeled by f, corresponding to the back-folded topological surface state along X̄–M̄ [see Fig. 1(d)]. Since the supercell size is much smaller than the effective Bohr radius, impurities cannot be regarded as truly isolated and the impurity bands become dispersive. In addition, further broadening of the bands arises from the hybridization of impurity states with Bloch states of the conduction band of SnTe. In fact, the impurity states are mainly composed of Sb and Sn contributions, as depicted in the Appendix A for both slab and bulk band structures. Band inversion between Sn and Te orbitals also persists in the doped material, as discussed and shown in Appendix A.

In the lower panels (b) and (c) of Fig. 3, magnetization maps along the x- and y-axis, M_x(k,E) and M_y(k,E), respectively, are shown. We omit the out-of-plane component M_z(k,E) since it is always negligible except for Sb1. As it could be expected, the in-plane spin orientations of the surface band remain unaltered with respect to the pristine case. Indeed, this is best seen in the PDOS(k_x,k_y) maps shown in the lower panels, computed at energies 20 meV below [panel (d)] and above [panel (e)] the location of the Dirac point.

The spin texture for all structures is shown in Fig. 4 in the form of M_x/y(k,E) maps. Note that pristine SnTe only shows in-plane spin polarization due to time-reversal and mirror symmetries.³⁴ The effect of folding on the spin texture of pristine SnTe is explicitly detailed in the Appendix A, showing additional M_x/y(k_x,k_y) maps near the Dirac point for all structures. A close inspection of the spin textures reveals that the spin orientation of the topological surface states always remains unchanged with respect to the pristine surface. More surprising is the behavior of the Sb bands: close to X̄ they behave as standard Rashba-split states, with a reduced separation in energy between a and b of only 10–20 meV while they present opposite magnetizations. In fact, this is the only region where out-of-plane magnetization, M_z(k,E), is not negligible. However, as they disperse towards [capital Gamma, Greek, macron], band a keeps its spin orientation but b undergoes an abrupt spin inversion and ends up with the same orientation as a. Note that the spin texture of the impurity states is not as robust as that of the Dirac cone, since it varies along the high-symmetry path and is dependent on the choice of atom projections.

Fig. 4 Magnetization maps (a) M_x(k,E) and (b) M_y(k,E) for the semi-infinite pristine and Sb-doped structure, projected on the first atomic layers. White and black shades account for positive and negative values of the magnetization M_α(k,E), respectively.

Given such ultrathin layers, the energy shift of the Dirac point can be related to quantum size effects. In order to clarify this point, we study the spatial distribution of the wave function for the impurity band, depicted in Fig. 5. A series of nodes appear right after the impurity and approximately each two atomic layers in the z-direction. When the impurity is located at an odd atomic layer, there is an even number of layers above it, so the wave function amplitude is almost zero at the surface layer. In contrast, when the impurity is placed in an even layer, there is an odd number of layers above it and the wave function is nonzero at the surface, therefore favoring hybridization with the surface states and a larger band gap. Since the derivative of the wave function is related to the kinetic energy operator, it is also expected for these impurity bands to move upwards, as shown in both slab and semi-infinite calculations in Fig. 2. Thus, we come to the conclusion that this energy shift can be also seen as a size or quantum confinement effect.

Fig. 5 Wave functions for one of the impurity bands at [capital Gamma, Greek, macron]. Yellow and blue colors represent positive and negative values, respectively. The colored rectangles next to each panel, which are two atomic layers wide, outline the alternation of regions with positive and negative values of the wave function in the z-direction. The layer where the impurity is located is highlighted in dark blue.

3. Continuum model

In order to have an insight into the effect of the impurity upon the Dirac cones and the behavior of the impurity band when the position of the donor Sb atom is varied, we resort to a simple continuum model of a TI. The SnTe family has two non-equivalent valleys, L₁ and L₂, in the Brillouin zone. However, a single valley approach successfully explains the existence of robust surface states on crystal surfaces such as (001), (110) or (111).⁶ Hence, in this work we focus on states close to one of the L points of the Brillouin zone and neglect other valleys since surface states are stable against gap opening by valley mixing. It consists of an isotropic two-band spin-full continuum Hamiltonian in a particle-hole symmetric scenario. In fact, this model has been already employed for TIs.^35,36 The Hamiltonian in momentum space is given by (we set ħ = 1 hereafter)where ε is a constant energy term, 1₄ is the 4 × 4 identity matrix, v is the Fermi velocity, m is the mass parameter related to the gap and the constant B takes into account the quadratic terms in momentum. In the former expression, α is the vector of Dirac matrices and the momentum-vector defined as α = (α_x,α_y,α_z). We choose the Dirac matrices expressed in terms of the Pauli matrices in the following basis: α_i ≡ σ_x⊗σ_i, β ≡ σ_z⊗σ₀. This simple model has the minimum elements needed to show topologically protected surface states if mB > 0 and the translational symmetry is broken by making a finite-size system.³⁷

Modifications of this model, including anisotropy and on-site momentum dependent terms, have been successfully proposed to study the Bi₂Se₃ family of materials.^36,38 In fact, the model has also been applied to describe the (001) surface states of SnTe including just the leading terms linear in momentum.^{6,10,15,17,34,39} Since our first goal within the continuum approach is to analyze the effect of an impurity in close proximity of an isolated surface, quadratic terms are indispensable to have a nontrivial topology in a vacuum–SnTe interface.³⁵ We expect to find the same shifting behavior as predicted by Araújo et al.¹⁷ and our DFT calculations. Thus, we perform first-order perturbation theory (PT) on the surface states obtained within this formalism with the aim of elucidating the effect of the impurity on the gapless surface states. Secondly, in order to describe the energy dependence of the impurity states with respect to the impurity position, we solve a simplified version of the Hamiltonian presented in eqn (1), summing up all terms in the perturbation series using a Green's function method. These two approaches allow us to investigate complementary aspects of the problem.

3.1. Perturbation theory approach

With the aim of comparing PT results with DFT calculations, we consider the case of a semi-infinite system which extends from z > 0, being infinite in the other perpendicular directions. The surface states are derived in Appendix B imposing Dirichlet boundary conditions at the plane z = 0. In this way we obtain two Kramers partners that only differ in the spinor part Φ⁰_±

ψ_± = A_s(e^−λ₁z − e^−λ₂z)exp[i(k_xx + k_yy)]Φ⁰_±. (2)

In the former expression, the inverses of the decay lengths are given by λ₁ and λ₂, which are functions of the in-plane momenta and A_s is a normalization factor (see Appendix B). The dispersion of the states is linear in the in-plane momenta, yielding the Dirac cones of the surface states inside the bulk gap E_± = ε ± vk_∥sign(B), where .

For the purpose of obtaining the corrections on the energy, we perform first-order PT, modeling the impurity with a Coulomb-like interaction

where z₀ is the depth where the impurity is located, and and are the effective Rydberg energy and Bohr radius in SnTe.^40,41 The former constants provide the correct units to the expression and scale the interaction due to the medium, as detailed in Appendix B.

The first-order correction in PT is given by

where ψ_± are the surface states as defined in eqn (2) for the impurity-free case. Due to the trivial structure of the correction in the spin space, both states, ψ₊ and ψ₋, show the same correction. Therefore, the cones are just shifted by the same amount and no gap opens.

Fig. 6 shows the displacement of the Dirac point as a function of the impurity depth z₀. Similarly to Fig. 2, a shift of the Dirac point as the impurity is placed farther from the surface is observed, in opposition to hole doping.²⁰ While DFT calculations are difficult to address due to their computational cost when the impurity is far from the surface, this minimal model allows us to predict a return to the pristine case for large z₀. Despite the aforementioned limitations, the model describes qualitatively well the phenomenon and it can be used as a simple approach to the effect of the impurity on the surface bands.

Fig. 6 Displacement of the Dirac point ΔE⁽¹⁾ as a function of the impurity depth z₀, obtained from PT calculations.

3.2. Green's function approach

In the previous subsection, the interaction of the electron with the impurity was assumed weak so that the first-order perturbation theory was a reasonable approximation. We can go beyond this limit, summing up all terms of the perturbation series by means of the Green's function approach. To this end, we start with eqn (1) and neglect the quadratic in momentum term (B = 0). The linear approximation (Dirac-like equation) is unable to explain the existence of surface states in a TI unless a normal semiconductor is attached, thus forming a heterojunction. In this way, there exists a band inversion along the direction normal to the surface and interface states can appear (see, e.g., ref. 35 and references therein).

The envelope functions of the electron states near the band extrema are then determined from the following Dirac-like Hamiltonian

where we assume aligned and same-sized gaps for simplicity. The origin of energy is set at the center of the gaps (ε = 0) and the mass term is given as m(z) = msign(z). The Hamiltonian (5) acts upon the envelope function χ(r), which is a four-component vector composed of the two-component spinors χ₊(r) and χ₋(r). It decays as exp(−|z|/d) along the direction normal to the surface, with d = v/m, and the interface dispersion relation is a single Dirac cone E(k_‖) = ±v|k_‖|.

The presence of the impurity, located at r₀ = (0,0,z₀) without loss of generality, breaks the translational symmetry in the surface and the in-plane momentum is no longer conserved. Expressing distance in units of d and energy in units of m, the Dirac-like equation for the envelope function can be cast in the form

The actual interaction potential V(r − r₀) will be replaced by a non-local separable pseudo-potential of the form

where λ is a coupling constant and ω(r − r₀) is a real function referred to as shape function hereafter. It is worth mentioning that this replacement is exact, in the sense that it is always possible to find a non-local separable pseudo-potential (or a sum of them) able to reproduce any set of given electronic states and, consequently, there is no theoretical limitation to the numerical accuracy with which physical results can be obtained.⁴² In spite of its seemingly more complicated form, eqn (7) is amenable to analytical solution for any arbitrary shape function.^43–48

The Green's function for the unperturbed problem (λ = 0) and the Green's function for the total Hamiltonian in terms of the transition operator are derived in Appendix B. Poles in the complex plane of the transition matrix, E = E_r − iΓ, where E_r and Γ are real magnitudes, yield the energies of the impurity states. Resonances within the gap are given by the condition |E_r| < 1. According to the model presented, an impurity embedded in a bulk semiconductor can support a truly bound state. Lengthy but straightforward calculations lead to the following energy for the bound state

in the limit ω(r − r₀) → δ(r − r₀) and k_c is a cut-off momentum that regularizes the δ function. Since only states within the gap correspond to bound states in this model (|E_b| < 1), there exists a minimum value of the product |λ|k_c to induce bound states. If |λ|k_c < 2π the impurity cannot bind electrons (λ < 0) or holes (λ > 0) in bulk semiconductors.

Fig. 7 shows the real and imaginary parts of the poles of the transition matrix as a function of the position of the impurity close to the heterojunction, when k_c = 100 and λ = −0.15. The real part E_r approaches the value in the bulk material E_b when increasing z₀, as expected. These impurity states move to higher energies as the distance from the surface increases, as in ref. 49. The alternating behavior found in DFT calculations, with different energy shifts of the impurity bands for even (odd) slabs, is not reflected in this model since simpler, exponentially-decaying wave functions are used here.

Fig. 7 Real and imaginary parts of the complex poles E = E_r − iΓ of the transition matrix as a function of the impurity position. The parameters of the interaction potential are λ = −0.15 and k_c = 100. Dashed line indicates the energy of the impurity bound state in the bulk semiconductor [see eqn (8)].

4. Conclusions

Doping a SnTe slab with a non-magnetic Sb impurity breaks mirror symmetry, which affects the well-known topologically-protected surface states of pristine SnTe. Additionally, impurity doping has an effect on the band structure: it creates impurity bands which hybridize with the aforementioned surface states. We have studied this doped system by resorting to first-principles calculations and continuum models. DFT calculations have been performed to study a finite structure, i.e., a slab with two surfaces, as well as a semi-infinite geometry with only one surface. We have analyzed the role of crystal symmetries for each case. Surface states are unaltered when a single substitutional impurity is used in a semi-infinite system, yielding a band structure with Dirac cones as well as spin-polarized impurity states. In contrast, a slab geometry even with 16 atomic layers presents coupling of states at opposite surfaces, which produces a small gap in the Dirac cones. Our calculations allow us to isolate the contribution of either the impurity or the first atomic layers of the surface to the electronic structure of the system, which is specially useful since we consider different impurity positions. Impurity states show an interesting alternating behavior depending on the location of the impurity inside the material, which can be regarded as an even–odd effect.

A continuum Hamiltonian is also used to describe the effect of doping on the Dirac cone, which shifts downwards. This result coincides qualitatively with our DFT calculations and allows us to predict the behavior of the system when the impurity is placed deeper in the slab. As for the hybridization between the topologically-protected surface states and the impurity states, a simplified version of the aforementioned Hamiltonian is used, showing a diminishing resonance as the impurity moves inwards.

The n-doped character of SnTe:Sb could make difficult to separate the contribution to the electronic properties of surfaces and bulk states. However, angle-resolved photoemission spectroscopy (ARPES) enables revealing the signatures of surface states (i.e. Dirac cone) even if the Fermi level lies well above the Dirac point, as was demonstrated by Chen et al. in the Bi₂Te₃ samples.⁵⁰ We are then confident that our predictions could be observed in ARPES experiments.

In conclusion, SnTe, the first discovered TCI, is an excellent platform to explore the interplay of symmetry and topology in this class of Dirac matter. A simple perturbation such as a substitutional impurity is shown to give rise to non-trivial spin textures. Besides its fundamental interest, tuning the Dirac cones of topological insulators can be of interest for transport and spintronic applications.

Conflicts of interest

There are no conflicts to declare.

Appendix A

Additional DFT results

A more detailed analysis of the orbital character and spin texture of pristine and impurity-doped SnTe is illustrated by the following figures. Fig. 8(a) displays the band structure along with the projected contributions of all chemical species (green for Sn, blue for Te and orange for Sb) for the slab structure. The band inversion at the [capital Gamma, Greek, macron] point is observed for both pristine and doped structures. The impurity states hybridize with the upwards Sn-like branch of the Dirac cone. Panel (b) confirms the band inversion in bulk pristine and Sb-doped SnTe.

Fig. 8 (a) Electronic structure of the 16-layer slabs for the undoped, Sb1 and Sb2 structures, along with the projections on Sn (green), Te (blue) and Sb (orange) p-orbitals. (b) Bulk band structure of the pristine and doped supercell. The color scheme is the same as in panel (a).

Fig. 9 depicts PDOS(k_x,k_y,E) and spin texture maps of SnTe above the Dirac point. Two pairs of Dirac cones appear, which are folded into a ×-shaped fourfold pattern at the Γ point for the c(4 × 4) supercell. The M_z(k,E) component is omitted due to the absence of out-of-plane spin texture.

Fig. 9 PDOS(k_x,k_y) and maps of the magnetization components M_x(k_x,k_y) and M_y(k_x,k_y) of pristine SnTe at 0.05 eV. Top panel shows the unfolded 1 × 1 unit cell, demonstrating the appearance of four helical edge states. Bottom panel shows the c(4 × 4) folded unit cell; the previous helical states arrange around Γ in a fourfold pattern. White and black shades account for positive and negative values of M_i, respectively.

Fig. 10 shows the persistence of the pristine SnTe spin texture when impurities are included. The helicity is unaltered independently of the impurity position, as mentioned in Fig. 3. Additionally, it can be seen in Fig. 4 that the spin texture changes sign upon moving above or below the Dirac point.

Fig. 10 PDOS(k_x,k_y) and maps of the magnetization components M_x(k_x,k_y) and M_y(k_x,k_y) at an energy 20 meV below (left panels) and above (right panels) the Dirac point. White and black shades account for positive and negative values of M_i, respectively.

Appendix B

Surface states in the quadratic Dirac equation

The surface states used in the analytic calculations are obtained considering a semi-infinite slab extended for z > 0 and with periodic boundary conditions in the x- and y-directions. DFT calculations use a periodical supercell of size L_x × L_y; in order to match the same conditions of the ab initio structure we consider periodic boundary conditions in x- and y-directions while the z coordinate is unbounded in performing the integrals of eqn (4).

Due to the breaking of translational symmetry, we solve the Hamiltonian (3) with k_z → −i∂_z, using the ansatz of the form

ψ ∼ e^−λz exp[i(k_xx + k_yy)]Φ; (9)

in other words, we assumed plane waves in the in-plane directions, an exponential decay form in the z-direction and a constant spinorial part Φ.

We found that two decay lengths satisfy the eigenvalue equation

By combining solutions of the form (9) in order to fulfill the hard-wall boundary conditions at z = 0, we obtain two surface states with linearly dispersive energy given by E_± = ε ± vk_∥sign(B) with wave functions

ψ₊ = A_s(e^−λ₁z − e^−λ₂z)exp[i(k_xx + k_yy)]Φ⁰₊, (11)

ψ₋ = A_s(e^−λ₁z − e^−λ₂z)exp[i(k_xx + k_yy)]Φ⁰₋, (12)

where Φ⁰_± are space-independent normalized spinors and the normalization constant depends on the momenta and model parameters

The numerical values of the parameters have been obtained by fitting to DFT data and are reported in Table 1.

Table 1 Fitting parameters of the continuum model to DFT data

SnTe

ε = −0.065 eV

m = 0.07 eV

v = 1.591 eV Å

B = 2.74 eV Ų

---

In eqn (3), replacing the electron free mass by the effective mass and the vacuum dielectric constant by that of the medium, the following relation for the effective Rydberg and Bohr radius can be easily obtained

where R_y and a_B are the well-known atomic Rydberg energy and Bohr radius and ε_r is the relative dielectric permittivity. The value of the permittivity in SnTe has been tabulated in multiple works, such as ref. 51 and 52 for a wide range of temperatures. In our DFT simulations, as already mentioned in the previous sections, we are in a room-temperature regime. Hence, the tabulated value of ε ≈ 40 at room temperature is the more convenient.^51,53

Green's function in the absence of impurity

When describing the envelope function χ(r), the momentum perpendicular to the growth direction k_‖ is conserved and the envelope function can be factorized as χ(r) = [small chi, Greek, tilde](z)exp(ir_‖·k_‖). It is understood that the subscript ‖ in a vector indicates the nullification of its z-component, namely r_‖ = (x,y,0) and k_‖ = (k_x,k_y,0).

The Green's function for the unperturbed problem (λ = 0) satisfies the following equation

The Green's function G₀(r,r′;E) = 〈r|Ĝ₀(E)|r′〉 is the matrix element of the unperturbed resolvent Ĝ₀(E) in coordinate representation. Let Ĝ(E) be the resolvent of the Hamiltonian , where is the operator associated to the non-local separable potential (7). Ĝ(E) and Ĝ₀(E) are related by the well-known equation (see ref. 54), where

is the transition operator. A pole of at a real (complex) energy E corresponds to a bound (resonance state). After some algebra we get the following transcendental equation for the poles of the transition operator where G₀(r,r′;E) can be expanded aswithand κ² = 1 + k_‖² − E².

The proper Green's function of the unperturbed band-inverted heterojunction can be obtained as follows. We introduce (18) in eqn (15) to obtain

The resolvent of ĥ₀ is Ĝ₀(E) = (E1₄ − ĥ₀)⁻¹. We now define the auxiliary resolvent ĝ₀(E) = (E²1₄ − ĥ₀²)⁻¹ so that Ĝ₀(E) = (E1₄ + ĥ₀)ĝ₀(E). Therefore, once ĝ₀(E) is known, we can easily obtain Ĝ₀(E). In coordinate representation

This equation can be easily solved regarding the term −2iβα_xδ(x) as an interaction potential and then using the Dyson equation. In doing so we obtain

Finally, recalling that Ĝ₀(E) = (E + ĥ₀)ĝ₀(E), we get eqn (19).

Once the basic equations are presented, we now need to take a particular shape function to perform the calculation. In what follows we consider

ω(r − r₀) = f(z − z₀)δ(r_∥), (23)

where f(z) is a top-hat function of width L and height 1/L, centered at z₀ and approaching the δ-function limit (L ≪ 1). The transcendental eqn (17) for the poles reduces to

det[1₄ − M(E)] = 0, (24a)

We have numerically solved eqn (24a) considering a finite bandwidth. This is equivalent to introducing an upper cutoff for the in-plane momentum, k_c. In this way, we prevent divergences when taking the limit L → 0, i.e. f(z − z₀) → δ(z − z₀).

Acknowledgements

We thank Gloria Platero for generously sharing her computational resources and Sergio Bravo and Álvaro Díaz-Fernández for helpful discussions. This work was supported by Ministerio de Economía y Competitividad, Spanish MCIN and AEI and the European Union under Grants PGC2018-097018-B-I00 (MCIN/AEI/FEDER, UE), PRE2019-088874 funded by MCIN/AEI/ 10.13039/501100011033 and by “ESF Investing in your future”, and PID2019-106820RB-C21.

References

1. M. Z. Hasan and C. L. Kane, Rev. Mod. Phys., 2010, 82, 3045.
2. X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys., 2011, 83, 1057.
3. A. Manchon, H. C. Koo, J. Nitta, S. M. Frolov and R. A. Duine, Nat. Mater., 2015, 14, 871.
4. Y. Ando and L. Fu, Annu. Rev. Condens. Matter Phys., 2015, 6, 361.
5. B. J. Wieder, B. Bradlyn, J. Cano, Z. Wang, M. G. Vergniory, L. Elcoro, A. A. Soluyanov, C. Felser, T. Neupert, N. Regnault and B. A. Bernevig, Nat. Rev. Mater., 2021, 7, 196.
6. T. H. Hsieh, H. Lin, J. Liu, W. Duan, A. Bansil and L. Fu, Nat. Commun., 2012, 3, 982.
7. P. Dziawa, B. J. Kowalski, K. Dybko, R. Buczko, A. Szczerbakow, M. Szot, E. Łusakowska, T. Balasubramanian, B. M. Wojek, M. H. Berntsen, O. Tjernberg and T. Story, Nat. Mater., 2012, 11, 1023.
8. Y. Tanaka, Z. Ren, T. Sato, K. Nakayama, S. Souma, T. Takahashi, K. Segawa and Y. Ando, Nat. Phys., 2012, 8, 800.
9. S.-Y. Xu, C. Liu, N. Alidoust, M. Neupane, D. Qian, I. Belopolski, J. Denlinger, Y. Wang, H. Lin and L. A. Wray, et al. , Nat. Commun., 2012, 3, 1.
10. J. Liu, W. Duan and L. Fu, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 88, 241303.
11. Z.-Y. Ye, H.-X. Deng, H.-Z. Wu, S.-S. Li, S.-H. Wei and J.-W. Luo, npj Comput. Mater., 2015, 1, 15001.
12. Y. Okada, M. Serbyn, H. Lin, D. Walkup, W. Zhou, C. Dhital, M. Neupane, S. Xu, Y. J. Wang, R. Sankar, F. Chou, A. Bansil, M. Z. Hasan, S. D. Wilson, L. Fu and V. Madhavan, Science, 2013, 341, 1496.
13. F. Schindler, A. M. Cook, M. G. Vergniory, Z. Wang, S. S. P. Parkin, B. A. Bernevig and T. Neupert, Sci. Adv., 2018, 4, eaat0346.
14. J. O. Dimmock, I. Melngailis and A. J. Strauss, Phys. Rev. Lett., 1966, 16, 1193.
15. H. Ozawa, A. Yamakage, M. Sato and Y. Tanaka, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 90, 045309.
16. J. Liu, T. H. Hsieh, P. Wei, W. Duan, J. Moodera and L. Fu, Nat. Mater., 2013, 13, 178.
17. A. L. Araújo, G. J. Ferreira and T. M. Schmidt, Sci. Rep., 2018, 8, 9452.
18. K. Hoang, S. D. Mahanti and M. G. Kanatzidis, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 115106.
19. X. Zhao, X. Zhang, T. Wang, C. Xia and S. Wei, J. Mater. Sci., 2018, 53, 15995.
20. T. M. Schmidt and G. Srivastava, Comput. Mater. Sci., 2020, 182, 109777.
21. Q. Qu, B. Liu, H. Liu, J. Liang, J. Wang, D. Pan and I. K. Sou, Nanoscale, 2021, 13, 18160–18172.
22. M. Yarmohammadi and K. Mirabbaszadeh, Phys. Chem. Chem. Phys., 2019, 21, 21633.
23. J. Wang, W. Powers, Z. Zhang, M. Smith, B. J. McIntosh, S. K. Bac, L. Riney, M. Zhukovskyi, T. Orlova, L. P. Rokhinson, Y.-T. Hsu, X. Liu and B. A. Assaf, Nano Lett., 2022, 22, 792–800.
24. C. Fang, M. J. Gilbert and B. A. Bernevig, Phys. Rev. Lett., 2014, 112, 046801.
25. L. T. Hoa, T. C. Phong and B. D. Hoi, Phys. Chem. Chem. Phys., 2020, 22, 13613–13621.
26. J. K. Hulm, C. K. Jones, D. W. Deis, H. A. Fairbank and P. A. Lawless, Phys. Rev., 1968, 169, 388.
27. J. M. Soler, E. Artacho, J. D. Gale, A. García, J. Junquera, P. Ordejón and D. Sánchez-Portal, J. Phys.: Condens. Matter, 2002, 14, 2745.
28. A. García, N. Papior, A. Akhtar, E. Artacho, V. Blum, E. Bosoni, P. Brandimarte, M. Brandbyge, J. I. Cerdá, F. Corsetti, R. Cuadrado, V. Dikan, J. Ferrer, J. Gale, P. García-Fernández, V. M. García-Suárez, S. García, G. Huhs, S. Illera, R. Korytár, P. Koval, I. Lebedeva, L. Lin, P. López-Tarifa, S. G. Mayo, S. Mohr, P. Ordejón, A. Postnikov, Y. Pouillon, M. Pruneda, R. Robles, D. Sánchez-Portal, J. M. Soler, R. Ullah, V. W.-Z. Yu and J. Junquera, J. Chem. Phys., 2020, 152, 204108.
29. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
30. R. Cuadrado and J. I. Cerdá, J. Phys.: Condens. Matter, 2012, 24, 086005.
31. N. Troullier and J. L. Martins, Phys. Rev. B: Condens. Matter Mater. Phys., 1991, 43, 1993.
32. J. Cerdá, M. A. Van Hove, P. Sautet and M. Salmeron, Phys. Rev. B: Condens. Matter Mater. Phys., 1997, 56, 15885.
33. E. T. R. Rossen, C. F. J. Flipse and J. I. Cerdá, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 87, 235412.
34. Y. J. Wang, W. F. Tsai, H. Lin, S. Y. Xu, M. Neupane, M. Z. Hasan and A. Bansil, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 87, 235317.
35. S. Shen, Topological Insulators: Dirac Equation in Condensed Matter, Springer, Berlin, 2017.
36. C.-X. Liu, X.-L. Qi, H. Zhang, X. Dai, Z. Fang and S.-C. Zhang, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 045122.
37. S.-Q. Shen, W.-Y. Shan and H.-Z. Lu, SPIN, 2010, 1, 33.
38. H. Zhang, C.-X. Liu, X.-L. Qi, X. Dai, Z. Fang and S.-C. Zhang, Nat. Phys., 2009, 5, 438.
39. F. Zhang, C. L. Kane and E. J. Mele, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 86, 081303.
40. R. Knox, Theory of excitons, Academic Press, Cambridge, 1963.
41. H. Ibach and H. Lüth, Solid-State Physics: An Introduction to Principles of Materials Science, Springer, Berlin, 2009.
42. M. Glasser, Surf. Sci., 1977, 64, 141.
43. B. W. Knight and G. A. Peterson, Phys. Rev., 1963, 132, 1085.
44. P. R. Sievert and M. L. Glasser, Phys. Rev. B: Solid State, 1973, 7, 1265.
45. S. López and F. Domínguez-Adame, Semicond. Sci. Technol., 2002, 17, 227.
46. F. Domínguez-Adame, B. Méndez, E. Maciá and M. A. González, Mol. Phys., 1991, 74, 1065.
47. F. Domínguez-Adame, E. Diez and A. Sánchez, Phys. Rev. B: Condens. Matter Mater. Phys., 1995, 51, 8115.
48. C. González-Santander, T. Apostolova and F. Domínguez-Adame, J. Phys.: Condens. Matter, 2013, 25, 335802.
49. V. G. Kantser and N. M. Malkova, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17074.
50. Y. L. Chen, J.-H. Chu, J. G. Analytis, Z. K. Liu, K. Igarashi, H.-H. Kuo, X. L. Qi, S. K. Mo, R. G. Moore, D. H. Lu, M. Hashimoto, T. Sasagawa, S. C. Zhang, I. R. Fisher, Z. Hussain and Z. X. Shen, Science, 2010, 329, 659.
51. H. Riedl, J. Dixon and R. Schoolar, Solid State Commun., 1965, 3, 323.
52. M. Cardona, J. Appl. Phys., 1965, 36, 2181.
53. S. Adachi, Optical Constants of Crystalline and Amorphous Semiconductors: Numerical Data and Graphical Information, Springer, Berlin, 1999.
54. E. N. Economou, Green's functions in quantum physics, Springer, Berlin, 2006.

---

Footnote

† Deceased.

---