††thanks<EMAIL_ADDRESS><EMAIL_ADDRESS>
# Low-energy modeling of three-dimensional topological insulator
nanostructures

Eduárd Zsurka  Peter Grünberg Institute (PGI-9), Forschungszentrum Jülich,
52425 Jülich, Germany  JARA-Fundamentals of Future Information Technology,
Jülich-Aachen Research Alliance, Forschungszentrum Jülich and RWTH Aachen
University, 52425 Jülich, Germany  Department of Physics and Materials
Science, University of Luxembourg, 1511 Luxembourg, Luxembourg  Cheng Wang
Peter Grünberg Institute (PGI-1), Forschungszentrum Jülich and JARA, 52425
Jülich, Germany  Julian Legendre  Department of Physics and Materials Science,
University of Luxembourg, 1511 Luxembourg, Luxembourg  Daniele Di Miceli
Department of Physics and Materials Science, University of Luxembourg, 1511
Luxembourg, Luxembourg  Institute for Cross-Disciplinary Physics and Complex
Systems IFISC (CSIC-UIB), E-07122 Palma, Spain  Llorenç Serra  Institute for
Cross-Disciplinary Physics and Complex Systems IFISC (CSIC-UIB), E-07122
Palma, Spain  Department of Physics, University of the Balearic Islands,
E-07122 Palma, Spain  Detlev Grützmacher  Peter Grünberg Institute (PGI-9),
Forschungszentrum Jülich, 52425 Jülich, Germany  JARA-Fundamentals of Future
Information Technology, Jülich-Aachen Research Alliance, Forschungszentrum
Jülich and RWTH Aachen University, 52425 Jülich, Germany  Thomas L. Schmidt
Department of Physics and Materials Science, University of Luxembourg, 1511
Luxembourg, Luxembourg  Philipp Rüßmann  Peter Grünberg Institute (PGI-1),
Forschungszentrum Jülich and JARA, 52425 Jülich, Germany  Institute for
Theoretical Physics and Astrophysics, University of Würzburg, 97074 Würzburg,
Germany  Kristof Moors  Peter Grünberg Institute (PGI-9), Forschungszentrum
Jülich, 52425 Jülich, Germany  JARA-Fundamentals of Future Information
Technology, Jülich-Aachen Research Alliance, Forschungszentrum Jülich and RWTH
Aachen University, 52425 Jülich, Germany

###### Abstract

We develop an accurate nanoelectronic modeling approach for realistic three-
dimensional topological insulator nanostructures and investigate their low-
energy surface-state spectrum. Starting from the commonly considered four-band
$\boldsymbol{\mathrm{k\cdot p}}$ bulk model Hamiltonian for the Bi2Se3 family
of topological insulators, we derive new parameter sets for Bi2Se3, Bi2Te3 and
Sb2Te3. We consider a fitting strategy applied to _ab initio_ band structures
around the $\Gamma$ point that ensures a quantitatively accurate description
of the low-energy bulk and surface states, while avoiding the appearance of
unphysical low-energy states at higher momenta, something that is not
guaranteed by the commonly considered perturbative approach. We analyze the
effects that arise in the low-energy spectrum of topological surface states
due to band anisotropy and electron-hole asymmetry, yielding Dirac surface
states that naturally localize on different side facets. In the thin-film
limit, when surface states hybridize through the bulk, we resort to a thin-
film model and derive thickness-dependent model parameters from _ab initio_
calculations that show good agreement with experimentally resolved band
structures, unlike the bulk model that neglects relevant many-body effects in
this regime. Our versatile modeling approach offers a reliable starting point
for accurate simulations of realistic topological material-based
nanoelectronic devices.

##  I Introduction

Topological insulators (TIs) are a novel class of materials that have garnered
substantial interest in the recent decades due to their possible application
in electronics, spintronics and quantum information processing [1, 2]. TIs are
characterized by the existence of topologically protected states at the
boundaries of a sample, which are protected against any local perturbations
that respect time-reversal symmetry. The Bi2Se3 family of materials, here
referring to Bi2Se3, Bi2Te3 and Sb2Te3, are three-dimensional (3D) time-
reversal-invariant TIs with a large inverted gap. They have a layered
structure, consisting of five-atom, or quintuple layers (QL) arranged along
the $\hat{z}$ direction (see Fig. 1a). The bulk electronic structure is
described by a nontrivial $Z_{2}$ topological invariant, which ensures the
existence of protected spin-nondegenerate surface states with massless Dirac-
like dispersion.

The Bi2Se3 family of 3D TIs have been first described in 2009, when their
topological properties were uncovered and a four-band
$\boldsymbol{\mathrm{k\cdot p}}$ Hamiltonian describing the bulk dispersion
was proposed [3]. Material parameters of this Hamiltonian were obtained using
perturbation theory [4, 5], which yields an accurate description of the
electronic band structure at $\Gamma$ at low energies. However, the degree to
which the vicinity of $\Gamma$ is captured by the obtained parameters varies
from case to case. When used for nanoelectronic device simulations, the model
has to accurately capture the entire region in momentum space over which the
Dirac cone of the topological surface states stretches out, while also
remaining well behaved at larger momenta such that unphysical electronic
states do not appear at low energies. This is not always guaranteed by
applying perturbation theory at $\Gamma$, which motivates us to derive a new
set of parameters that can give a good quantitative description of the Bi2Se3
family of materials. In this work, we obtain the parameters of the four-band
$\boldsymbol{\mathrm{k\cdot p}}$ Hamiltonian for Bi2Se3, Bi2Te3 and Sb2Te3 by
an alternative method. We fit the model to _ab initio_ band structures in such
a way that the vicinity of $\Gamma$ is accurately considered up to
sufficiently large momenta and all the relevant features of the band structure
(e.g., the topology) are taken into account.

With the newly obtained material parameters, we analyze the low-energy
spectrum of experimentally relevant nanostructures. We consider the effects of
band structure anisotropy and electron-hole asymmetry, which modify the
dispersion of the surfaces states on surfaces with different orientations.
While usually ignored for resolving the transport properties related to Dirac
surface states in nanostructure systems [6, 7, 8, 9, 10], we find that the
low-energy spectrum can be significantly affected by this anisotropy and
asymmetry for some of the materials under consideration here.

|  | Bi2Se3 [3] | Bi2Se3 [4] | Bi2Se3 [5] | Bi2Se3 fit | Bi2Te3 [4] | Bi2Te3 [5] | Bi2Te3 fit | Sb2Te3 [4] | Sb2Te3 [5] | Sb2Te3 fit
---|---|---|---|---|---|---|---|---|---|---|---
$A_{0}$ | $[\text{eV\AA}]$ | 4.1 | 3.33 | 2.51 | 4.33 | 2.87 | 4 | 4.40 | 3.4 | 3.7 | 3.89
$B_{0}$ | $[\text{eV\AA}]$ | 2.2 | 2.26 | 1.83 | 1.94 | 0.3 | 0.9 | 0.55 | 0.84 | 1.17 | 1.69
$C_{0}$ | $[\text{eV}]$ | -0.0068 | -0.0083 | 0.048 | -0.28 | -0.18 | -0.12 | -0.014 | 0.001 | 0.02 | 0.10
$C_{1}$ | $[\text{eV\AA}^{2}]$ | 1.3 | 5.74 | 1.41 | 1.46 | 6.55 | 2.67 | 1.65 | -12.39 | -14.2 | -6.48
$C_{2}$ | $[\text{eV\AA}^{2}]$ | 19.6 | 30.4 | 13.9 | 22.81 | 49.68 | 154.35 | 29.47 | -10.78 | -6.97 | -4.26
$M_{0}$ | $[\text{eV}]$ | -0.28 | -0.28 | -0.17 | -0.30 | -0.3 | -0.3 | -0.26 | -0.22 | -0.18 | -0.21
$M_{1}$ | $[\text{eV\AA}^{2}]$ | 10 | 6.86 | 3.35 | 6.00 | 2.79 | 9.25 | 4.62 | 19.64 | 22.12 | 19.32
$M_{2}$ | $[\text{eV\AA}^{2}]$ | 56.6 | 44.5 | 29.35 | 70.38 | 57.38 | 177.23 | 72.80 | 48.51 | 51.28 | 63.91
Bulk gap | [meV] | 560 | 344 | 280 | 472 | - | 155 | 303 | 135 | 155 | 303
$\lambda_{z}$ | $[\text{\AA}]$ | 9.01 | 12.83 | 9.87 | 6.0 | - | 19.68 | 15.8 | 36.28 | 28.99 | 21.48
$\zeta_{\text{DP}}$ |  | 0.57 | 0.92 | 0.71 | 0.62 | 1.67 | 0.64 | 0.67 | 0.18 | 0.18 | 0.31
$\Delta E_{\text{DP}}$ | $[\text{meV}]$ | 61 | 43 | 9 | 24 | 445 | 175 | 12 | 90 | 91 | 56
Table 1: Parameters for the $\boldsymbol{\mathrm{k\cdot p}}$ bulk model
Hamiltonian (1) for the three studied materials: Bi2Se3, Bi2Te3 and Sb2Te3.
The parameters taken from Refs. [4, 5] are obtained through perturbation
theory applied at $\Gamma$. We introduce a new set of parameters obtained by
fitting to _ab initio_ band structure data around $\Gamma$. Additionally, the
bulk gap, the penetration depth $\lambda_{z}$ along $z$ of the surfaces states
in the $x\text{-}y$ plane, the relative position $\zeta_{\text{DP}}$ of the
Dirac point of a slab parallel to the $x\text{-}y$ plane with respect to the
bulk gap at $\Gamma$ (see Eq. (5)), and the energy difference $\Delta
E_{\text{DP}}$ between the Dirac points of surfaces states parallel to the
$x\text{-}y$ and $x\text{-}z$ planes, are also included (see Sec. II for
details).

When the thickness of a nanoribbon approaches a few QLs, the surface states
localized on top and bottom surfaces can hybridize, leading to a gap opening
in the surface-state dispersion at $\Gamma$ [11, 12, 13, 14]. In this case, we
refer to the system as being in the _thin-film limit_. If the hybridization of
the surface states is accompanied by an inversion of the surface-state
spectrum, the system enters a quantum spin Hall insulator (QSHI) regime. Such
thin-film geometries have attracted significant interest, being suitable for
studying QSHI edge channels [15, 16], the quantum anomalous Hall effect [17,
18, 19, 20, 21, 22, 23, 24, 25, 26, 27], and topological superconductivity
[28, 29, 30, 31, 32].

According to early theoretical work on 3D TI thin-films, the gap at $\Gamma$,
which here we refer to as the hybridization gap $\Delta E_{\text{hyb}}$, shows
an oscillatory behavior between a QSHI and a normal insulator (NI) state when
the thickness of the thin-film is varied [11, 12, 14]. However, more recent
results suggest that many-body effects arising in the thin-film limit modify
the oscillations and the size of the hybridization gap, giving a better
agreement with experimental results [33, 34, 35, 36]. To describe this limit,
we employ an effective thin-film model that captures solely the surface-state
dispersion. We extract the material parameters of the thin-film model by
fitting to the surface-state spectra obtained from _GW_ calculations of thin
films, which take into account the relevant many-body effects, for thicknesses
ranging from 2 to 6 QL.

This paper is structured as follows. In Sec. II, we give an overview of the
models used to describe 3D TIs in the bulk, at the surface, and in the thin-
film limit, and also discuss the material parameters. In Sec. III, we present
the dispersion of quasi-one-dimensional nanostructures and analyze the effect
of anisotropy and electron-hole asymmetry. In Sec. IV, we treat the thin-film
limit using the bulk model and compare the results to experimental findings.
We also provide thickness-dependent material parameters for the effective
thin-film model. Finally, in Sec. V, we interpret our findings and discuss
other aspects that may be relevant for accurate nanoelectronic device modeling
with 3D TIs.

##  II Models

###  II.1 Bulk model

The low-energy bulk electronic structure of the Bi2Se3 family of materials
around $\Gamma$ ($\boldsymbol{k}=\boldsymbol{0}$) can be described using a
four-band model, where only the valence and conduction bands that are
responsible for the band inversion are considered. The bulk model Hamiltonian
can be written in the following form [3, 4],

$\displaystyle H_{\text{bulk}}(\boldsymbol{k})$ $\displaystyle=$
$\displaystyle\epsilon(\boldsymbol{k})+\mathcal{M}(\boldsymbol{k})\tau_{z}$
(1)
$\displaystyle+\,A_{0}(k_{y}\sigma_{x}-k_{x}\sigma_{y})\tau_{x}+B_{0}k_{z}\tau_{y},$
$\displaystyle\epsilon(\boldsymbol{k})$ $\displaystyle=$ $\displaystyle
C_{0}+C_{1}k_{z}^{2}+C_{2}\big{(}k_{x}^{2}+k_{y}^{2}\big{)},$
$\displaystyle\mathcal{M}(\boldsymbol{k})$ $\displaystyle=$ $\displaystyle
M_{0}+M_{1}k_{z}^{2}+M_{2}\big{(}k_{x}^{2}+k_{y}^{2}\big{)},$

with $\sigma_{i}$, $\tau_{i}$ ($i=x,y,z$) the Pauli matrices for the spin and
orbital degree of freedom and model parameters
$A_{0},B_{0},C_{0},C_{1},C_{2},M_{0},M_{1},M_{2}\in\mathbb{R}$ that can be
obtained from _ab initio_ calculations or perturbation theory. This bulk model
describes an insulator only when $|C_{1}|<|M_{1}|$ and $|C_{2}|<|M_{2}|$,
avoiding the closing of the band gap at large values of $|\boldsymbol{k}|$
[37], which is necessary to avoid the appearance of unphysical states at low
energies upon applying the finite-difference method. The band inversion and
consequently the topological properties are determined by $M_{0},M_{1}$ and
$M_{2}$ – only when $M_{0}M_{1}<0$ and $M_{0}M_{2}<0$ are satisfied, Eq. (1)
describes a topologically nontrivial system [4]. $C_{0},C_{1}$ and $C_{2}$ are
responsible for the electron-hole asymmetry, while $A_{0}$ and $B_{0}$ can be
interpreted as the group velocities (up to a factor of $\hbar$) of the surface
states on surfaces orthogonal to any in-plane ($x$-$y$) direction, and to the
$\hat{z}$ direction, respectively. The anisotropy of the band structure is
captured by different values for the corresponding in-plane and out-of-plane
terms $A_{0}\neq B_{0}$, $C_{1}\neq C_{2}$ or $M_{1}\neq M_{2}$. Finite values
of $M_{1}$ and $M_{2}$ prevent the fermion doubling problem, hence the
Hamiltonian (1) can be safely discretized on a lattice without acquiring
unphysical Dirac points (at e.g. $k_{x,y,z}=\pm\pi/a$ if one considers a cubic
lattice with lattice constant $a$ [38]), making the model suitable for
modeling the low-energy spectrum of TI nanostructures with arbitrary shapes
[9].

###  II.2 Material parameters

Our aim is to obtain a model that accurately describes the topological surface
states and can be discretized on a lattice using the finite-difference method
[39], without the appearance of unphysical states in the bulk band gap. As the
distinctive feature of 3D TIs is the Dirac cone of the topological surface
states, the region in $k$-space around $\Gamma$ where this Dirac cone appears
is of central importance. Thus, the material parameters used in Eq. (1) should
yield a dispersion that reliably describes the bulk bands, up to the momenta
where the Dirac cone merges with the bulk bands. The first full set of
parameters introduced for the Bi2Se3 family of materials in Ref. [4] were
obtained using perturbation theory, giving an accurate description of the low
energy states only very close to $\Gamma$ (up to
$k_{x,z}<0.04~{}\text{\AA}^{-1}$). In a later work, a
$\boldsymbol{\mathrm{k\cdot p}}$ perturbation approach applied to _ab initio_
calculations yielded a good qualitative description of the conduction (CB) and
valence bands (VB) of Bi2Se3 and Sb2Te3 [5]. However, for Bi2Te3, in the
region in $k$-space where the Dirac cone appears, the obtained dispersion has
a significantly smaller band gap than the one observed in _ab initio_ band
structures (see Figs. 1d,e). In our fitting procedure the eigenvalues of the
bulk model Hamiltonian (1) are fitted to _ab initio_ band structures of bulk
Bi2Se3, Bi2Te3 and Sb2Te3, with the band structure of Bi2Te3 shown in Fig. 1b.
More details on the _ab initio_ calculations are given in Appendix A. To
obtain the most accurate fit, we vary the region in $k$-space over which the
CB and VB are considered (imposing a minimal extent of the region to
accurately capture the Dirac cone in the bulk gap), and we maximize the
coefficient of determination $R^{2}$ of the fit. The conditions
$|C_{1}|<|M_{1}|$, $|C_{2}|<|M_{2}|$, $M_{0}M_{1}<0$ and $M_{0}M_{2}<0$ (see
Sec. IIA for more details) were enforced on the fitted parameters of the bulk
model. The resulting parameters are shown in Table 1 alongside the parameters
of Refs. [3, 4, 5]. For more details on the fitting procedure, see Appendix B.
We also evaluate the size of the band gap as the difference between the
minimum of the CB and the maximum of the VB. In Figs. 1d,e, we show the result
of our fit for Bi2Te3 to the relevant bands, together with the dispersion
obtained with parameters taken from Refs. [4, 5]. The band gap of the fitted
dispersion is much closer to that of the _ab initio_ calculation, than the
value obtained through perturbation theory.

Figure 1: (a) Crystal structure of the Bi2Se3 family of materials. (b) Band
structure for Bi2Te3. (c) Brillouin zone with space group $R\bar{3}m$. The
blue hexagon is the 2D Brillouin zone of the projected (1,1,1) surface, and
the high-symmetry points $\Gamma$, K and M are labeled. The dispersion along
(d) $k_{x}$ and (e) $k_{z}$ of the bulk model (orange) that is fitted to the
conduction and valence band of Bi2Se3 around the $\Gamma$ point, obtained from
_ab initio_ calculations (black dots), compared to the dispersion evaluated
with the material parameters taken from Ref. [4] (green) and Ref. [5] (cyan).
The upper limits of the fit in momentum space $k_{\parallel,z}^{\text{max}}$
(see Appendix B), are shown with gray dashed lines. Dispersion in $k_{x}$ of
the surface state of Bi2Te3, with the surface in the (f) $x\text{-}y$ and (g)
$x\text{-}z$ plane, evaluated analytically (dashed lines) and numerically
(solid lines), the hatched areas indicate the projection of the bulk bands.
The surface-state bands appear in the band gap and intersect each other at the
Dirac point, forming a Dirac cone (shaded region, marking the analytical
expression). Figure 2: (a),(b) Schematics of the nanostructures under
consideration: (a) a nanowire with $W\sim d$ and (b) a nanoribbon with $W\gg
d$. (c),(d) The spectrum of (c) a nanowire with $W=d=28$ nm and (d) a
nanoribbon with $W=50$ nm and $d=6$ nm, with the fitted material parameters
for Sb2Te3. The Dirac cone of the surface state in the $x\text{-}y$
($x\text{-}z$) plane (see Figs. 1f,g) is shown with blue (red) shading.
(e),(f),(g) The wave function density $|\Psi(y,z)|^{2}$ over the nanoribbon
cross section of states localized in different Dirac cones [see corresponding
symbol in (d)], indicated by the grayscale colormap (arbitrary units). Where
both Dirac cones overlap (purple shading), the states are delocalized over the
whole perimeter of the nanoribbon as shown in (g).

###  II.3 Surface-state model

When confined to a semi-infinite geometry, with a surface in the $x\text{-}y$
plane at $z=0$, the bulk model Hamiltonian (1) yields a gapless surface-state
spectrum, described by the following effective Hamiltonian [37, 13, 7],

$\begin{split}&H_{\text{surf}}^{\hat{z}}(k_{x},k_{y})=C_{0}-C_{1}M_{0}/M_{1}\\\
&\quad-\text{sgn}(M_{1})\sqrt{1-(C_{1}/M_{1})^{2}}A_{0}(k_{x}\sigma_{y}-k_{y}\sigma_{x}).\end{split}$
(2)

The wave function profile perpendicular to the $x\text{-}y$ surface of the
$k_{x}=k_{y}=0$ surface state has the following form when the 3D TI is
confined to the $z>0$ region [37, 4]

$\displaystyle\chi(z)$ $\displaystyle=$
$\displaystyle\begin{pmatrix}c_{1}&-c_{1}&c_{2}&c_{2}\end{pmatrix}^{T}(e^{q_{z}^{+}z}-e^{q_{z}^{-}z}),$
(3) $\displaystyle q_{z}^{\pm}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{2}\sqrt{\dfrac{B_{0}^{2}}{M_{1}^{2}-C_{1}^{2}}}\pm\sqrt{\dfrac{1}{4}\dfrac{B_{0}^{2}}{M_{1}^{2}-C_{1}^{2}}+\dfrac{M_{0}}{M_{1}}},$
(4)

with two independent parameters $c_{1}$ and $c_{2}$ (up to normalization). The
wave function extends into the bulk with a characteristic penetration depth
$\lambda_{z}=\text{max}\\{1/\mathfrak{R}(q_{z}^{+}),1/\mathfrak{R}(q_{z}^{-})\\}$
[9], which is listed in Table 1 for the different sets of material parameters.
In Fig. 1f the solutions of the Hamiltonian (2) (dashed lines) are shown,
along with the numerically evaluated spectrum of a semi-infinite slab (black
lines), which we obtain by using a version of the bulk Hamiltonian (1)
discretized on a lattice with confinement along $z$ and translational
invariance along $x$ and $y$. This solution is for a surface orthogonal to the
$\hat{z}$ direction, and analogous solutions can be obtained for surfaces with
other orientations. In Fig. 1g the solutions for a surface over the
$x\text{-}z$ plane at $y=0$ is given (dashed lines), together with the
numerical result (black lines). For simplicity, we set $C_{0}=0$ in all
calculations, since this term only yields an overall shift of the spectrum in
energy.

Note that in general, while barely noticeable in Figs. 1f,g, there can be a
shift in energy between the Dirac points (DPs) of the surface states, i.e.,
the energy where the surface-state bands cross each other at $\Gamma$. In the
Dirac-like dispersion given by Eq. (2), the DP can be defined as
$E_{\text{DP}}^{\hat{z}}=C_{0}-M_{0}C_{1}/M_{1}$. Surfaces orthogonal to the
$\hat{x}$ or $\hat{y}$ direction will host surface states with the DP
positioned at $E_{\text{DP}}^{\hat{x}/\hat{y}}=C_{0}-M_{0}C_{2}/M_{2}$. We
define the difference between the two DPs as $\Delta
E_{\text{DP}}=|E_{\text{DP}}^{\hat{x}/\hat{y}}-E_{\text{DP}}^{\hat{z}}|$. For
the studied materials (Bi2Se3, Bi2Te3, Sb2Te3), we obtain $\Delta
E_{\text{DP}}=24,12$ and $56~{}\text{meV}$ from the fitted parameters,
respectively, as also shown in Table 1. When the electron-hole asymmetry is
ignored ($C_{1}=C_{2}=0$), or the parameters are considered to be isotropic
($C_{1}=C_{2}$ and $M_{1}=M_{2}$) one naturally obtains $\Delta
E_{\text{DP}}=0$ ($E_{\text{DP}}^{\hat{z}}=E_{\text{DP}}^{\hat{x}/\hat{y}}$).
However, as the material parameters in Table 1 show, anisotropy and electron-
hole asymmetry can be significant.

Another important consequence of anisotropy and electron-hole asymmetry is
that the DP is not centered in the middle of the bulk gap. This can also be
seen in _ab initio_ calculations of slab geometries or angle-resolved
photoemission spectroscopy measurements [40, 41, 42]. In Table 1 we evaluate
the relative position of the DP for a surface in the $x\text{-}y$ plane as

$\zeta_{\text{DP}}=\frac{E_{\text{DP}}^{\hat{z}}-E_{\text{VB}}(\Gamma)}{E_{\text{CB}}(\Gamma)-E_{\text{VB}}(\Gamma)}=\frac{1}{2}\bigg{(}1-\text{sign}(M_{0})\frac{C_{1}}{M_{1}}\bigg{)},$
(5)

where $E_{\text{VB(CB)}}(\Gamma)$ is the energy of the VB (CB) at $\Gamma$. A
DP centered between the VB and CB at $\Gamma$ yields $\zeta_{\text{DP}}=0.5$,
while a DP at the top (bottom) of the VB (CB) at $\Gamma$ corresponds to
$\zeta_{\text{DP}}=0~{}(1)$.

###  II.4 Thin-film model

In the thin-film limit, with surfaces in the $x\text{-}y$ plane at $z=0$ and
$z=d$ (see Fig. 2b), tunneling between the surface states on the top and
bottom surfaces can open a finite _hybridization gap_ in the Dirac cone of the
surface states. In this limit, the dispersion of the system can be captured
using a low-energy thin-film model Hamiltonian [14, 13, 37],

$\displaystyle H_{\text{tf}}(k_{x},k_{y})$ $\displaystyle=$ $\displaystyle
E_{0}-Dk_{\parallel}^{2}+\hbar v_{\text{F}}(k_{y}\sigma_{x}-k_{x}\sigma_{y})$
(6) $\displaystyle+\,(\Delta/2-Bk_{\parallel}^{2})\sigma_{z}\tau_{z},$

with $\sigma_{i}$ ($i=x,y,z$) still the Pauli matrices for spin,
$k_{\parallel}^{2}\equiv k_{x}^{2}+k_{y}^{2}$, and $\tau_{z}$ acting on a
different subspace from the one before in Eq. (1), with eigenvalues $\pm$
representing a hyperbola index that distinguishes between the doubly-
degenerate surface-state solutions of Eq. (6). The Hamiltonian of Eq. (6) is
equivalent (ignoring the term $-Dk_{\parallel}^{2}$) to the four-band
effective model for a two-dimensional QSHI proposed by Bernevig, Hughes and
Zhang [43], which has been shown to capture the behavior of the bulk model of
Eq. (1) in the thin-film limit [12]. The hybridization gap of the surfaces
states is $\Delta E_{\text{hyb}}=2\Delta$ in this model, while the gap is
trivial if $\Delta B<0$, and the system is in the nontrivial inverted regime
if $\Delta B>0$, yielding a QSHI state [43, 14]. A non-zero value of $D$
results in electron-hole asymmetric surface states and $v_{\text{F}}$ is the
Fermi velocity of the surface states. It should be noted that $|D|<|B|$ is
required, otherwise there is no band gap at large $|\mathbf{k}|$ [14, 37].

##  III Nanostructures

In this section, we consider the bulk model Hamiltonian (1) to study the
dispersion of two relevant 3D TI nanostructures: nanowires, with approximately
equal width and thickness ($W\sim d$), and nanoribbons, for which $W\gg d$,
shown schematically in Figs. 2a,b. We use Kwant [44] to obtain the spectra of
the nanostructures, and Adaptive [45] for efficient parameter sampling. Note
that we describe thicknesses in terms of QL and nm interchangeably, as 1 QL
$\approx$ 1 nm for the materials under consideration.

Experimentally, 3D TI films down to a few-QL thickness have been achieved in
all three materials [40, 46, 41, 42, 47, 48, 49] and nanoribbons with widths
down to 50 nm were realized [50, 51]. Here we consider nanostructures with
negligible hybridization between surface states on opposing sides (for the
thin-film limit where hybridization becomes relevant, see Sec. IV). In Fig.
2c, we present the dispersion of a nanowire with square cross section and
$W=d=28$ nm long edges (black lines). In Fig. 2d, we also present the
dispersion of a nanoribbon with the same perimeter as the nanowire, but a much
larger width-to-height ratio, $W=50$ nm and $d=6$ nm. Here we have chosen
Sb2Te3 because it has the largest value of $\Delta E_{\text{DP}}$ for the
fitted parameters. The obtained spectra qualitatively resemble a conventional
Dirac cone with confinement quantization [6, 52]. However, there is a clear
difference between the dispersions of the nanowire and the nanoribbon. The
differences can be attributed to the effect of the surfaces of the
nanostructures that are oriented in different directions, which we explain
below.

To understand the effect of the different sides, in Figs. 2c,d, we overlay the
Dirac cones of the top/bottom (left/right side) surfaces, centered around
their respective DPs, with a blue (red) shading. The states appearing in the
spectrum of the nanostructures can be divided into three groups: states with
energies in the blue regions, extending over the top and bottom surfaces (Fig.
2e); states with energies in the red regions that are localized on the side
surfaces (Fig. 2f); and states that wrap around the whole perimeter of the
nanowire in the regions where both Dirac cones overlap (Fig. 2g). As the side
surfaces of the nanoribbon have a much smaller area than the top and bottom
surfaces, the side surfaces do not host as many states as in the case of the
nanowire. Hence, the spectrum of the nanoribbon resembles more closely a
quantized Dirac cone as expected for the top and bottom surface.

A relevant energy scale in the nanostructure is the spacing of subbands
originating from the confinement of the surface states to the finite perimeter
$P$ of its cross section. In the case of a Dirac dispersion, the spacing of
the subbands can be approximated by $2\pi v_{\text{F}}/P$. For the nanoribbon
considered in Fig. 2d, one obtains $2\pi v_{\text{F}}/P=17~{}\text{meV}$ (here
we consider an isotropic Fermi velocity $\hbar v_{\text{F}}\approx
3~{}\text{eV\AA}$). The spacing of the subbands has to be compared to the
effect of the surfaces of the nanoribbon being oriented in different
directions, which can be quantified using the energy difference of DPs $\Delta
E_{\text{DP}}$, which is 56 meV in the case of Sb2Te3. From this observation,
we can deduce that the difference in DP energy on different sides may even
affect the surface-state spectrum of the smallest attainable nanostructures.

##  IV Thin-film limit

We first investigate the thin-film limit in Sec. IV.1, using a discretized
version of the bulk model Hamiltonian (1) and we compare the results to
experimental findings, with an emphasis on the topology and size of the
hybridization gap. Second, in Sec. IV.2, we consider a different approach to
capture the physics of the thin-film limit quantitatively. Material parameters
for the thin-film model Hamiltonian (6) are obtained by fitting the model to
band structure data of _GW_ calculations that accurately describe the thin-
film limit [33, 34, 53].

###  IV.1 Treatment with the bulk model

When described using the bulk model of Eq. (1), it was shown that the
hybridization gap in the thin-film limit oscillates in size as the thickness
is varied. The closing and subsequent reopening of hybridization gap occurs at
certain critical thicknesses $d_{\text{c}n}$ (see Fig. 3b). The low-energy
physics of the system at a thickness close to $d_{\text{c}n}$ can be described
by the BHZ model [12, 43], which implies that the oscillation of the
hybridization gap is also accompanied by topological phase transitions, with
the system alternating between a NI and a QSHI phase.

In a 3D TI thin film of thickness $d$ the interference of two transverse wave
functions located on the top and bottom surfaces (orthogonal to $\hat{z}$),
given by Eqs. (3) and (4), will close the hybridization gap if the thickness
matches the critical value $d_{\text{c}n}=nl_{\text{c}}$, where
$l_{\text{c}}\simeq\pi/\mathfrak{Im}(q_{z}^{\pm})$ is the period of the
oscillation [13]. When the out-of-plane Fermi velocity $B_{0}$ is set to zero,
the period of the oscillation is equal to
$l_{\text{c}}(B_{0}=0)=\pi\sqrt{|M_{0}/M_{1}|}$. As shown in Fig. 3 for
Bi2Te3, when $B_{0}$ takes on a non-zero value, the value of $l_{\text{c}}$
increases from $l_{\text{c}}(B_{0}=0)$. However, if $B_{0}\gtrsim
2\sqrt{-M_{0}(M_{1}^{2}-C_{1}^{2})/M_{1}}$, the imaginary part of
$q_{z}^{\pm}$ goes to zero, the oscillatory behavior vanishes and the system
is a NI for all thicknesses.

In Fig. 4, for Bi2Se3, Bi2Te3 and Sb2Te3, the numerically evaluated
hybridization gap $\Delta E_{\text{hyb}}$ of the surface states at $\Gamma$ is
given as a function of the film thickness, where we use a discretized version
of the bulk model Hamiltonian (1). For the different parameters sets of Eq.
(1), we also indicate the topological phase for every value of $d$ with blank
(filled) shading below the curve for the NI (QSHI) phase. We compare the size
of the hybridization gap of the bulk model Hamiltonian (1) with values
determined experimentally from angle-resolved photoemission spectroscopy
measurements [40, 41, 42]. For the topology of the gap we consult results
obtained with the _GW_ method that accurately capture the many-body effects of
the surface states [33, 34, 53]. The size of the hybridization gap as
determined from the _GW_ calculations is also shown. The exact values of
$\Delta E_{\text{hyb}}$ determined with the different methods can be found for
2-6 QL thicknesses in Table 2, with a blank (grey) background indicating that
the thin film is a NI (QSHI).

Figure 3: (a) The real and imaginary part of $q_{z}^{\pm}$ (see Eq. (4)) as a
function of the out-of-plane Fermi velocity $B_{0}$. The real part is
inversely proportional to the penetration depth of the surface state, while
the imaginary part of $q_{z}^{\pm}$ determines the topological phase at a
given thickness $d$. (b) When described using the bulk model (1) the
topological phase of the 3D TI slab oscillates between a NI and QSHI as the
thickness $d$ is varied. The period of this oscillation increases as a
function of $B_{0}$ and diverges above a critical value of $B_{0}$ such that
the slab is always in the NI state. Figure 4: Size of the hybridization gap
$\Delta E_{\text{hyb}}$ at the $\Gamma$ point for thicknesses $2$ nm $<d<7$ nm
for (a) Bi2Se3, (b) Bi2Te3 and (c) Sb2Te3. The gap is extracted from the
spectra of the surface states in a slab geometry described using the bulk
model (1), with parameters taken from Refs. [3, 4, 5] and our fits. At a given
thickness the system is a QSHI (NI) if the area below the curve is colored
(blank). The hybridization gap obtained in the $GW$ calculations [33, 34] are
included, with the QSHI states marked with a circle. For reference,
experimentally determined gaps at the Dirac point are also shown [40, 41, 42].

According to experimental measurements and _GW_ calculations, it is expected
that Bi2Se3 is a NI in the thin-film limit [33, 36, 35]. The parameter set
from Ref. [3] and the fitted parameter set do not capture accurately this
behavior, rather suggesting that the thin film is in the QSHI above a
thickness of 3 QL. The parameter set listed in Refs. [4, 5] gives a better
description of Bi2Se3 thin films, yielding a NI for all thicknesses with
hybridization gaps which are systematically smaller than the experimentally
determined values, with Ref. [4] giving the closest values. Remarkably, the
hybridization gaps for Bi2Te3 match the experimentally determined values up to
4 QL very well when considering our parameter set obtained via fitting,
whereas the topology of the gap is also captured up to 5 QL. In contrast, the
same cannot be said of the parameters set of Ref. [5], while the parameter set
of Ref. [4] yields a gapless dispersion. In the case of Sb2Te3, all parameters
sets capture the topology of the hybridization gap accurately, but the size of
the hybridization gap deviates significantly from experimentally measured
values. From the obtained hybridization gaps, the thin film is NI for 2 and 3
QL and QSHI for 4, 5 and 6 QL, as predicted by _GW_ calculations. While
reasonable agreement is retrieved in certain cases, the results shown in Fig.
4 indicate that parameter sets, obtained with both perturbative and fitting
approaches, cannot be used to describe the thin-film limit reliably for the
different materials.

###  IV.2 Treatment with the thin-film model

As highlighted in the previous section, it is evident that the surface-state
dispersion obtained when confining the bulk model (1) to a slab geometry does
not reliably capture the electronic and topological properties of thin films
observed in experiment. However, recent theoretical works have shown that
employing the _GW_ method in determining the properties of thin films yields
band structure data that is in excellent agreement with angle-resolved
photoemission spectroscopy measurements [33, 34]. By fitting the low-energy
effective model (6) to this band structure data, we can obtain material
parameters for the Bi2Se3 family of materials in the thin-film limit. Note
that, by considering model Hamiltonian (6), the physics of the side surfaces,
as discussed in the previous section (e.g., the difference in DP energy) is no
longer captured.

The material parameters were obtained for thicknesses varying between 2 and 6
QL. We imposed the constraint $|D|<|B|$ such that the obtained parameters
yield a gapped dispersion. It has been shown that one can obtain parameters
corresponding to a QSHI or a NI from fitting to the same surface-state
spectrum without any discernible difference [33]. Thus, we take the liberty of
imposing the constraint $\Delta B>0$ when the $GW$ results suggest that the
system is a QSHI and conversely impose $\Delta B<0$ when the system is a NI
[34]. The obtained parameters are shown in Figs. 5a,b,c,d and listed in Table
3. In Fig. 5e the _GW_ band structure of 5 QL thick Sb2Te3 thin-film is
presented, together with the fitted dispersion, and in Fig. 5f the dispersion
of a nanoribbon with the same thickness and width $W=100$ nm is shown. At this
thickness, the Sb2Te3 thin-film is a QSHI, thus edge states appear in the
dispersion of the nanoribbon with energies in the hybridization gap.

##  V Discussion

Comparing the values in Tables 1 and 2, it is clear that there is substantial
variation between all the parameter sets and their properties that are highly
relevant for the low-energy spectrum (e.g., bulk gap, DP positioning and
asymmetry). Hence, it can be of crucial importance to consider an appropriate
parameter set that is tailored to the specific TI material for accurately
modeling the nanoelectronic properties of TI nanostructure-based devices.

In this work, we focus on the accurate low-energy description of 3D TI
nanostructures, including the thin-film limit with a hybridization gap, and
the impact of anisotropy and electron-hole asymmetry. However, there are more
aspects that are not under consideration in this work while being relevant for
nanoelectronic device modeling. Aside from a thickness-dependent hybridization
gap, there is also a thickness-dependent energy shift of the spectrum,
captured by thin-film model parameter $E_{0}$. For this shift, the $GW$
calculations that we considered for our model fits do not line up with
experimental data, while other many-body calculations match experimental
values better [36]. However, such a shift can be easily taken into account by
adjusting the Fermi level in the nanostructure simulations accordingly.
Another important effect can be observed in 3D TI films grown by molecular
beam epitaxy. The top surface is usually exposed to vacuum and the bottom
surface is attached to a substrate, breaking the inversion symmetry along the
$\hat{z}$ direction. In this case, a structure inversion asymmetry term can
account for the influence of the substrate [13, 54]. Such effects can also be
revealed through asymmetric electrostatic gating [10, 55]. While not under
consideration here, we can also introduce such terms in a straightforward
manner in our models.

|  |  $N_{\text{QL}}$
---|---|---
|  |  2 |  3 |  4 |  5 |  6
$\Delta E_{\text{hyb}}~{}[\text{meV}]$ |  Bi2Se3 [3] |  137 |  gray!5040 |  gray!5022 |  - |  \-
|  Bi2Se3 [4] |  191 |  58 |  22 |  10 |  4
|  Bi2Se3 [5] |  133 |  42 |  15 |  5 |  -
|  Bi2Se3 fit |  13 |  gray!5016 |  - |  - |  -
|  Bi2Se3 _GW_ [33] |  238 |  83 |  34 |  (16) |  (7)
|  Bi2Te3 [4] |  - |  - |  - |  - |  -
|  Bi2Te3 [5] |  gray!5091 |  gray!50138 |  54 |  35 |  gray!5025
|  Bi2Te3 fit |  gray!50161 |  50 |  6 |  gray!5018 |  11
|  Bi2Te3 _GW_ [34, 53] |  gray!50 152 |  26 |  (14) |  gray!50 (3) |  gray!50 (0)
|  Sb2Te3 [4] |  536 |  15 |  gray!50 99 |  gray!50 68 |  gray!50 9
|  Sb2Te3 [5] |  747 |  156 |  gray!50 27 |  gray!50 57 |  gray!50 41
|  Sb2Te3 fit |  587 |  78 |  gray!50 65 |  gray!50 73 |  gray!50 30
|  Sb2Te3 _GW_ [34, 53] |  254 |  60 |  gray!50 (4) |  gray!50 (16) |  gray!50 (12)
Table 2: The size of the hybridization gap $\Delta E_{\text{hyb}}$ of Bi2Se3,
Bi2Te3 and Sb2Te3 thin films, described using the bulk model Hamiltonian (1)
for different sets of material parameters and varying number of quintuple
layers $N_{\text{QL}}$. For each case, the topological phase is indicated by a
gray (blank) shading if the thin film is a QSHI (NI). The results of the _GW_
calculations are also given. The gaps in parentheses should not be taken at
face value, as the accuracy of the _GW_ calculations does not allow a safe
statement when the band gap becomes too small [53].  Figure 5: (a),(b),(c),(d)
The parameters of the effective thin-film model (6) for film thicknesses from
2 to 6 QLs determined by fitting to spectra obtained with the _GW_ method [33,
34]. The parameters are obtained by imposing $|D|<|B|$, and $\Delta\cdot
B<(>)\,0$ when the thin film is a NI (QSHI). (e) Band structure of a 5 QL
thick Sb2Te3 thin-film obtained with the _GW_ method (continuous black line),
and the fitted dispersion (red dashed line). The shaded region corresponds to
the hybridization gap $\Delta E_{\text{hyb}}$. (f) The spectrum of a Sb2Te3
nanoribbon with $W=100$ nm and $d=5$ QL. The parameters were taken from Table
3.

In addition to offering a good starting point for nanoelectronic device
modeling, the models considered in this work are also suitable for the study
of hybrid devices that include superconductivity, since an extension to a
Bogoliubov-de Gennes framework follows naturally [56, 57]. The effect of
proximity-induced superconducting pairing in 3D TIs has received considerable
interest [58, 59, 60, 61, 62, 56, 63, 64], with the effective thin-film model
(6) being considered for the study of (proximitized) magnetically doped 3D TI
nanoribbons [30, 29, 27, 65]. In such systems, typically electron-hole
asymmetry is neglected ($D=0$). It has been shown, however, that electron-hole
asymmetry can play an important role in the transport properties of
magnetically doped 3D TI films [54].

Here, we focus on obtaining suitable model parameters for Bi2Se3, Bi2Te3 and
Sb2Te3 by considering four-band models and a parameter fitting strategy
instead of perturbation theory applied at $\Gamma$. Alternatively, however, it
has already been shown that an eight-band model Hamiltonian obtained from
perturbation theory can also give an accurate low-energy description of all
three materials [5]. However, as that approach is also expected to become
unreliable in the thin-film limit and is computationally more demanding, our
approach offers some distinct advantages for efficient and accurate
nanoelectronic device modeling.

##  VI Conclusion

We model the low-energy electronic spectrum of 3D TI nanostructures (e.g.,
nanowires and nanoribbons) based on Bi2Se3, Bi2Te3 and Sb2Te3 in detail. We
use the commonly considered four-band (bulk and thin-film) model Hamiltonians
and derive new parameter sets by fitting to _ab initio_ band structure data.
Our fitting strategy is tailored to accurately capture the (in general,
anisotropic and electron-hole asymmetric) low-energy electronic structure of
the Dirac surface states, while avoiding any unphysical behavior that may
arise when the Hamiltonian is discretized on a lattice. We studied the
accuracy of the obtained fitted material parameters in the thin-film limit,
when the surface states hybridize through the bulk, by using a discretized
version of the bulk Hamiltonian. We have found that our fitting method yields
material parameters that capture the size and topology of the hybridization
gap in Bi2Te3 remarkably well. However, both our new and existing sets of
parameters cannot reliably describe the thin-film limit for all thicknesses
and materials under consideration. Hence, we resort to a thin-film model with
material parameters extracted from the surface-state spectra of thin-film _GW_
calculations. With our new parameter sets, the considered models provide a
suitable framework for simulating the low-energy spectrum and corresponding
properties (e.g., topology, transport) of 3D TI-based nanoelectronic devices
with a broad range of applications.

###### Acknowledgements.

We thank Peter Krüger for providing the band structure data of the _GW_
calculations in the thin-film limit. We also thank Jonas Buchhorn, Malcolm
Connolly, Stefanos Dimitriadis, Dennis Heffels, Jan Karthein, Jalil Abdur
Rehman, Thomas Schäpers, Michael Schleenvoigt, Peter Schüffelgen, Kaycee
Underwood and Max Vaßen-Carl for useful discussions. This work is supported by
the QuantERA grant MAGMA (by the German Research Foundation under grant
491798118, by the National Research Fund Luxembourg under Grant No.
INTER/QUANTERA21/16447820/MAGMA, and by MCIN/AEI/10.13039/501100011033 and the
European Union NextGenerationEU/PRTR under project PCI2022-132927), and by
Germany’s Excellence Strategy – Cluster of Excellence Matter and Light for
Quantum Computing (ML4Q) EXC 2004/1 – 390534769. K.M. and P.R. acknowledge the
financial support by the Bavarian Ministry of Economic Affairs, Regional
Development and Energy within Bavaria’s High-Tech Agenda Project “Bausteine
für das Quantencomputing auf Basis topologischer Materialien mit
experimentellen und theoretischen Ansätzen” (Grant No. 07 02/686 58/1/21 1/22
2/23), and K.M. acknowledges the financial support by the Quantum Future
project ‘MajoranaChips’ (Grant No. 13N15264) within the funding program
Photonic Research Germany.

##  Appendix A _ab initio_ calculations

In our density functional theory (DFT) calculations we use the full-potential
relativistic Korringa-Kohn-Rostoker Green function method (KKR) [66] as
implemented in the JuKKR code [67]. Our calculations are carried out for the
experimental crystal structures for Bi2Te3 [68], Sb2Te3 [69] and Bi2Se3 [70]
and we parametrize the exchange correlation functional using the local density
approximation (LDA) [71] because a comparison of LDA and generalized gradient
approximation (using the PBE functional [72]) resulted in gap sizes that
reproduce the experimentally observed gap better. We employ Lloyd’s formula
[73] to correct for the error arising from the finite $\ell_{\text{max}}=3$
cutoff in the angular momentum expansion of the space filling Voronoi cells
around the atomic centers, where the exact (i.e. full-potential) description
of the atomic shapes is taken into account [74, 75]. The DFT calculations are
orchestrated using the AiiDA-KKR plugin [76] to the AiiDA infrastructure for
automated FAIR data provenance tracking [77]. Our results are uploaded to the
Materials cloud archive [78] and the JuKKR and AiiDA-KKR codes are published
as open-source software [67, 79].

##  Appendix B Fitting procedure

The material parameters for the bulk model (Eq. (1)) and the effective thin-
film model (Eq. (6)) for the three studied materials are obtained by fitting
to _ab initio_ band structures.

Figure 6: (a),(c) The dispersion resulting from the fit (orange) to _ab initio_ band structures (black dots), for (a) Bi2Se3 and (c) Sb2Te3, compared with dispersion evaluated with the material parameters taken from Ref. [4] (green) and Ref. [5] (cyan). The upper limits of the fit in momentum space $k_{\parallel,z}^{\text{max}}$ are shown with gray dashed lines. (b),(d) A map of the $R^{2}$ of the fits for (b) Bi2Se3 and (d) Sb2Te3, obtained for different values of $k_{\parallel}^{\text{max}}$ and $k_{z}^{\text{max}}$. A black cross marks the selected pair of ($k_{\parallel}^{\text{max}}$, $k_{z}^{\text{max}}$) values. | $N_{\text{QL}}$ | $E_{0}~{}[\text{eV}]$ | $D~{}[\text{eV\AA}^{2}]$ | $\Delta~{}[\text{eV}]$ | $B~{}[\text{eV\AA}^{2}]$ | $v_{\text{F}}~{}[\text{eV\AA}]$
---|---|---|---|---|---|---
Bi2Se3 | 2 | 0.121 | -16.14 | 0.239 | -17.56 | -0.048
3 | 0.043 | -13.94 | 0.082 | -15.82 | 1.697
4 | 0.018 | -13.31 | 0.034 | -16.35 | 1.920
5 | 0.008 | -13.23 | 0.014 | -16.51 | 2.010
6 | -0.002 | -13.06 | 0.006 | -16.75 | 2.046
Bi2Te3 | 2 | 0.077 | -28.41 | -0.153 | -29.14 | 2.463
3 | 0.013 | -27.52 | 0.027 | -28.10 | 0.876
4 | 0.001 | -29.55 | 0.003 | -30.06 | 1.301
5 | 0.002 | -28.50 | -0.006 | -28.97 | 1.340
6 | 0.000 | -28.36 | -0.000 | -28.81 | 1.232
Sb2Te3 | 2 | 0.127 | -21.70 | 0.254 | -26.29 | 0.482
3 | 0.030 | -15.45 | 0.064 | -16.04 | 2.921
4 | 0.002 | -13.39 | -0.007 | -13.89 | 2.952
5 | 0.008 | -12.65 | -0.018 | -13.10 | 2.870
6 | 0.006 | -13.38 | -0.012 | -13.83 | 2.887
Table 3: Parameters of the thin-film Hamiltonian (6), for Bi2Se3, Bi2Te3 and
Sb2Te3, for thicknesses ranging 2-6 QL.

In the case of the bulk model, the eigenvalues of the Hamiltonian were fitted
to the relevant bands from $\Gamma\equiv 0$ up to a certain
$k_{\parallel}^{\text{max}}$ value in $k_{x}$ and $k_{y}$, and
$k_{z}^{\text{max}}$ in $k_{z}$, with a sampling interval of
$0.003~{}\text{\AA}^{-1}$ along the three axes (for reference, the distance to
other symmetry points shown in Fig. 1c is $0.882~{}\text{\AA}^{-1}$ for
$\Gamma\textnormal{-}\text{L}$, $0.903~{}\text{\AA}^{-1}$ for
$\Gamma\textnormal{-}\text{F}$, and $0.329~{}\text{\AA}^{-1}$ for
$\Gamma\textnormal{-}\text{Z}$). For each fit corresponding to a pair of
$(k_{\parallel}^{\text{max}},k_{z}^{\text{max}})$ values, the $R^{2}$ of the
fit was evaluated (see Figs. 6b,d). The best fit is obtained by choosing the
one for which
$R^{2}(k_{\parallel}^{\text{max}},k_{z}^{\text{max}})>0.999\cdot\text{max}(R^{2})$
and also $k_{\parallel}^{\text{max}}+k_{z}^{\text{max}}$ is maximal, while
imposing certain constraints on $k_{z,\parallel}^{\text{max}}$, as explained
below. In Figs. 6a,c, we show the result of our fits for Bi2Se3 and Sb2Te3 to
the relevant bands, together with the dispersion obtained with parameters
taken from Refs. [4, 5]. In Figs. 6b,d, maps of the $R^{2}$ values obtained
for different pairs of $k_{\parallel}^{\text{max}}$ and $k_{z}^{\text{max}}$
are shown, with the optimal value indicated with a black cross. For the three
materials the maximal value of $R^{2}$ was 0.9891, 0.9868 and 0.9458, for
Bi2Se3, Bi2Te3 and Sb2Te3, respectively.

We impose a lower bound of
$k_{z,\parallel}^{\text{max}}>0.045~{}\text{\AA}^{-1}$ such that we capture at
least the same extent as can be accurately captured with a perturbative
approach [52]. In order to have a fit that covers the whole region in
$k$-space where the surface states can be found inside the bulk gap, we impose
an additional constraint. The points in momentum space where the Dirac cone of
the surface states joins the bulk-band energies can be approximated by

$\displaystyle
k_{\parallel}^{\text{surf}}=\frac{\text{max}_{i\in\\{0,1\\}}\\{E_{i}+(-1)^{i}(C_{0}-C_{1}M_{0}/M_{1})\\}}{A_{0}\sqrt{1-(C_{1}/M_{1})^{2}}},$

for the in-plane dispersion,
$E_{0}=-\text{max}_{\mathbf{k}}\\{E_{\text{VB}}(\mathbf{k})\\}$ and
$E_{1}=\text{min}_{\mathbf{k}}\\{E_{\text{CB}}(\mathbf{k})\\}$. For
completeness, the dispersion of the surface states with other orientations
also have to be considered, e.g., orthogonal to the $\hat{y}$ direction. In
this case, the Dirac cone will join the bulk-band energies at

$\displaystyle
k_{\parallel}^{\text{surf}}=\frac{\text{max}_{i\in\\{0,1\\}}\\{E_{i}+(-1)^{i}(C_{0}-C_{2}M_{0}/M_{2})\\}}{A_{0}\sqrt{1-(C_{2}/M_{2})^{2}}},$
$\displaystyle
k_{z}^{\text{surf}}=\frac{\text{max}_{i\in\\{0,1\\}}\\{E_{i}+(-1)^{i}(C_{0}-C_{2}M_{0}/M_{2})\\}}{B_{0}\sqrt{1-(C_{2}/M_{2})^{2}}},$

for the dispersion in $k_{x}$ and $k_{z}$, respectively. Using values found in
literature, we found that
$k_{\parallel}^{\text{surf}},k_{z}^{\text{surf}}<0.2~{}\text{\AA}^{-1}$, which
is contained in the region of $k$-space over which we sample the _ab initio_
band structures. Hence, we can capture the surface states inside the bulk gap
by imposing the conditions
$k_{\parallel}^{\text{max}}>k_{\parallel}^{\text{surf}}$ and
$k_{z}^{\text{max}}>k_{z}^{\text{surf}}$, where
$k_{\parallel,z}^{\text{surf}}$ are evaluated using the _ab initio_ data and
the fitted parameters.

For the effective surface-state model, a similar procedure was used, with the
difference that the fit was performed for $k_{x}$ and $k_{y}$ up to
$k_{\parallel}^{\text{max}}$, and the fit with
$R^{2}>0.999\cdot\text{max}(R^{2})$ and maximal $k_{\parallel}^{\text{max}}$
was selected.

## References

  * Hasan and Moore [2011] M. Z. Hasan and J. E. Moore, Annu. Rev. Condens. Matter Phys. 2, 55 (2011).
  * Breunig and Ando [2022] O. Breunig and Y. Ando, Nat. Rev. Phys. 4, 184 (2022).
  * Zhang  _et al._ [2009] H. Zhang, C.-X. Liu, X.-L. Qi, X. Dai, Z. Fang, and S.-C. Zhang, Nat. Phys. 5, 438 (2009).
  * Liu  _et al._ [2010a] C.-X. Liu, X.-L. Qi, H. Zhang, X. Dai, Z. Fang, and S.-C. Zhang, Phys. Rev. B 82, 045122 (2010a).
  * Nechaev and Krasovskii [2016] I. A. Nechaev and E. E. Krasovskii, Phys. Rev. B 94, 201410 (2016).
  * Bardarson  _et al._ [2010] J. H. Bardarson, P. W. Brouwer, and J. E. Moore, Phys. Rev. Lett. 105, 156803 (2010).
  * Brey and Fertig [2014] L. Brey and H. A. Fertig, Phys. Rev. B 89, 085305 (2014).
  * Bardarson and Ilan [2018] J. H. Bardarson and R. Ilan, Transport in topological insulator nanowires, in  _Springer Series in Solid-State Sciences_ (Springer International Publishing, 2018) p. 93–114.
  * Moors  _et al._ [2018] K. Moors, P. Schüffelgen, D. Rosenbach, T. Schmitt, T. Schäpers, and T. L. Schmidt, Phys. Rev. B 97, 245429 (2018).
  * Ziegler  _et al._ [2018] J. Ziegler, R. Kozlovsky, C. Gorini, M.-H. Liu, S. Weishäupl, H. Maier, R. Fischer, D. A. Kozlov, Z. D. Kvon, N. Mikhailov, S. A. Dvoretsky, K. Richter, and D. Weiss, Phys. Rev. B 97, 035157 (2018).
  * Linder  _et al._ [2009] J. Linder, T. Yokoyama, and A. Sudbø, Phys. Rev. B 80, 205401 (2009).
  * Liu  _et al._ [2010b] C.-X. Liu, H. Zhang, B. Yan, X.-L. Qi, T. Frauenheim, X. Dai, Z. Fang, and S.-C. Zhang, Phys. Rev. B 81, 041307(R) (2010b).
  * Shan  _et al._ [2010] W.-Y. Shan, H.-Z. Lu, and S.-Q. Shen, New J. Phys. 12, 043048 (2010).
  * Lu  _et al._ [2010] H.-Z. Lu, W.-Y. Shan, W. Yao, Q. Niu, and S.-Q. Shen, Phys. Rev. B 81, 115407 (2010).
  * Zhang  _et al._ [2015] S.-B. Zhang, H.-Z. Lu, and S.-Q. Shen, Sci. Rep. 5, 13277 (2015).
  * Leis  _et al._ [2022] A. Leis, M. Schleenvoigt, K. Moors, H. Soltner, V. Cherepanov, P. Schüffelgen, G. Mussler, D. Grützmacher, B. Voigtländer, F. Lüpke, and F. S. Tautz, Adv. Quantum Technol. 5, 2200043 (2022).
  * Yu  _et al._ [2010] R. Yu, W. Zhang, H.-J. Zhang, S.-C. Zhang, X. Dai, and Z. Fang, Science 329, 61 (2010).
  * Chang  _et al._ [2013] C.-Z. Chang, J. Zhang, X. Feng, J. Shen, Z. Zhang, M. Guo, K. Li, Y. Ou, P. Wei, L.-L. Wang, Z.-Q. Ji, Y. Feng, S. Ji, X. Chen, J. Jia, X. Dai, Z. Fang, S.-C. Zhang, K. He, Y. Wang, L. Lu, X.-C. Ma, and Q.-K. Xue, Science 340, 167 (2013).
  * Kou  _et al._ [2013] X. Kou, M. Lang, Y. Fan, Y. Jiang, T. Nie, J. Zhang, W. Jiang, Y. Wang, Y. Yao, L. He, and K. L. Wang, ACS Nano 7, 9205 (2013).
  * Wang  _et al._ [2015a] J. Wang, B. Lian, and S.-C. Zhang, Phys. Scripta 2015, 014003 (2015a).
  * Yasuda  _et al._ [2017] K. Yasuda, M. Mogi, R. Yoshimi, A. Tsukazaki, K. S. Takahashi, M. Kawasaki, F. Kagawa, and Y. Tokura, Science 358, 1311 (2017).
  * Sun  _et al._ [2019] H. Sun, B. Xia, Z. Chen, Y. Zhang, P. Liu, Q. Yao, H. Tang, Y. Zhao, H. Xu, and Q. Liu, Phys. Rev. Lett. 123, 096401 (2019).
  * Tokura  _et al._ [2019] Y. Tokura, K. Yasuda, and A. Tsukazaki, Nat. Rev. Phys. 1, 126 (2019).
  * Yasuda  _et al._ [2020] K. Yasuda, T. Morimoto, R. Yoshimi, M. Mogi, A. Tsukazaki, M. Kawamura, K. Takahashi, M. Kawasaki, N. Nagaosa, and Y. Tokura, Nat. Nanotechnol. 15, 831–835 (2020).
  * Wang  _et al._ [2021] P. Wang, J. Ge, J. Li, Y. Liu, Y. Xu, and J. Wang, Innovation 2, 100098 (2021).
  * Qiu  _et al._ [2022] G. Qiu, P. Zhang, P. Deng, S. K. Chong, L. Tai, C. Eckberg, and K. L. Wang, Phys. Rev. Lett. 128, 217704 (2022).
  * Atanov  _et al._ [2024] O. Atanov, W. T. Tai, Y.-M. Xie, Y. H. Ng, M. A. Hammond, T. S. Manfred Ho, T. H. Koo, H. Li, S. L. Ho, J. Lyu, S. Chong, P. Zhang, L. Tai, J. Wang, K. T. Law, K. L. Wang, and R. Lortz, Cell Rep. Phys. Sci. 5, 101762 (2024).
  * Wang  _et al._ [2015b] J. Wang, Q. Zhou, B. Lian, and S.-C. Zhang, Phys. Rev. B 92, 064520 (2015b).
  * Zeng  _et al._ [2018] Y. Zeng, C. Lei, G. Chaudhary, and A. H. MacDonald, Phys. Rev. B 97, 081102(R) (2018).
  * Chen  _et al._ [2018] C.-Z. Chen, Y.-M. Xie, J. Liu, P. A. Lee, and K. T. Law, Phys. Rev. B 97, 104504 (2018).
  * Mandal  _et al._ [2022] P. Mandal, N. Taufertshöfer, L. Lunczer, M. P. Stehno, C. Gould, and L. W. Molenkamp, Nano Lett. 22, 3557 (2022).
  * Uday  _et al._ [2023] A. Uday, G. Lippertz, K. Moors, H. F. Legg, A. Bliesener, L. M. C. Pereira, A. A. Taskin, and Y. Ando 10.48550/arXiv.2307.08578 (2023), arXiv:2307.08578 [cond-mat.mes-hall] .
  * Förster  _et al._ [2015] T. Förster, P. Krüger, and M. Rohlfing, Phys. Rev. B 92, 201404 (2015).
  * Förster  _et al._ [2016] T. Förster, P. Krüger, and M. Rohlfing, Phys. Rev. B 93, 205442 (2016).
  * Wang  _et al._ [2019] Z. Wang, T. Zhou, T. Jiang, H. Sun, Y. Zang, Y. Gong, J. Zhang, M. Tong, X. Xie, Q. Liu, C. Chen, K. He, and Q.-K. Xue, Nano Lett. 19, 4627 (2019).
  * Liu  _et al._ [2023] J.-n. Liu, X. Yang, H. Xue, X.-s. Gai, R. Sun, Y. Li, Z.-z. Gong, N. Li, Z.-K. Xie, W. He, X.-Q. Zhang, D. Xue, and Z.-H. Cheng, Nat. Commun. 14, 4424 (2023).
  * Shen [2012] S.-Q. Shen,  _Topological Insulators: Dirac Equation in Condensed Matters_ (Springer Berlin Heidelberg, 2012).
  * Nielsen and Ninomiya [1981] H. Nielsen and M. Ninomiya, Nucl. Phys. B 185, 20 (1981).
  * Datta [2005] S. Datta,  _Quantum Transport: Atom to Transistor_ (Cambridge University Press, 2005).
  * Zhang  _et al._ [2010] Y. Zhang, K. He, C.-Z. Chang, C.-L. Song, L.-L. Wang, X. Chen, J.-F. Jia, Z. Fang, X. Dai, W.-Y. Shan, S.-Q. Shen, Q. Niu, X.-L. Qi, S.-C. Zhang, X.-C. Ma, and Q.-K. Xue, Nat. Phys. 6, 584 (2010).
  * Liu  _et al._ [2012] Y. Liu, G. Bian, T. Miller, M. Bissen, and T.-C. Chiang, Phys. Rev. B 85, 195442 (2012).
  * Jiang  _et al._ [2012] Y. Jiang, Y. Wang, M. Chen, Z. Li, C. Song, K. He, L. Wang, X. Chen, X. Ma, and Q.-K. Xue, Phys. Rev. Lett. 108, 016401 (2012).
  * Bernevig  _et al._ [2006] B. A. Bernevig, T. L. Hughes, and S.-C. Zhang, Science 314, 1757 (2006).
  * Groth  _et al._ [2014] C. W. Groth, M. Wimmer, A. R. Akhmerov, and X. Waintal, New J. Phys. 16, 063065 (2014).
  * Nijholt  _et al._ [2019] B. Nijholt, J. Weston, J. Hoofwijk, and A. Akhmerov, Adaptive: parallel active learning of mathematical functions (2019).
  * Zhang  _et al._ [2011] J. Zhang, C.-Z. Chang, Z. Zhang, J. Wen, X. Feng, K. Li, M. Liu, K. He, L. Wang, X. Chen, Q.-K. Xue, X. Ma, and Y. Wang, Nat. Commun. 2, 574 (2011).
  * Vidal  _et al._ [2013] F. Vidal, M. Eddrief, B. Rache Salles, I. Vobornik, E. Velez-Fort, G. Panaccione, and M. Marangolo, Phys. Rev. B 88, 241410 (2013).
  * Landolt  _et al._ [2014] G. Landolt, S. Schreyeck, S. V. Eremeev, B. Slomski, S. Muff, J. Osterwalder, E. V. Chulkov, C. Gould, G. Karczewski, K. Brunner, H. Buhmann, L. W. Molenkamp, and J. H. Dil, Phys. Rev. Lett. 112, 057601 (2014).
  * Neupane  _et al._ [2014] M. Neupane, A. Richardella, J. Sánchez-Barriga, S. Xu, N. Alidoust, I. Belopolski, C. Liu, G. Bian, D. Zhang, D. Marchenko, A. Varykhalov, O. Rader, M. Leandersson, T. Balasubramanian, T.-R. Chang, H.-T. Jeng, S. Basak, H. Lin, A. Bansil, N. Samarth, and M. Z. Hasan, Nat. Commun. 5, 3841 (2014).
  * Kong  _et al._ [2010] D. Kong, J. C. Randel, H. Peng, J. J. Cha, S. Meister, K. Lai, Y. Chen, Z.-X. Shen, H. C. Manoharan, and Y. Cui, Nano Lett. 10, 329 (2010).
  * Xiu  _et al._ [2011] F. Xiu, L. He, Y. Wang, L. Cheng, L.-T. Chang, M. Lang, G. Huang, X. Kou, Y. Zhou, X. Jiang, Z. Chen, J. Zou, A. Shailos, and K. L. Wang, Nat. Nanotechnol. 6, 216–221 (2011).
  * Zhang and Vishwanath [2010] Y. Zhang and A. Vishwanath, Phys. Rev. Lett. 105, 206601 (2010).
  * Förster [2016] T. Förster,  _Topological insulators ${Bi}_{2}{Se}_{3}$, ${Bi}_{2}{Te}_{3}$, and ${Sb}_{2}{Te}_{3}$: Electronic and topological properties of surfaces and thin films from DFT and GW calculations_, Phd thesis, University of Münster (2016).
  * Lu  _et al._ [2013] H.-Z. Lu, A. Zhao, and S.-Q. Shen, Phys. Rev. Lett. 111, 146802 (2013).
  * Rosenbach  _et al._ [2022] D. Rosenbach, K. Moors, A. R. Jalil, J. Kölzer, E. Zimmermann, J. Schubert, S. Karimzadah, G. Mussler, P. Schüffelgen, D. Grützmacher, H. Lüth, and T. Schäpers, SciPost Phys. Core 5, 017 (2022).
  * Sitthison and Stanescu [2014] P. Sitthison and T. D. Stanescu, Phys. Rev. B 90, 035313 (2014).
  * Heffels  _et al._ [2023] D. Heffels, D. Burke, M. R. Connolly, P. Schüffelgen, D. Grützmacher, and K. Moors, Nanomaterials 13, 723 (2023).
  * Fu and Kane [2008] L. Fu and C. L. Kane, Phys. Rev. Lett. 100, 096407 (2008).
  * Stanescu  _et al._ [2010] T. D. Stanescu, J. D. Sau, R. M. Lutchyn, and S. Das Sarma, Phys. Rev. B 81, 241310 (2010).
  * Potter and Lee [2011] A. C. Potter and P. A. Lee, Phys. Rev. B 83, 184520 (2011).
  * Cook and Franz [2011] A. Cook and M. Franz, Phys. Rev. B 84, 201105 (2011).
  * Cook  _et al._ [2012] A. M. Cook, M. M. Vazifeh, and M. Franz, Phys. Rev. B 86, 155431 (2012).
  * Chiu  _et al._ [2016] C.-K. Chiu, W. S. Cole, and S. Das Sarma, Phys. Rev. B 94, 125304 (2016).
  * Bai  _et al._ [2022] M. Bai, X.-K. Wei, J. Feng, M. Luysberg, A. Bliesener, G. Lippertz, A. Uday, A. A. Taskin, J. Mayer, and Y. Ando, Commun. Mater. 3, 20 (2022).
  * Burke  _et al._ [2024] D. Burke, D. Heffels, K. Moors, P. Schüffelgen, D. Grützmacher, and M. R. Connolly, Phys. Rev. B 109, 045138 (2024).
  * Ebert  _et al._ [2011] H. Ebert, D. Ködderitzsch, and J. Minár, Rep. Prog. Phys. 74, 096501 (2011).
  * Rüßmann  _et al._ [2022] P. Rüßmann, D. A. Silva, D. Bauer, P. Baumeister, M. Bornemann, J. Bouaziz, S. Brinker, J. Chico, P. Dederichs, B. Drittler, F. D. Santos, M. dos Santos Dias, N. Essing, G. Géranton, I. Klepetsanis, A. Kosma, N. Long, S. Lounis, P. Mavropoulos, E. M. Tapia, C. Oran, N. Papanikolaou, E. Rabel, B. Schweflinghaus, N. Stefanou, A. Thiess, R. Zeller, B. Zimmermann, and S. Blügel, Zenodo 10.5281/zenodo.7284738 (2022).
  * Nakajima [1963] S. Nakajima, J. Phys. Chem. Solids , 479 (1963).
  * Ullner [1968] H.-A. Ullner, Annalen der Physik 21, 45 (1968).
  * Wyckoff [1964] R. W. G. Wyckoff,  _Crystal Structures (2 ed.)_ (J. Wiley and Sons, 1964).
  * Vosko  _et al._ [1980] S. H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. 58, 1200 (1980).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Zeller [2004] R. Zeller, J. Phys.: Condens. Matter 16, 6453 (2004).
  * Stefanou  _et al._ [1990] N. Stefanou, H. Akai, and R. Zeller, Comput. Phys. Commun. 60, 231 (1990).
  * Stefanou and Zeller [1991] N. Stefanou and R. Zeller, J. Phys.: Cond. Matter 3, 7599 (1991).
  * Rüßmann  _et al._ [2021] P. Rüßmann, F. Bertoldo, and S. Blügel, npj Comput. Mater. 7, 13 (2021).
  * Huber  _et al._ [2020] S. P. Huber, S. Zoupanos, M. Uhrin, L. Talirz, L. Kahle, R. Häuselmann, D. Gresch, T. Müller, A. V. Yakutovich, C. W. Andersen, F. F. Ramirez, C. S. Adorf, F. Gargiulo, S. Kumbhar, E. Passaro, C. Johnston, A. Merkys, A. Cepellotti, N. Mounet, N. Marzari, B. Kozinsky, and G. Pizzi, Sci. Data 7, 300 (2020).
  * Zsurka  _et al._ [2024] E. Zsurka, C. Wang, J. Legendre, D. D. Miceli, L. Serra, D. Grützmacher, T. L. Schmidt, P. Rüßmann, and K. Moors, Materials Cloud Archive 2024.X, 10.24435/materialscloud:mx-bn (2024).
  * Rüßmann  _et al._ [2023] P. Rüßmann, D. Antognini Silva, R. Aliberti, J. Bröder, H. Janssen, R. Mozumder, M. Struckmann, J. Wasmer, and S. Blügel, Zenodo 10.5281/zenodo.3628251 (2023).