Bi2Se3 topological insulator quantum wires

We study the effect of 3D topological insulator contributions to the band structure and wavefunctions of quasi-one-dimensional electron systems. Our model for this system consists of an effective Hamiltonian derived previously in the literature from the crystal symmetries of Bi2Se3. We find that in wires whose face lies in the plane formed by the y and z crystal directions and whose width is around 25 quintuple layers or more, the bands nearest to the gap are non-monotonic; we show that this has implications for the conductance of the wire. In addition, we observed increasing penetration depth of surface states with increasing wavenumber of the propagating mode. We believe these results have qualitative relevance to the family of 3D topological insulators whose crystal structure is characterised by the space group D 3 d 5 , and that the work done here contributes to the wider field of the study of conductance in topological insulators.


Introduction
Topological insulators (TIs) are a recently identified quantum state of matter with bulk states characteristic of semiconductors or insulators and conducting gapless surface states characteristic of a Dirac metal with strong Rashba-like spin-orbit coupling [1][2][3][4]. Recent advances in the growth of thin-film ultra high quality topological insulator material promise to enrich the physics of quantum wires through the introduction of a number of topology-related phenomena. Such thin films of TIs have been the subject of recent study for the purposes of determining the relationship between the bulk and surface states [5][6][7]. Quantum wires fabricated from gated or etched thin-film TI material will have two inequivalent perpendicular surfaces. The physics of these systems combines quantised surface conduction of Dirac electrons or holes, with Klein tunnelling between inequivalent surfaces, spin-orbit coupling effects and electron-electron interactions. In this paper we will not consider electron-electron interaction effects directly, other than by admitting that there will be some form of lateral effective potential. We will concentrate on a specific, well known, TI Bi 2 Se 3 that has bulk states with a semiconducting gap arising from single-particle band structure effects and three inequivalent surfaces in the x, y and z directions [8][9][10][11][12]. Our results are therefore specific to this material but they will be observable to some degree in any TI material, especially those where interaction effects are weak and the material has a single Dirac cone. Examples of similar approaches to that taken in this work can be found in [13,14], amongst many others.
Topological insulator wires have seen much study experimentally and theoretically. The Aharonov-Bohm effect in topological nanowires was carefully studied in [15]. Bi 2 Se 3 nanowires have been studied as field-effect transistors [16]; there has also been considerable interest in the presence of Majorana modes in TI wires [17]. More recently, work has been done specifically in the area of disordered quantum wires (see [18][19][20]). The model Hamiltonian used throughout the paper, an effective four-band continuum Hamiltonian previously used to capture the fundamental physics of surface states in TIs, has also seen considerable study [21][22][23][24][25][26][27][28]. This model is presented in section 2. Also in section 2, we quantise 2 components of the wavevector k, discuss the discretization of the problem to make it suitable for numerical computation and formulate the corresponding Hamiltonian and spin operators. We present our results in section 3, where we explain the physical origin of our dispersion relations and eigenstates. The band structure is used to calculate the ballistic conductance of a long Q1DES. Conclusions are made in section 4.

Model
For this work, we use the Hamiltonian presented by Liu et al [29], which has its origin in a tight binding model of the hybridized and spin-split Bi and Se p-orbitals, as well as the symmetries of the crystal. Here we retain only terms up to k 2 ; thus we arrive at the following effective Hamiltonian: 2 . The constants used in these defintions are defined in table 1. Γ i (   i 1 5) represent the 4×4 Dirac matrices, and V(y, z) is zero in the simulation domain and infinite at its edges. We present the derivation and atomic origin of this Hamiltonian in appendix A. Figure 1 illustrates the device that we apply this Hamiltonian to. In order to apply the Hamiltonian(1) to a finite system instead of the infinite system limit, the substitutions k z →−i∂ z and k y →−i∂ y are made in the Schrödinger equation representing quantisation in the z and y directions respectively. The stationary states of the system are denoted by F( ) y z , , and the components with respect to the basis of the Hamiltonian by {f i (y, z) }   ( ) i 1 4 which depend only on y and z by virtue of the systemʼs translational invariance in the x direction. Equation (3) then becomes  where x x 0 the index j=y, z and the matrices  ( ) k x 0 and  ( ) i j may be found in equations (28) to (32). We solve equation (4) on the surface  defined w.l.o.g. by x=0, and   y L 0 ,   z L 0 . We impose Dirichlet boundary conditions that Φ=0 on  ¶ , the boundary of  , since the potential V is set to be infinite outside the wire. The surface is split up into (n+1)×(n+1) lattice sites, separated by a distance h=L/(n+1) in both y and z. We then substitute the z-and y-derivatives for (central) finite difference operators with step length h [30] (defined in the appendix)to produce n 2 sets of 4 coupled equations for {f i } at each site lattice site of the form , , 1 4 and the 4 ×4 matrices  m n , can be found in equations (33) to (35).
if both μ and ν are nonzero. The boundary conditions on  ¶ give us that Φ 0, j =Φ n+1,j =Φ i,0 =Φ i,n+1 =0. The equations may then be written as an ( n 4 2 -dimensional) eigenvalue problem = † is a sparsely populated Hermtian band matrix, which gives rise to real eigenenergies as required. The z index (1 to n) of the ith element of F is given by 1and the y index by + (⌊ ⌋ ) i n mod 4 , 1. Numerical methods were used to find the eigenvalues and eigenvectors of the Hamiltonian matrix H [31], which we denote as E i and F ( ) i respectively. The same procedure can be carried out with the substitutions k x →−i∂ x and k y →−i∂ y such that the wire lies instead along theẑ direction, and the surface over which we solve, ¢, lies in the x-y plane. We still obtain n 2 sets of 4 simultaneous equations similar to (6), but the matrices are now given by  ¢ m n , . The states are labelled as F¢ ñ | ( ) i with energies ¢ E i . We expect the states to exhibit more symmetry due to the symmetry between k x and k y in the model Hamiltonian(2) as  k is invariant under the transformation « k k x y . Naturally, we should also be able to substitute k x →−i∂ x and k z →−i∂ z to produce results similar to the case where we quantise in y and z, but in this work we restrict ourselves to the y-z and x-y planes.
Within this model we also consider the action of the total angular momentum operatorsˆ( ) J x yz . Projecting the total momentum operators onto the eigenbasis of the Hamiltonian 2 ) which are also the eigenstates ofĴ z , we obtain , where {σ i } are the Pauli spin matrices and Ä is the matrix direct (Kronecker) product. Using the fact that , where † denotes Hermitian conjugation, upon transition from the continuous to the discrete case, the matrix representations of the total momentum operators for our n We find that each state is doubly degenerate such that any linear combination is also an eigenstate with the same eigenenergy. To find the states of definite angular momentum á ñ J i we lift the degeneracy by adding a very small electric field across the wire, large enough to break the degeneracy but small enough not to affect the energies significantly. The perturbation to the Hamiltonian takes the form with V 0 small (of the order 10 −9 eV), and a factor of 1 2 so that the degeneracy is lifted differently in y and z.

Dispersion relations
The band structure for various cross sections having quantised in theŷ andẑ can be found in figure 2. For smaller values of L, the bands approximately take the form of hyperbolae. For larger cross-sections, however, we see some anticrossing near the gamma point induced by spin-orbit coupling (SOC) which is characteristic of TIs. This is consistent with the ab initio calculations of Zhang et al [32] for Bi Se 2 3 . The SOC shifts the electron subbands down in energy and the hole subbands up, forcing the subbands to anticross. Without anticrossing, the band structure would presumably be a set of intersecting hyperbolae. This results in nonmonotonic bands with points of inflexion where  = k d d 0 . At these inflexion points the electrons and holes have a divergent effective mass Far away from ò=0, the bands return to being approximately hyperbolic. Even the bands near ò=0 are linear away from k x =0 where the electron and hole subbands are well separated. For 15 QL in figure 3(b) we see some 'unexpected' electron subbands with a lower curvature (at least as k x → 0, that is). These do not correspond to surface states-they are bulk states which so happen to have a low enough energy to appear within the pictured energy range.
We quote the result derived by Brey and Fertig [33] who give an analytic expression for the dispersion relation under the assumption that , such that the Hamiltonian exhibits particle hole symmetry For k x =0 and L z =L y =L equation (10) reduces to the simple expression which defines the analytic energy gap ΔE a . We compare the energy gaps observed in our numerically calculated dispersion relations against this analytic result. For our numerical results we define ΔE n ≡E n+1 −E n such that the topmost hole subband is given an index of n=0 and hence the band gap is given by ΔE 0 . We see from figure 2 that equation (11) gets the salient features of the dispersion relations correct. Away from ò=0, the bands are approximately hyperbolae. For each value of L the ratio ΔE 0 /ΔE a is of order ( ) 1 and the energy gaps ΔE n , at least for larger values of L, remain approximately constant.
When quantising instead in x and y we produce a different set of dispersion relations, which can be found in figure 3, as we expected since the crystal structure and thus (2) are not isotropic. The bands are, for all values of L, approximately hyperbolic, but only for larger values of L are the states equally spaced in energy, and hence consistent with the analytic results derived by Brey and Fertig in [33]. The electron and hole subbands are sufficiently well separated in energy for there to be no anticrossing. For smaller cross sections the dispersion relation is interrupted by low-energy bulk states.
As we traverse the band structure upwards, i.e.the order of the electron subbands N=1, 2, 3, K, we see that the number of nodes in F | | 2 in the horizontal (ŷ) direction increases as N−1, that is 0, 1, 2, K. This is a well known result in quantum mechanics. For each state in figure 4  x we see an increasing number of nodes but in the vertical (ẑ) direction. As for the electron subbands there exists a degenerate partner related via symmetry with opposite spin.

Surface states
We plotted the wavefunctions with Matlab using Delaunay triangulation (delaunay(x, y)) and interpolated shading (Figures 5, 6). For every cross section the states near the Γ point and the zero in energy are surface states.
Quantising in x and y, the states have a smaller penetration depth in theẑ direction than in theŷ direction, but the wavefunction along the y edges (z=0, L) has a greater magnitude. At the corners of the domain, 'constructive interference' between the two surfaces at the y and z edges causes a peak in each corner. Taking a quadrant ( defined by, say, y>L/2 and z>L/2) of the domain we may deduce whether an electron is more likely to be found along a z or a y surface by splitting the quadrant in two along its diagonal (the line y=z) and evaluating  F ∬ | | r d 2 2 . For 25 QL it turns out that it is more probable for an electron to be found along the y edge (53%), while for 40 QL the z edge is more likely.
When quantising instead in k x and k y , the surface states exhibit a much greater degree of symmetry, but otherwise contain the same features. For 8 QL (and below), the states from both surfaces overlap to produce a non-zero amplitude at the centre of the wire. To find the corresponding penetration depth we must discretise an inner product along a thin strip of surface at the edge in order to calculate the expectation value of y from our states F ( ) i , namely  where Φ is the unnormailsed wavefunction and h=L/(n+1) is the distance between adjacent lattice sites. A similar expression exists for á ñ z with the substitution ih→jh. This expression in most appropriate for states which possess little symmetry about the centre of the cross section. However, a more appropriate measure of penetration depth for a more symmetrical wavefunction is the difference between the width L and the root mean square deviation from the centre, º -á -ñ ℓ (ˆ) L y L 2 y 2 and å å å å where again the relationship for ℓ z is obtained by making the substitution ih→jh. For a symmetrical surface state equation (12) would give a penetration depth of zero. For the bottommost electron subband, (N=1), we look at the effect of increasing k x on the wavefunctions. Counter-intuitively, the penetration depth á ñ y (these states are not symmetric) increases with increasing k x . If k x were representative of the systemʼs total kinetic energy, one would expect the states to exhibit greater curvature (∇ 2 Φ), and thus a smaller penetration depth. Resorting to classical reasoning, it is possible that increasing the systemʼs momentum in thex direction (k x ) results in increased penetration depth as SOC may be represented by an effective magnetic field, which gives rise to a Lorentz force = - ( ) e F v B , which acts towards the centre of the wire. We plot á ñ y against k x in figure 7 for 25 and 50 QL.

Spin and angular momentum
The expectation value of the angular momentum operatorsĴ x andĴ y , á ñ J x and á ñ  eigenvalues having the same magnitude. With inversion k kthe energy ordering of these two states is reversed(when the degeneracy is lifted). The degeneracy can be observed in figures 8 and 9, where each band is symmetric under inversion k x → −k x . Had we not lifted the degeneracy of these states using a linear electric field, we would have obtained random linear combinations of both degenerate states. This degeneracy is an example of Kramerʼs pairs.

Ballistic conduction
Following Moroz and Barnes [34], we now use our dispersion relations to calculate the Ballistic Conductance G of a long Q1DES at low temperature (Figures 10, 11).
where  F is the Fermi energy and M(ò F ) is the number of occupied electron subbands. M(ò F ) is given by where θ(x) is the unit step function and  ( ) n i , s min represents the energy of the ith minimum of the nth subband. We assume that equation (14) holds true which necessitates conservation of (charge) current j. Moroz and Barnes in [34] proved that current is indeed conserved in the presence of finite SOC. To find G from equation (14) one needs only to draw a horizontal line for each Fermi energy and count the number of intersections with the band structure which gives the value of G in units of G 0 .
We found in section 3.1 that for larger widths L the dispersion relations exhibit points of inflexion. The current passing through a 1D subband is not affected by the diverging effective mass when there is no scattering mechanism [34]. However between the two points of inflexion there exist regions where for k x >0 and, equivalently, v x >0 for k x <0. This does significantly affect G(ò) since both forward and backward travelling electron modes contribute and we see a corresponding peak in G(ò). As we move further away from ò=0 where SOC becomes increasingly negligible, the points of inflection disappear along with the regions over which the group velocity is negative and no peak in G(ò) is observed. This is also the case for smaller values of L where the effect of SOC is not sufficiently strong.

Conclusions
We have found using numerical techniques the band structure and wavefunctions for Bi 2 Se 3 TI quantum nanowires of square cross section. Quantising in y and z we observed anticrossing non-monotonic bands for larger widths than previous analytic analyses; far away from the gap, the bands return to being approximately hyperbolic. Quantising instead in the more symmetrical x and y we observed well separated hyperbolic, monotonic bands, consistent with previous literature [33].  We observed the change in the wavefunction with electron subband order for some fixed k and found it exhibits increasing numbers of nodes with increasing subband order, as might be expected from standard quantum mechanics. For all possible quantisations we found surface states near the Γ point, whose penetration depth increases with increasing wavenumber. We observed constructive interference between surface states at the device corners, and showed that the electronʼs 'preferred' edge (namely, the one with the higher probability density) changes with wire dimension in the case where the wire cross-section directions are along the y-and z-directions of the crystal.
In addition to investigating the energetic structure and penetration depth of the surface states of the TI, we have also discussed the impact on the ballistic conductance of a Q1DES. We used the conservation of current in the presence of finite SOC to calculate the ballistic conductance G of a long TI nanowire directly from the energy spectrum, using the Landauer formula, and observed peaks in G(ò) due to points of inflexion in the bandstructure. We made some predictions about the location and magnitude of these peaks compared to a usual  quantised conductance measurement. The work done in this paper, particularly on the subject of conduction, is relevant to potential experimental probes such as those in [35,36].
It is important to note that the Liu Hamiltonian, used through this paper, is derived from crystal symmetries with parameters from · k p theory fitted to ab initio calculations. Thus, results derived using such a Hamiltonian should carry over to a lesser or greater extent to any crystal with similar symmetry. This gives our results some  is not met [29]).

Acknowledgments
The work of A N was supported by the EPSRC. The authors would like to extend their thanks to Oliver Hart for calculations that supported this paper.

Appendix A. Derivation of the Liu Hamiltonian
The following is a summary of work done by Liu et al [29]. Bi 2 Se 3 has a layered structure, with each layer consisting of Be or Si atoms arranged in hexagonal lattices which are stacked in the A-B-C-A-B-C-Kconfiguration. The unit cell contains five atoms total with two equivalent Se atoms (labelled Se1 and Se1′) and two equivalent Be atoms (Be1 and Be1′). Each set of 5 layers is known as a quintuple layer (QL) with » Å 1 QL 9.5 . Strong chemical bonds hold together atoms within one QL, while adjacent QLs are held together by weak Van der Waals forces, so only the bonding within one QL is considered first. The hybridised atomic states associated with Bi and Bi′ ( ñ a |B , ¢ ñ a |B respectively), and Se and Se′ ( ñ a |S , ¢ ñ a |S respectively) are combined to form states with parity eigenvalues ±1 (upper index) where α=p x , p y , p z . añ -|P1 , and añ + |P2 , are too far away from the Fermi energy to give a significant contribution to wavefunctions and so are neglected from this point onwards. We now introduce a basis formed of the eigenstates of the angular momentum operatorL z and include spin s =   , , which now includes spin-orbit coupling, is expressed in the basis of the 12 remaining states. The matrix is block diagonal and so may be diagonalised by individually diagonalising each block. We require the solution to the eigenvector problem defined by the matrix and ΔE Λ =(E Λ,x −E Λ,z −λ Λ /2)/2. Since ΔE Λ and λ Λ are real, u Λ ± and  L v are real also. We obtain the new states which diagonalise the Hamiltonian with SOC included: L ñ = L ñ + L ñ . These states were found to be closest to the Fermi energy and hence give the most significant contribution to the wavefunction.

Appendix B. Submatrices of the Hamiltonian
We define the following matrices