Bulk and interface quantum states of electrons in multi-layer heterostructures with topological materials

In this article we describe the bulk and interface quantum states of electrons in multi-layer heterostructures in one dimension, consisting of topological insulators (TIs) and topologically trivial materials. We use and extend an effective four-band continuum Hamiltonian by introducing position dependence to the eight material parameters of the Hamiltonian. We are able to demonstrate complete conduction-valence band mixing in the interface states. We find evidence for topological features of bulk states of multi-layer TI heterostructures, as well as demonstrating both complete and incomplete conduction-valence band inversion at different bulk state energies. We show that the linear kz terms in the low-energy Hamiltonian, arising from overlap of pz orbitals between different atomic layers in the case of chalcogenides, control the amount of tunneling from TIs to trivial insulators. Finally, we show that the same linear kz terms in the low-energy Hamiltonian affect the material’s ability to form the localised interface state, and we demonstrate that due to this effect the spin and probability density localisation in a thin film of Sb2Te3 is incomplete. We show that changing the parameter that controls the magnitude of the overlap of pz orbitals affects the transport characteristics of the topologically conducting states, with incomplete topological state localisation resulting in increased backscattering.


Introduction
Topological insulators (TIs) are a quantum state of matter identified a little over ten years ago [1]. They are commonly distinguished from trivial insulators by conducting states associated with surfaces or interfaces; the dispersion relation for these conducting states generally obeys the Dirac equation [2,3]. Key to this classification is the fact that the existence of these states is guaranteed by the bulk properties of the topological insulator relative to the vacuum or more recognisable 'trivial' insulator on the other side of the interface [4,5].
The electronic structure of TI junctions has been studied previously both theoretically and experimentally, often with the objective of designing spintronics-type devices or bandstructure engineering towards different states of matter [6][7][8].
In this work we study the quantum states and dispersion relations of multiple TI junctions in a single heterostructure, by modifying an effective four-band continuum Hamiltonian previously used to capture the fundamental physics of surface states in TIs [9][10][11][12][13][14][15][16]. We apply this to a number of structures with alternating topological insulator and trivial material layers, including conductors. Because of quantum-mechanical tunneling by interface states in TIs to adjacent trivial material layers, which also results in hybridisation between interface states in proximity to each other, we have interface states that are not fully localised inside one material. This requires us to extend the effective four-band continuum Hamiltonian to take into account the position dependence of the eight material parameters it includes.
We consider a set of substitutions in the four-band Hamiltonian that allow us to write down a system of position-dependent Schrödinger equations connecting the four bands described by the Hamiltonian. The substitutions we use are similar to those performed in the effective mass approx imation for semiconductor junctions-they connect low-energy continuum Hamiltonians for one or more bands derived using k · p perturbation theory to a system of Schrödinger equations that contain position-dependent carrier masses. These substitutions are derived from making use of an envelope function approximation [17,18]. We demonstrate how to numerically solve the system of positiondependent Schrödinger equations we obtain by making use of the finite difference approximation. Our position-dependent effective Hamiltonian applies to real systems and materials, but in this work we illustrate it by making use of a series of simulated materials defined only by the eight material parameters of the Hamiltonian, as well as making use of values of those parameters that are sometimes used in the description of the low-energy behaviour of commonly studied 3D TIs Bi 2 Se 3 and Sb 2 Te 3 [13][14][15]19].
We aim to demonstrate the expected interface states between topological and trivial insulators, and we show that these states naturally tunnel into adjacent trivial insulator layers and hybridise with other interface states. We show the evolution of trivial bulk states into the spin-split interface states with changing topological phase. We also investigate the band mixing and inversion in both interface and bulk states, and look for possible topological features of bulk states. Additionally, we investigate the effect that the linear k z term in the low-energy Hamiltonian, which represent p z orbitals overlaps between different atomic layers in chalcogenide crystal structures, has on a material's ability to support localised interface states regardless of its topological phase. We investigate this effect using model parameters corresponding to real materials.
This paper is organised as follows. Section 2 describes the position-dependent effective Hamiltonian we use and discusses our numerical approach to solving our effective position-dependent Hamiltonian. In section 3 we show the electronic states and dispersion relations of multilayer heterostructures of topological and trivial materials as a function of changing material parameters, demonstrating a series of new physical effects. Section 4 contains our conclusions.

Effective Hamiltonian
The Hamiltonian H used for bulk materials here is a standard 4-band model, based on that used by Liu et al [13] which can be written as (ignoring the k 3 terms) where I is the 4 × 4 identity matrix, The constants (C 0 , C 1 , C 2 , M 0 , M 1 , M 2 , A, B) in these definitions are material-dependent parameters; the interpretation of these parameters arises from k · p theory and their specific values for a particular material can be extracted from ab initio calculations or experiments. Broadly speaking, C(k) represents the overall elliptical dispersion relation, M 0 is the bulk gap at the Γ point which can be positive or negative, and A and B are Rashba-type spin-orbit interaction parameters. M 1 and M 2 describe the effective masses of the quasiparticles formed by the hybridised atomic orbitals at k = 0, which form the basis set in which this Hamiltonian is written. This Hamiltonian describes topologically distinct phases for M 0 > 0 and M 0 < 0 [13,15,20] (see also appendix B), and this is the parameter we use in our simulations to continuously change our structures from TI to non-TI and back 1 .
Since we are going to be looking only at k x dispersions in this paper, we will neglect all the k y terms in what follows. These can easily be added back in and treated the same way as k x terms. It is important to note that this Hamiltonian has a spin degeneracy in its energy bands. This is due to its overall timereversal invariant nature, a key feature of Z 2 -invariant TIs [4].
We are interested in obtaining solutions to Hamiltonian (1) in systems that are bulk-like in two of three dimensions (x, y), and contain a series of interfaces in the third dimension (z) represented by infinite 2D planes in the x and y dimensions. We follow the literature (see for example [21][22][23]) in performing the standard substitution in equation (1) to obtain the time-independent Schrödinger equation In order to simulate the changing of material parameters across an interface, we will need to extend equation (1). The approach to this problem pursued here is based on the envelope function approximation (EFA), commonly in use in traditional semiconductor physics [24,25]. In semiconductors, the EFA leads to the effective mass approximation, which has been succesfully used to analyse semiconductor heterostructures by means of a single-band position-dependent Schrödinger equation [17], where the effective mass of the carriers varies with position. In line with the effective mass approximation in semiconductors, our position-dependent equation can be obtained from effective k · p Hamiltonians such as equation (1) by substituting k → −i∇, and writing in the position dependence for our material parameters (C 0 (z), C 1 (z), C 2 (z), M 0 (z), M 1 (z), M 2 (z), A(z), B(z)). We will also need to make changes to preserve the hermiticity of the model-these are discussed in the next section.

Numerical approach
We need to rewrite equation (4) after performing the substitution in equation (3) such that it is Hermitian in regions where the material parameters are changing. This approach is similar to that historically taken in other materials e.g. wurtzite semiconductor heterostructures [26], and has also been published recently elsewhere specifically for TIs [27]. For this we use the literature standard form of the kinetic energy (see [17,18,24,25]). The relevant terms that need to be rewritten are where the curly braces indicate the usual anticommutation relation. We now integrate equation (4) with respect to z from a point − just to the left of the interface at z = 0 to a point ε just to the right of the interface. Using equations (5) and (3) where In order to solve it numerically, we would like to turn equation (6) into a discretised equation. For the purposes of integration we begin by discretising the position variable z such that where i is our new lattice index variable and ∆z is the lattice spacing. We introduce the wavefunction values Ψ i at each point We perform the integral on the right-hand side of equation (6) by setting = ∆z 2 to obtain where the superscript (0) denotes that H 0 uses the material parameters for site 0, which here is the interface site. Because we assume that the junction is symmetric, the interface site has material parameter values equal to the average of the material parameter values either side of the interface. We evaluate the integral on the left-hand side of equation (6) by parts and evaluate it at − ∆z 2 and ∆z 2 by interpolation using the approximations This finally allows us to rewrite equation (6) as where The superscript (i) denotes the value of the material parameter for site i. Note that equation (12) is exactly the same as a normal second-order finite difference equation if no change in material parameters occurs across the interface-see appendix A. Equation (12) can also be extended to cover the application of a constant magnetic field in the y-direction by making use of the Peierls substitution k x → k x + e B y z, where B y is the field strength. Before our results we briefly cover the numerical detail of our simulations. Other than in section 3.6, we modelled 15 nm heterostructures using 150 lattice points (1 Å lattice spacing), with interface sites at lattice points 31 and 119. The results in section 3.6 were obtained using 100 lattice points to simulate 10 nm, with no interface sites as we simulated a single slab of Sb 2 Te 3 . Thicknesses were chosen to minimise the amount of coupling between interface or edge states, such as that

Parameter Value
described in [12]-this can be seen wherever we plot bandstructures with Dirac cones, where the gap in the cone that would exist due to coupling is small enough to be invisible.

Emergence of spin-polarised heterostructure interface states tunneling into adjacent trivial insulators
We have stated previously that the Hamiltonian (1) has distinct topological phases. These phases are controlled by the material parameter M 0 , which is equal to the Γ-point energy gap between the two electronic states nearest to the Fermi level. For topologically non-trivial materials, this parameter is negative, and it indicates a crossing over of bands has occured. In Bi 2 Se 3 , this is the gap between the lowest bonding orbital of the hybridised Bi p-orbitals (labelled |P1 + − , ±1/2 ), and the highest antibonding orbital of the hybridised Se p-orbitals (labelled |P2 − + , ±1/2 )-see [13]. In a multi-layer heterostructure, we can observe the emergence of the interface states by changing the parameter M 0 in some part of the structure. Beginning with a heterostructure of topologically-trivial materials, each with a positive gap M 0 , we would expect to see a change in total probability distribution of the lowest conduction and highest valence subbands as M 0 decreases and becomes negative in part of the heterostructure. Figure 1 shows how they change from bulklike states that resemble closely an electron trapped in a potential well, to states localised towards the material interfaces. Also in figure 1, we show how the subbands of the full heterostructure change as the gap M 0 closes and becomes negative. The conduction and valence bands move towards each other until they meet, and as they do they lose their parabolic character and adopt a Dirac-like linear dispersion relation about the crossing point-see figures C1 and C2 in appendix C. Figure 1 also shows how the resulting interface states are spin-polarised-once the gap M 0 becomes negative enough, there is almost no overlap between the spin-up and spin-down interface states. This is one of the reasons for the lack of backscattering mechanisms for electrons travelling along these interface states-in order for such an event to occur, both the crystal momentum and the spin of the electron must be reversed. In addition, it is clear that even before the change of topological phase there is some separation of the spin channels-this is due to the presence of the spin-orbit interaction, which is controlled by the parameters A and B in the Hamiltonian (1).

Evidence of complete mixing of conduction and valence orbitals in the topological interface state
The states presented in figure 1 are shown as probability distributions in space over all four basis states of the Hamiltonian. However, we would not necessarily expect the electronic occupation of each basis state of the Hamiltonian to be equal in each case.
One important point about the topological interface states in a finite-size device such as the ones studied in this work, is that they are an example of one of many subbands that are formed in a finite device. These subbands are distinct from the bulk bands that are described at the start of section 3.1, labelled |P1 + − , ±1/2 and |P2 − + , ±1/2 . However, this section, and section 3.4, are concerned with the four components of a particular subband; those four components can be traced back to the bulk states |P1 + − , ±1/2 and |P2 − + , ±1/2 by examining the form of the numerical Hamiltonian (12).
In figure 2, we plot the occupation of each basis orbital as a function of position in the heterostructure, for one of the two highest valence orbitals. We find that, for the topologically trivial cases where the gap M 0 is positive throughout the structure, the electrons in the highest valence band occupy primarily the |P2 − + , ±1/2 states. We find this to be an expected resultbefore the band inversion caused by the spin-orbit interaction, the P2 orbitals represent the highest energy orbital that is still below the Fermi level [13]. As the gap M 0 is made negative for the part of the heterostructure indicated on figure 2, we see for the highest valence orbital that the |P1 + − , ±1/2 states begin to match the |P2 − + , ±1/2 states in terms of occupation. It is clear that the surface states are formed of a mixture of the two sets of states, as the spin orbit interaction in the crystal brings them together and inverts them, changing the topological phase of the material.

Evidence of topological features of heterostructure bulk states
We turn our attention now to the bulk states just above and just below the bands that form the interface states. We define here bulk states as those states which are associated with bands that remain parabolic or quartic throughout our parameter sweep; these are states with probability distributions which are spread further away from the material interface. These states also experience some changes in probability distribution and orbital composition as the relevant part of the heterostructure changes topological phase. Figures 3 and 4 plot the evolution of the bulk valence bands closest to the interface states. We can see that at the highlighted point k x = 0.2nm −1 on the bandstructure, the probability distribution of the resulting state shifts with the change in topological phase from a bimodal distribution to one with a single maximum. Both the initial and the final states shown are reminiscent of the states of an electron in a potential well, but we have moved from two nodes to one node. The quanti zed electron in a potential well has states such that each extra node-each additional maximum in the probability distribution-represents an additional unit of energy. We can see that where we have lost the lowest energy, zero-node distribution state in figure 1 as it forms into an interface state, we have gained a new zero-node distribution state in figures 3 and 4. However, the shapes of these distributions are not identical. We can see from figures 3 and 4 compared to figure 1 that the zero-node bulk state has an additional kink towards the edge of the heterostructure, near to the interface (highlighted in the figure). It appears likely that this bulk state is affected by the topological phase trans ition in this area.
It is also clear from figures 3 and 4 that there is some spin polarisation across the heterostructure in these bulk states. Again, this is due to the spin orbit interaction in these materials, and the non-zero crystal momentum of these states. We would expect to see similar effects in all of the states created by our simulations.

Evidence of complete and incomplete band inversion in heterostructure bulk states of different energies relative to the gap
We can also demonstrate bulk band inversion within the heterostructure by examining the composition of the bulk state in terms of the Hamiltonian basis states-the |P1 + − , ±1/2 bonding states and the |P2 − + , ±1/2 antibonding states, as discussed previously. Figure 5 shows the probability distribution of one of the highest energy valence bulk states, plotting the magnitude of each of the four components. We can see that much like in the case of the interface state that arises from the valence band, the state starts of dominated by the |P2 − + , ±1/2 antibonding states. As the gap M 0 of the structure is changed, we find that occupation of the |P1 + − , ±1/2 bonding states increases, and by the time the structure has reached a negative enough gap, the |P1 + − , ±1/2 bonding states become the dominant component of this valence bulk state. While this image is not included here, we find the same kind of inversion from |P1 + − , ±1/2 bonding states to |P2 − + , ±1/2 antibonding states takes place in the lowest parts of the conduction band.
However, further away from the gap the inversion is less complete, and in fact no inversion occurs in the trivial insulator parts of the heterostructure. We illustrate this phenomenon in figure 6. Inside the topological insulator, the bulk state contains a mixture of the |P1 + − , ±1/2 bonding and |P2 − + , ±1/2 antibonding states, while inside the trivial insulators the bulk state contains only |P2 − + , ±1/2 antibonding states. This fits with the expected ordering of the bonding and antibonding states in trivial and TIs. It is also to be expected that the further away a band is from the gap, the less affected it will be by the topological crossover of the bonding and antibonding states.

Demonstrating control of interface state tunneling into adjacent trivial insulator by p z -orbital interactions
In the devices studied so far, we have been focusing on the effect of changing the primary topological parameter, the band gap M 0 . We have been studying the topological phase transition controlled by this parameter. However, there are more subtle effects taking place at the TI-trivial insulator boundary that are controlled by other parameters in the Hamiltonian. In this section and the next we consider the parameter B, which is found from k · p theory to be [13]  In the absence of the spin-orbit interaction, this parameter is zero [25]. From equation (14) and by comparison with Hamiltonian (1), we can see that the B parameter is responsible for mixing |P1 + − , ±1/2 bonding and |P2 − + , ±1/2 antibonding states. Having seen in previous sections that the TI-trivial insulator interface state is a mixed state of the bonding and antibonding states, we would expect to see some impact on the interface state from variations of the B parameter.
In figure 7 we plot the total probability density of one of the interface states across a trivial-TI-trivial heterostructure for four different values of B. We find in this figure that the B parameter controls the depth of tunneling of the interface state into the trivial insulator, separately from the bandgap of the trivial insulator. Decreasing the B parameter in the trivial insulator relative to the topological insulator decreases the tunneling depth of the interface state into the trivial insulator. Similarly, increasing the B parameter in the trivial insulator will increase the amount of tunneling.
On an atomic level, the overlap in equation (14) is dominated by the interaction of the p z orbitals of Bi and Se respectively (in the case of Bi 2 Se 3 , or their equivalent atoms in other chalcogenides and alloys of chalcogenides). These will be affected by the atomic lattice spacing and the relevant outmost shells of the atoms in question (6p for Bi, 5p for Sb and Te, 4p for Se). We would expect, for example, that the overlap with the 5p orbitals of Sb atoms would be smaller than the 6p orbitals of Bi, which have probability density spread much further away from the nucleus. We explore this prediction in the next section.

Evidence of incomplete localisation of surface states of Sb 2 Te 3 due to reduced p z -orbital interactions
So far we have used models with particle-hole symmetry to examine effects inside TI-trivial insulator heterostructures. Here, we employ parameters for a known material, Sb 2 Te 3 , to identify new features in its surface states. We established in the previous section that the atomic p z orbital overlap between the different layers of the chalcogenide affects the formation of the interface state. This overlap is controlled in Hamiltonian (1) by the parameter B. [13] gives material parameter values for each of three materials, which we reproduce in table 2.
By comparing the entries for M 0 and B for the two materials, we can see that the spin-orbit coupling that inverts the bonding and antibonding bands is weaker in Sb 2 Te 3 than in Bi 2 Se 3 . With the lattice vectors of Sb 2 Te 3 being similar to those of Bi 2 Se 3 [19,28], one possible suggestion for the weaker p z overlap is that the much closer distribution of the Sb 5p z orbital compared to the Bi 6p z orbital, as well as the reduced spin orbit interaction due to decreased mass, is not offset by the increase in ionic mass and orbital diffusion by moving from the Se ion and 4p z orbital to the Te ion and 5p z orbital. Figure 8 shows the low-energy bandstructure of Sb 2 Te 3 as reproduced by Hamiltonian (1) using the parameters from table 2, and the probability distribution and spin expectation value of the surface state that is indicated on the bandstructure. We can see that the surface state is no longer confined around the edge, but has nodes further into the core of the TI slab. In fact, there are spin up and spin down nodes overlapping in the centre of the structure, which creates the possibility for backscattering in the edge states. As discussed above, we believe this is the result of the weaker overlap of Sb and Te 5p z orbitals relative to the same overlaps in Bi 2 Se 3 . This is an important result, because it suggests that the effectiveness of topologically protected surface conduction can be undermined by atomic-level features separate to the topological phase of the material. We investigate this further in the next section.

Backscattering measurements
We make use of our simulated wavefunctions from solving equation (1) to calculate the effect on the backscattering in interface states due to the change in the B parameter, which represents the dipole matrix element between the p z orbitals of the different atomic layers of the chalcogenide. We consider two modes of longitudinal transport in the x − y plane, through a cross-section of a device such as the one illustrated in figure 7 which is finite in the z-dimension. For a Fermi level close to the Dirac point, conduction takes place through two sets of spin-locked interface states (see figures 1 and 8). Each spin-locked set has a forward and backward moving mode associated with it (+k x or −k x ). Scattering between the different spin states is not possible [1], so we focus on two modes-forward and backwardwith the same spin.
To measure the backscattering, we follow [29] in making use of the Green's function in the form where z and z are defined at the terminals y = 0 and y = 0, v is the electron velocity, and Ψ(z) is the subband which solves equation (4). From the Green's function we make use of the Fisher-Lee relation where q and p label device terminals, to find the backscattering matrix element q = p (in terms of equation (15), y = y ). Figure 9(b) shows the magnitude of the backscattering matrix element s pp for the interface state (labelled |s i |) plotted as a fraction of the magnitude of the backscattering matrix element for the bulk state (labelled |s b |), against the value of the parameter B, for the device shown in figure 9(a), using material parameters in table 1. We can see clear evidence of an increase in backscattering as the value of B decreases, though this is clearly not linear.

Conclusion and outlook
In this article we have investigated the bulk and interface quantum states of electrons in multi-layer heterostructures of topological and trivial materials. We made use of an effective four-band continuum Hamiltonian by treating the eight material parameters of the Hamiltonian as position-dependent.
We have shown that, as expected, spin-polarised interface states form in the heterostructure, and tunnel a finite distance into trivial insulators. We have presented evidence of complete mixing of conduction and valence bands in interface states. We have shown that bulk states in TI heterostructures exhibit features not present in trivial insulator heterostructures, and that there exists a varying degree of conduction and valence band mixing in bulk states according to energy. We have shown that the p z orbital overlap between neighbouring atomic layers in a chalcogenide can affect the tunneling of the interface state into a neighbouring trivial insulator. Finally, we present evidence that the surface states of Sb 2 Te 3 may not be completely localised at the edges, due to a lower amount of p z orbital overlap in that crystal relative to the more closely studied TI of the same family, Bi 2 Se 3 . We showed that this could potentially impact the transport through the surface states due to creating the possibility of backscattering.
Our system relied on realistic physical parameters drawn from fitting of the Hamiltonian (1) to ab initio calculations in Bi 2 Se 3 and Sb 3 Te 3 [13]. Other materials could be fitted using the eight parameters of the model, and our analysis could be compared to other studies of those materials. In addition, [13] presents extensions of Hamiltonian (1) that take into account the hexagonal warping of the Dirac cone highly prominent in materials such as Bi 2 Te 3 [19], which can also be included in our analysis by extending the finite difference approximation to higher order derivatives. It would be very interesting to see ab initio calculations that model a sequence of alloys with relatively similar material properties but changing B param eter (p z orbital dipole matrix element) in order to compare to our results-this is something we are interested in pursuing in a further publication. The approach we have outlined here is powerful enough to apply to multiple interfaces across various different materials as required, and thus should be of interest to anyone in the field looking for a first approach to modelling a heterostructure of mixed topological phase.