Quantum confinement of the Dirac surface states in topological-insulator nanowires

The non-trivial topology of three-dimensional topological insulators dictates the appearance of gapless Dirac surface states. Intriguingly, when made into a nanowire, quantum confinement leads to a peculiar gapped Dirac sub-band structure. This gap is useful for, e.g., future Majorana qubits based on TIs. Furthermore, these sub-bands can be manipulated by a magnetic flux and are an ideal platform for generating stable Majorana zero modes, playing a key role in topological quantum computing. However, direct evidence for the Dirac sub-bands in TI nanowires has not been reported so far. Here, using devices fabricated from thin bulk-insulating (Bi1−xSbx)2Te3 nanowires we show that non-equidistant resistance peaks, observed upon gate-tuning the chemical potential across the Dirac point, are the unique signatures of the quantized sub-bands. These TI nanowires open the way to address the topological mesoscopic physics, and eventually the Majorana physics when proximitized by an s-wave superconductor.


Supplementary Note 3. Universal conductance fluctuations
In mesoscopic systems universal conductance fluctuations (UCF) arise from a specific distribution of scattering sites. Our devices fall into this regime, which raises the question whether type I peaks could also be explained by UCF. However, type I peaks are much too regular to be caused by conductance fluctuations. UCF occur without any clear periodicity as a function of gate voltage or magnetic field. In contrast, the features of type I that we identify with the sub-band crossings occur in a regular fashion, i.e. at regularly spaced gate voltages, as shown by Fig. 2d. Secondly, the amplitude of type I features strongly differs from what would be expected for typical UCF. For UCF, the amplitude of individual peaks is random. In contrast, our type I oscillations occur with a largely uniform amplitude. Both of these facts speak against UCF as the origin of type I peaks. Furthermore, for every type I peak at gate voltage ∆V G we observe a symmetry-related peak close to −∆V G (see Fig. 2d). This is expected from the presence of an approximate particle-hole symmetry of the Dirac cone of (Bi 1x Sb x ) 2 Te 3 [4] and is not consistent with UCFs for which no such symmetry is expected.
Another hallmark of UCFs is that they are strongly affected by small magnetic fields. If features of type I in our gating curves were UCF, they should rapidly change when applying a small magnetic field. Using device 6 (fabricated in the same way to devices 1-5) multiple gate sweeps were performed at constant magnetic fields of 0.1 and 0.5 T. Fig. 2 shows the average of these sweeps (see Supplementary Note 4 below). While there is finite differences in the background (which partly consists of residual type II fluctuations that may well be field-dependent) the main features that we identify with the sub-band crossings remain robust upon changing the magnetic field. Thus, we can clearly rule out UCFs as the origin of type I peaks, both due to the regularity of their amplitude and periodicity, as well as their independence on small magnetic fields. As discussed in the main text, two distinct types of oscillations are seen in the R(V G ) curves: Large oscillations coming from sub-band crossings (type I, ∼ 3 kΩ) and smaller, time-dependent oscillations (type II, ∼ 0.5 kΩ). The type-I oscillations are reproducible and most of the peaks are well pronounced. On the other hand, as shown examplarily for device 5 in Fig. 3, the resistance of the devices changes with a timescale of a few minutes. The magnitude of these jumps is of the order of the type II features in the R(V G ) curves. We speculate that each jump of the resistance is associated with a reconfiguration of the scattering site distribution and thus changes the background on which we expect to observe the intrinsic gate voltage or magnetic field dependence of the surface state transport. Since the gating-curves were measured in about 80 minutes, this process also takes place during our gate sweeps. The fact that type II features consequently differ between each individual sweep opens a way to suppress their conductivity contribution by performing a time average or, more practically, an average of multiple sweeps. Comparing the individual gate sweeps with the average (see Fig. 2a-c of the main text) clearly shows that large features (type I) are preserved upon averaging while minor peaks (type II) in the individual curves are suppressed. As such, averaging of multiple sweeps is a useful tool to deterministically extract the dominant conductivity features due to the Dirac sub-band structure. Even without performing gate voltage sweep averages, type I features are also found in each individual sweep as discussed for two additionally fabricated devices in the following. For these devices we did not perform averaging at the time of the data acquisition. Yet, they also show full gate-tunability through the Dirac point and their individual gating curves exhibit features similar to those of devices 1-3 discussed in the main text ( Supplementary Fig. 4). By comparison to the theory, major peaks can be indexed for these devices as well, even without performing a time-consuming averaging over a large number of gate sweeps. However, the type-II oscillations can make the identification of the type-I peaks complicated, and in such rare cases it is useful to adopt two criteria for the identification of the type-I peaks: amplitude and regularity. Namely, we identify peaks that stand out in terms of the amplitude and appear semi-regularly to be type-I. For example, the identification of type-I peaks in device 5 ( Supplementary Fig. 4b) followed these criteria. Nevertheless, we emphasize that in most cases (such as device 1 shown in Fig. 2a) the identification of type-I peaks is evident. The intrinsic type-I peaks identified this way remain robust after averaging multiple gate sweeps as mentioned above and turn out to agree well with the theoretical calculations as shown in the main text. In comparison to theory, sub-band crossing peaks can be identified even though type II features are also present in these non-averaged, single gating curves.

Supplementary Note 6. Geometric and quantum capacitances
In general, the gate voltage required to provide the nanowire with additional charge Q can be expressed as where C G is the geometric capacitance of the device and C Q is the quantum capacitance, or, equivalently, as eV G = e 2 C G n(µ) − µ, where e < 0 is the electron charge, µ the chemical potential, and n(µ) the electron density. In our experiments, a change of the gate voltage by 30 V leads to changes of the chemical potential by less then 100 meV. Therefore, the effect of the geometric capacitance dominates by more than two orders of magnitude and we can safely approximate V G ≈ 1 C G en(µ). The geometric capacitance C G was found to be 2 -3.6 pF/m in our experiments. Note that C G depends on materials properties and the device geometry, and one needs to numerically calculate it for each situation [5]. For the present gate geometry, as presented in the section Supplementary Note 13, a value of C G ≈ 15 pF/m would be expected. However, our nanowires are not metallic as is typically assumed in simulations, and hence the result of such numerical simulations should be taken as an upper bound for C G . In addition, the effective capacitance is further affected by parasitic capacitances due to the leads, as well as by the finite length of the nanowire. For these reasons we treat the value of the capacitance as a phenomenological parameter.

Supplementary Note 7. Hamiltonian of a disordered TI nanowire
The Hamiltonian for the surface of a cylindrical TI nanowire of radius R w depends only on momentum around the wire,p φ =L φ /R w = −iv F ∂ φ /R w , and momentum along the wire,p x = −i ∂ x , and is given by [6] This can be simplified by applying a spinor-rotation by φ about the x-axis, i.e., the unitary transformation U x (φ) = exp(−iφσ x /2). The result is the Hamiltonian In the presence of a flux Φ threaded along the wire, Importantly, since spinor rotations are 4π periodic, the eigenfunctions of this Hamiltonian must obey anti-periodic boundary conditions, resulting from U x (φ) = −U x (φ + 2π). As such the wave function can be written in the general form ψ (k) = e i(kx+ φ) ξ, where the anti-periodic boundary conditions require the total angular momentum to satisfy = ±1/2, ±3/2, . . . and ξ is a spinor encoding the spin-structure of the surface state within this rotated system.
Due to the finite geometry along the circumference of the wire the eigenenergies of the Hamiltonian are quantised The band edge of each band is for η = 0 given by ε , Disorder is modelled by impurity potentials located at random positions where r and r i are 2d coordinates on the surface of the wire.

Supplementary Note 8. Density of states and charge density
The free retarded Green's function -which is a matrix -has the form where we take the matrix inverse. Here δ acts as some overall background scattering rate due degrees of freedom (e.g. impurity bands) not described by our Hamiltonian. It has the effect of broadening the divergence of the density of states at the bottom of each sub-band.
The local Green's function of each resolved sub-band is given by From the Green's function we can obtain the local density of states (at a fixed position at the surface of the wire) which, as should be expected, tracks the 2d Dirac density of states ρ 2d (µ) = |µ|/2πv 2 F 2 .
From the density of states we can also obtain the charge density along the wire. It is given by Here n(µ) is zero at the Dirac point and counts the total charge density per length of wire. For µ > 0 this means counting the contribution of all bands above the Dirac point up to the chemical µ. Correspondingly, for µ < 0 counting the charge of all holes contributed by bands below the Dirac point, in other words n(µ) is negative for µ < 0. It is shown as a grey line in Fig. 2d of the main text.

Supplementary Note 9. T-matrix approximation
The full disorder averaged Green's function is where Σ(µ) is the full self-energy from all scattering processes, including disorder and the other scattering processes modelled by the factor iδ (see Supplementary Note 8).
To approximate the self-energy contribution due to disorder scattering we use the full T-matrix calculated at first order in impurity density n imp (i.e. the non-crossing approximation). Within this approximation the self-energy, which is also a matrix, is [7] [Σ] αβ (ω) = n imp u 0 For small δ the ∆(µ, ), from Eq. (6), only has a substantial imaginary part for |µ| > |ε |, this means that a conductivity channel only contributes when the corresponding sub-band has an occupation of electrons (holes) for the upper (lower) Dirac cone. Further the off-diagonal components of Σ are exactly zero for diagonal impurities since the off-diagonal contributions of positive and negative cancel.
The self-energy can also be calculated self-consistently by replacing µ → µ+Σ 11 and ε → ε +Σ 12 , where Σ 11 and Σ 12 are the diagonal and off-diagonal components, respectively. Self-consistency has a similar impact on the self-energy as δ, that is, the breadth of the peaks in scattering rate become broadened in proportion to the self-energy itself. Ultimately this means that peaks at high sub-band index become "washed out" when the self-energy is of the same order of magnitude as the sub-band spacing. This is not of relevance for the experimental situation discussed here where the mean free path is several times the wire radius, but would limit the number of peaks seen in thicker or dirtier wires where the ratio of mean free path and radius can be substantially smaller.

Supplementary Note 10. Conductivity calculation
The conductivity is found using the Kubo formula to calculate the current-current correlation function in linear response. The current operator along the wire is given by j x = e∂H/∂k = −i v F eσ y . We neglect vertex corrections which previous studies have shown only have a quantitative effect on overall conductivity [7]. Within this scheme the DC conductivity is given by where in general the contribution σ is given by where Γ = ImΣ is the (spin-resolved) scattering rate,μ = µ + ReΣ 11 ,ε = ε + ReΣ 12 , and the factors ξ = −(µ + iΓ 11 ) 2 + (ε + iΓ 12 ) 2 and λ = −ε 2 − iε Γ 12 + iμΓ 11 +μ 2 ensure the channel only contributes to conductivity when |µ| |ε |. The approximation in Eq. (12) is valid for small self-energies and diagonal impurities. v F Rw ∼ 1, the curve is highly asymmetric. For u 0 v F Rw 1, the particle-hole symmetry is restored but instead of peaks one obtains dips. We find that only weak u 0 fit the experimental data. This is expected because the wavefunction of the conduction channel wraps around the circumference of the wire, effectivly diluting the impurity by a factor 1/(2πR w ). To be precise: Weak scattering can be defined by u 0ρ 1, wherē ρ ∼ 1 2πv F Rw is the typical density of states in the system for small . Assuming that the width of the scattering potential is of the order of the lattice constant a, the weak scattering limit is reached when the amplitude, V 0 = u 0 /a 2 , of the scattering potential fulfils V 0 v F a 2πRw a . Since in our experiments the circumference is many lattice constants, Rw a ∼ 100, the weak scattering limit is generically realised. Therefore, for the plots in Fig. 2 and Fig. 3 of the main text we use u 0 = 0.1v F R w , n imp = 20/R 2 w and δ = 0.025v F /R w , see also method section. This reproduces the order of magnitude ∼ 5 nm/kΩ of the fluctuations in wire conductivity found experimentally. While the main contribution to the peaks at the sub-band edges arised from the diverging density v F Rw = 20, n imp = 0.1/R 2 w . Strong scattering is capable of turning resistivity peaks at sub-band edges into resistivity dips, however such a situation is unlikely (see text). of states, this effect is further enhance by a matrix-element effect closely related to the topological protection of TI surface states. Here the spin orientation is locked to the propagation direction and scattering from k to −k is prohibited by time-reversal symmetry. In our nanowires, k y = /R w is quantized, and the 1D bands are doubly degenerate as for each k x = k, two transverse momenta, ±k y , are possible. For k x k y , however, the 2D vectors k = (k x , ±k y ) and k = (−k x , ±k y ) are almost antiparallel, leading to a suppression of the scattering rate by the factor (k y /k x ) 2 . This matrix element effect (encoded in the 2 × 2 matrix structure of the T-matrix) substantially suppresses backscattering among occupied channels with /R w k x relative to the scattering to a newly opened channel with a large k y .

Supplementary Note 12. Electron-hole puddles
Several physical effects have not been taken into account by our theoretical treatment. For example, we do not reproduce the type II fluctuations which likely arise from extra charge traps in the experimental setup. We also assume scattering from local impurities, but there will be also long-range potentials arising from randomly distributed positive and negative charges in our compensation doped samples of (Bi 1−x Sb x ) 2 Te 3 . In bulk systems, these lead to the formation of electron-hole puddles [8,9] both in the bulk and in the surface states [10] of TIs. Estimates from Ref. 11 suggest that bulk puddles are absent in TI nanowires with a radius of less than ∼ 100 nm. But also in this case, one can expect substantial fluctuations of the chemical potential on the surface of the TI. For surface states of bulk samples of BiSbTeSe 2 , fluctuations of the local chemical potential with an amplitude of 20 meV on a length scale of 50 nm have been measured by scanning tunneling spectroscopy [10]. While the size of these fluctuations will be reduced for nanowires, surface puddles can be expected to be of quantitative importance (perhaps, also for type II fluctuations).

Supplementary Note 13. Electrostatic modeling and anisotropy effects
In Supplementary Notes 7-11 above, we investigated disorder effects assuming an idealized setup with cylindrical symmetry. Now we investigate the role of anisotropy effects arising both from the hexagonal shape of the wire and of the electrostatic environment. As shown in Supplementary  Fig. 6a, wires with diameters between 29 and 43 nm are positioned on a 255 nm thick dielectric SiO 2 substrate, see methods for details. When in our experiment the chemical potential is adjusted via a flat bottom gate below the dielectric substrate, the induced charge around the wire will not be uniform and will depend on both the substrate/gate geometry and also the approximately hexagonal wire shape. These geometric effects can introduce changes to both the sub-bands and the charging physics of the wires. We will show that Klein-tunnelling physics means that the locations of the sub-band minima in terms of chemical potential are largely unchanged and so the underlying theoretical signature of confinement: An increased scattering at the bottom of each sub-band and a corresponding resistivity peak, should remain robust for most set-ups. In contrast, charging effects could be important such that the mapping from induced charge density to chemical potential is influenced by the geometry of the system. For our experiment, however, we will argue that these effects are below the resolution of the experiment and so the cylindrical wire description above remains a valid approximate model of our experimental set-up. Such effects could, however, be important in future experiments.
We would like to stress that, as discussed above, the following electrostatic modeling is not fully quantitative as it ignores, for example, parasitic capacitances as arising, e.g., from the gold contacts, see Fig. 1d of the main text, and the finite length of the wire. Furthermore, our modeling assumes metallic surfaces of the topological insulator and there are effects arising from the finite penetration depth of the surface state, the experimental uncertainty in the radius of the wire and its surface chemistry, and also possible effects arising from electron-hole puddles, see Supplementary Note 12, which may act as extra charge reservoirs.
Supplementary Note 14. Influence of geometry on energies of sub-band minima We start by addressing the influence of the geometric effects on the sub-band minima and associated resistivity maxima. An important question for our experiment and physical interpretation is the consequence the induced charge anisotropy has for the confinement of surface states around the wire. It turns out that the Dirac nature of the TI surface states leads to Klein tunnelling of the state through any smooth surface potential and that, as a result, the energies of sub-band minima and associated peaks in resistivity are only very weakly sensitive to any charge anisotropy.
The charge anisotropy induced by bottom gating can be modelled as an angular dependent chemical potential µ(φ) =μ + δµ(φ), whereμ is the averaged chemical potential around the wire and δµ(φ) the angular variation of the chemical potential. (An estimate of µ(φ) calculated using an electrostatic model and the full influence on sub-bands, which is important for the charge density, can be found below).
Of relevance to our study is the impact this angular dependent potential has on the quantisation of the surface states. In the cylindrical case the location of the quantisation peaks is set by the energies of the bottom of the sub-bands i.e. at k = 0. Adding an inhomogeneous chemical potential µ(φ) to the Hamiltonian of the surface at the sub-band bottom gives This Hamiltonian still has quantum confined eigenstates around the wire, given exactly by with = ± 1 2 , ± 3 2 , ± 5 2 , ... and the spinors ξ ± are eigenspinors of the Pauli-matrix σ x . Although the states ψ ± ( , k = 0) are no longer simultaneous eigenstates of the angular momentum operator L z = −iv F ∂ φ R W , they still have eigen-energies ε ,± = ± v F | |, the same as those of a wire with an isotropic surface chemical potentialμ. Importantly this means that the locations where the chemical potential are tuned to a sub-band edge, µ 0 = ε ,± , will be very close to those of an idealised circular wire without any inhomogeneity.
Physically, the absence of an effect at k = 0 is due to Klein tunnelling of 1d chiral modes around a wire through the additional potential δµ(φ), this potential only induces an angular dependence of the mode's phase and does not alter its energy at k = 0 [12]. Away from k = 0 the Klein-tunneling is not perfect and anistotropies can cause small shifts in the sub-band minima to finite k, most prominently for low angular momenta sub-bands, see Supplementary Fig. 6c and explanation below. Strikingly, however, the perfect Klein-tunnelling at k = 0 and strong but imperfect Kleintunnelling for small k means that there is very little change in the energies where sub-band minima occur and associated resistivity peaks at these energies, even if that anisotropy is much larger than the sub-band gap. Although Klein-tunnelling results in extremely small changes in the energies of sub-band minima due to geometric effects, the charge density -as changed in the experiment -will be dependent on these effects due to changes at large k where Klein-tunnelling is less effective and suppressed [12]. To assess the influence of these effects on density we first need to obtain an estimate of the induced charge and resulting anisotropy in chemical potential µ(φ) =μ + δµ(φ), this can be approximated using a classical electrostatic analysis. As in Ref. [12], the inhomogeneity of charge density can be estimated by numerically solving the Laplace equation for our set-up, ∇ ε(x, y)∇Φ(x, y) = 0, as shown in Supplementary Fig. 6. We use ε = 3.5 for the dielectric constant of SiO 2 and the Dirichlet conditions Φ(x, y) = V G on the gate and Φ(x, y) = 0 on the surface of the TI nanowire, which assumes the wire is perfectly metallic and hexagonal (with rounded corners, see below). The induced charge on the nanowire is proportional to the electric field normal to the wire, ∼ ∇Φ(x, y) ·r. To obtain an estimate for µ(φ), we can use the 2D Dirac relation n ∝ µ 2 . The ratio of the chemical potential to its average is shown Supplementary Fig. 6b in terms of the polar angle around the wire. We see that the flat bottom gate leads to an imbalance in chemical potential around the wire and that the hexagonal shape results in sharp features at the corners of the wire. This numerical estimate assumes that the surface of the nanowire is perfectly metallic, with no penetration of the surface state into the wire. The influence of the finite penetration of the surface state into the wire will be most apparent at the corners, which cannot anyway be treated as perfectly sharp, and at the bottom of the wire, where the wire touches the dielectric. To model this we assume a small distance of 1 nm between dielectric and bottom of the wire and round the edges of the hexagon by the equivalent of 0.2 nm, which is much smaller than estimates for the penetration depth of the surface state and so likely overestimates the sharpness of the voltage spikes at the corners. of the sub-band minima in terms of energy apart from a small splitting of the lowest angular momenta sub-band = ±1/2 away from k = 0, which results from the fact δµ (1) is largest. Since this sub-band is anyway not visible in our experiment the theoretical discussion of transport in the main text is essentially unchanged by geometric effects. In particular, the theoretical plot Fig. 3 of the main text, which shows that resistivity peaks are equally spaced as a function of chemical potential for a circular wire, remains valid for our geometry apart from very small shifts in peak positions and the onset of a finite conductivity at a slightly smaller chemical potential. Despite this, we see in Supplementary Fig. 6c that for large k the bands and their net velocity are affected by the anisotropies arising from µ(φ). Thus the charge density n(μ) is modified, which is the quantity controlled in our experiment by the gate voltage. To estimate the influence of the anisotropy on the position of the peaks as a function of charge density -rather than chemical potential -we show in panel Supplementary Fig. 6d the charge density at which the bottom of a sub-band occurs and a resistivity peak is expected (to be compared to Fig. 2d of our paper). The plot is calculated self-consistently by tuning the chemical potential to a given sub-band minimum, then recalculating all bands for the applied potential. This allows the recalculation of densities and applied voltages in an iterative procedure until convergence is reached. Supplementary Fig. 6d shows the resulting effect is small for all three of our experimental devices, which differ by the amount of gate voltage (V 0 G = 30.5, 10, and 18.5 V for device 1, 2, and 3, respectively) needed to reach charge neutrality. For negative voltages the effects are completely negligible for all three devices but for positive voltages which add to V 0 G there is a small effect visible, most pronounced for device 1 as V 0 G is largest in this case. More precisely, we find that the change of density ∆n is smaller than 0.2 2π R W which corresponds to a change of gate voltage of ≈ 4.5 V for the largest peak shift of device 1, i.e. = 7/2 for ∆V G > 0. Actually, this trend is consistent with our data (the circular-wire fits predict the position of the peak at ∆V G ≈ 24 V but it occurs experimentally at ∆V G ≈ 27 V). Although this might indicate the change of charge density is just visible in our experiment, one should also take into account experimental noise, the substantial uncertainties in our electrostatic modeling, and further small perturbations to the sub-band structure, e.g. due to non-linearities of the Dirac spectrum. For device 2 the effect of anisotropies is smaller compared to device 1, while device 3 is similar to device 1. In device 2 we find a maximal shift in density of 0.12 2π R W or ≈ 2.5 V for the largest shift in a peak, i.e. for = 7/2 of the upper Dirac cone, which occurs at ∆V G = 28 V. We therefore conclude that the effects of gate-, dielectric-and shape-induced anisotropies on charge density are small for our experimental set-up and do not affect our modelling of the experiment. In particular the presence of quantum-confined surface states in a bulk-insulating quantum wire identified by resistivity peaks with a super-linear dependence on gate voltage remains valid our set-up, but geometric effects may be relevant for future devices, in particular when a large gate voltage is necessary for accessing the Dirac point.