Quantum-Hall to Insulator Transition in Ultra-low-carrier-density Topological Insulator Films and a Hidden Phase of the Zeroth Landau Level

TL;DR

This paper reports the observation of a quantum-Hall-to-insulator transition in ultra-low-carrier-density topological insulator films, revealing a new universality class and a hidden phase at the zeroth Landau level, advancing understanding of topological surface states.

Contribution

It introduces ultra-low-carrier-density Sb₂Te₃ films enabling exploration of the zeroth Landau level and uncovers a novel quantum-Hall-to-insulator transition with a distinct universality class.

Findings

Observation of quantum-Hall-to-insulator transition near the zeroth Landau level

Discovery of a hidden phase at the zeroth Landau level

Identification of a new universality class for the transition

Abstract

A key feature of the topological surface state under a magnetic field is the presence of the zeroth Landau level at the zero energy. Nonetheless, it has been challenging to probe the zeroth Landau level due to large electron-hole puddles smearing its energy landscape. Here, by developing ultra-low-carrier density topological insulator Sb$_2$Te$_3$ films, we were able to reach an extreme quantum limit of the topological surface state and uncover a hidden phase at the zeroth Landau level. First, we discovered an unexpected quantum-Hall-to-insulator-transition near the zeroth Landau level. Then, through a detailed scaling analysis, we found that this quantum-Hall-to-insulator-transition belongs to a new universality class, implying that the insulating phase discovered here has a fundamentally different origin from those in non-topological systems.

Tables

Table 1. Table S1: Ti vapor pressure as a function of temperature provided by Veeco or on line.

Pressure (Torr)	10-8	10-7	10-6	10-5	10-4	10-3
Temperature ($°$C)	1062	1137	1227	1327	1442	1577

Table 2. Table S2: Some of the Ti fluxes at lower temperatures along with corresponding vapor pressure and calculated doping level which were used for the samples in this work. For calculating Ti flux at different temperatures, we first start by fitting the vapor pressure data as a function of temperature, with a polynomial ln ⁡ ( P ) = a 1 + a 2 ​ T + a 3 ​ T 2 𝑃 subscript 𝑎 1 subscript 𝑎 2 𝑇 subscript 𝑎 3 superscript 𝑇 2 \ln(P)=a_{1}+a_{2}T+a_{3}T^{2} . The fitting parameters a 1 = − 59.6 , a 2 = 0.051 formulae-sequence subscript 𝑎 1 59.6 subscript 𝑎 2 0.051 a_{1}=-59.6,a_{2}=0.051 , and a 3 = − 1.11 × 10 − 5 subscript 𝑎 3 1.11 superscript 10 5 a_{3}=-1.11\times 10^{-5} can thus be extracted. Also for an ideal gas, it can be shown that Φ i = Φ j ⋅ P i P j ​ ( T j T i ) 1 / 2 subscript Φ 𝑖 ⋅ subscript Φ 𝑗 subscript 𝑃 𝑖 subscript 𝑃 𝑗 superscript subscript 𝑇 𝑗 subscript 𝑇 𝑖 1 2 \Phi_{i}=\Phi_{j}\cdot\frac{P_{i}}{P_{j}}\left(\frac{T_{j}}{T_{i}}\right)^{1/2} = Φ j ⋅ a 1 + a 2 ​ T + a 3 ​ T 2 a 1 + a 2 ​ T + a 3 ​ T 2 ​ ( T j T i ) 1 / 2 ⋅ subscript Φ 𝑗 subscript 𝑎 1 subscript 𝑎 2 𝑇 subscript 𝑎 3 superscript 𝑇 2 subscript 𝑎 1 subscript 𝑎 2 𝑇 subscript 𝑎 3 superscript 𝑇 2 superscript subscript 𝑇 𝑗 subscript 𝑇 𝑖 1 2 \Phi_{j}\cdot\frac{a_{1}+a_{2}T+a_{3}T^{2}}{a_{1}+a_{2}T+a_{3}T^{2}}\left(\frac{T_{j}}{T_{i}}\right)^{1/2} where Φ i subscript Φ 𝑖 \Phi_{i} and Φ j subscript Φ 𝑗 \Phi_{j} are source fluxes (atoms/ c ​ m 2 ⋅ s ⋅ 𝑐 superscript 𝑚 2 𝑠 cm^{2}\cdot s ), P i subscript 𝑃 𝑖 P_{i} and P j subscript 𝑃 𝑗 P_{j} are vapor pressures of the source (Torr) for the source temperatures, T i subscript 𝑇 𝑖 T_{i} and T j subscript 𝑇 𝑗 T_{j} (K), respectively. Therefore, by knowing the flux Φ j subscript Φ 𝑗 \Phi_{j} at a specific temperature (usually it is measured at higher temperatures where you have reasonably measurable flux by the quartz crystal microbalance (QCM) system), Φ i subscript Φ 𝑖 \Phi_{i} can be calculated. Here, for each doped sample, we systematically vary Ti flux (temperature) while keeping Sb at a constant flux/temperature. The doping level is based on Φ Ti Φ Ti + Φ Sb × 100 % subscript Φ Ti subscript Φ Ti subscript Φ Sb percent 100 \frac{\Phi_{\text{Ti}}}{\Phi_{\text{Ti}}+\Phi_{\text{Sb}}}\times 100\% . Φ Sb subscript Φ Sb \Phi_{\text{Sb}} ≈ \approx 1.2 × 10 13 ​ c ​ m − 2 ​ s − 1 1.2 superscript 10 13 𝑐 superscript 𝑚 2 superscript 𝑠 1 1.2\times 10^{13}cm^{-2}s^{-1} at 395 ° ° \degree C and Φ Ti subscript Φ Ti \Phi_{\text{Ti}} ≈ \approx 6 × 10 12 ​ c ​ m − 2 ​ s − 1 6 superscript 10 12 𝑐 superscript 𝑚 2 superscript 𝑠 1 6\times 10^{12}cm^{-2}s^{-1} at 1550 ° ° \degree C.

Ti temperature ($°$C)	Vapor pressure ($\mu$Torr)	Flux (1010cm-2S-1)	Doping level (%)
1325	9.26	8.95	0.7
1340	12.7	12.3	1.0
1344	13.9	13.3	1.1
1350	15.7	15.1	1.3
1360	19.4	18.5	1.5

Equations

$$H(\mathbf{k})=ϵ(\mathbf{k})\mathbb{I}+v_F(k_x{\Gamma }_1+k_y{\Gamma }_2)+v_3{k}_z{\Gamma }_3+\mathcal{M}(\mathbf{k}){\Gamma }_4$$

$${\Gamma }_s={\tau }_1\otimes {\sigma }_s,{\Gamma }_4={\tau }_3\otimes {\mathbb{I}}_2,$$

$$D_j=p_j-e{A}_j,$$

$$a=\frac{{\ell }_B}{\sqrt{2}\hslash }(D_x+i{D}_y),a^{†}=\frac{{\ell }_B}{\sqrt{2}\hslash }(D_x-i{D}_y),$$

$$\displaystyle H(k_{z})=\widetilde{\epsilon}\ \mathbb{I}_{4}+\left(\begin{array}[]{cccc}\widetilde{\mathcal{M}}&v_{3}k_{z}&0&v_{F}\frac{\sqrt{2}}{\ell_{B}}a\\ v_{3}k_{z}&-\widetilde{\mathcal{M}}&v_{F}\frac{\sqrt{2}}{\ell_{B}}a&0\\ 0&v_{F}\frac{\sqrt{2}}{\ell_{B}}a^{\dagger}&\widetilde{\mathcal{M}}&-v_{3}k_{z}\\ v_{F}\frac{\sqrt{2}}{\ell_{B}}a^{\dagger}&0&-v_{3}k_{z}&-\widetilde{\mathcal{M}}\end{array}\right),$$

$$|{\Psi }_N⟩=({\alpha }_{N-1}|N-1⟩,{\beta }_{N-1}|N-1⟩,{\alpha }_N|N⟩,{\beta }_N|N⟩{)}^T,$$

$$\displaystyle H(k_{z},N)=\left(\begin{array}[]{cccc}\widetilde{\mathcal{M}}_{N-1}^{+}&v_{3}k_{z}&0&v_{F}\frac{\sqrt{2N}}{\ell_{B}}\\ v_{3}k_{z}&\widetilde{\mathcal{M}}_{N-1}^{-}&v_{F}\frac{\sqrt{2N}}{\ell_{B}}&0\\ 0&v_{F}\frac{\sqrt{2N}}{\ell_{B}}&\widetilde{\mathcal{M}}_{N}^{+}&-v_{3}k_{z}\\ v_{F}\frac{\sqrt{2N}}{\ell_{B}}&0&-v_{3}k_{z}&\widetilde{\mathcal{M}}_{N}^{-}\end{array}\right)$$

$$\left[-i{\sigma }_3{\tau }_1{v}_3{\partial }_z+{\tau }_3(M_0-M_1{\partial }_z^2)\right]|\Psi ⟩=0.$$

$$\displaystyle\left(\begin{array}[]{cc}M_{0}-M_{1}\partial_{z}^{2}&\mp iv_{3}\partial_{z}\\ \mp iv_{3}\partial_{z}&-(M_{0}-M_{1}\partial_{z}^{2})\end{array}\right)|\psi_{p}\rangle=0$$

$$H_{\text{surf}}=v_F{\eta }_3\otimes (k_x{\sigma }_2-k_y{\sigma }_1)+{\Delta }_t(k_{\perp }){\eta }_1\otimes {\mathbb{I}}_2,$$

$${\xi }_{\pm }^{-1}=\text{Re}\left[\frac{v_3\pm \sqrt{v_3^2+4M_0{M}_1}}{2M_1}\right].$$

$${\Delta }_0\approx M_1\frac{L}{{\xi }^3}{e}^{-L/\xi },{\Delta }_2\approx M_2\frac{L}{\xi }{e}^{-L/\xi }.$$

$$\displaystyle H(-i\partial_{z},N=0)=\left(\begin{array}[]{cc}\widetilde{\cal M}&iv_{3}\partial_{z}\\ iv_{3}\partial_{z}&-\widetilde{\cal M}\end{array}\right)$$

$$H_{\text{surf}}$$

$$\displaystyle=\left(\begin{array}[]{cccc}\Delta_{z}&-i\omega_{0}a^{\dagger}&{\Delta_{0}}&0\\ i\omega_{0}a&-\Delta_{z}&0&{\Delta_{0}}\\ {\Delta_{0}}&0&\Delta_{z}&i\omega_{0}a^{\dagger}\\ 0&{\Delta_{0}}&-i\omega_{0}a&-\Delta_{z}\end{array}\right)$$

$$|0,+⟩$$

$$|0,-⟩$$

$$|t_0⟩$$

$$|b_0⟩$$

$$|t_N^{+}⟩$$

$$|t_N^{-}⟩$$

$$|b_N^{+}⟩$$

$$|b_N^{-}⟩$$

$$H_N={\omega }_0\sqrt{N}\mathbb{I}\otimes {\alpha }_3+{\Delta }_0{\eta }_1\otimes {\alpha }_1+{\Delta }_z\mathbb{I}\otimes {\alpha }_1,$$

$$\pm {\epsilon }_N=\pm \sqrt{N{\omega }_0^2+({\Delta }_z\pm {\Delta }_0{)}^2}.$$

$$H=H_{\text{surf}}+V_{\text{rand}}$$

$$V_{\text{rand}}=\int d\text{r}v(\text{r}){\psi }^{†}(\text{r})\psi (\text{r}),$$

$$v(\text{r})=\sum _{j=1}^{N_{\text{imp}}}\frac{u_j}{2\pi d^2}\exp (-|\text{r}-{\text{R}}_j{|}^2/2d^2),$$

$${\Gamma }^2=8\pi u^2\frac{N_{\text{imp}}}{({\ell }_B^2+2d^2)L^2},$$

$$\overline{V(\text{r})}=0,\overline{V(\text{r})V({\text{r}}^{\prime })}=\frac{W^2}{12}\delta (\text{r}-{\text{r}}^{\prime }),$$

$$\mathcal{G}(\text{k},i{\omega }_n)=\frac{1}{i{\omega }_n-h(\text{k})-\Sigma (i{\omega }_n)}.$$

$$\tau =\underset{{\omega }_n\to 0}{\lim }1/\text{Im}[\Sigma ({\omega }_n)].$$

$$\Sigma (i{\omega }_n)=$$

Full text

See pages 1 of draft_hq See pages 2 of draft_hq See pages 3 of draft_hq See pages 4 of draft_hq See pages 5 of draft_hq See pages 6 of draft_hq See pages 7 of draft_hq See pages 8 of draft_hq See pages 9 of draft_hq See pages 10 of draft_hq See pages 11 of draft_hq See pages 12 of draft_hq See pages 13 of draft_hq See pages 14 of draft_hq See pages 15 of draft_hq See pages 16 of draft_hq See pages 17 of draft_hq See pages 18 of draft_hq See pages 19 of draft_hq See pages 20 of draft_hq See pages 21 of draft_hq See pages 22 of draft_hq

Supporting Information for

Quantum-Hall to Insulator Transition in Ultra-low-carrier-density Topological insulator Films and a Hidden Phase of the Zeroth Landau Level

Maryam Salehi, Hassan Shapourian, Ilan Thomas Rosen, Myung-Geun Han,

Jisoo Moon, Pavel Shibayev, Deepti Jain, David Goldhaber-Gordon, Seongshik Oh∗

∗Corresponding author. Email: ohsean@physics.rutgers.edu

I Supplementary Figures (Experiment)

II Supplementary Text (Theory)

In this section, we start by briefly discussing the low energy effective model of topological insulators and Dirac surface states and their Landau level spectrum in the presence of a strong magnetic field. Next, we explain the effect of random impurity potential on the energy spectrum and transport measurements.

II.1 Bulk effective model of topological insulators

We use the four-band Dirac Hamiltonian to study the low energy properties of Sb2Te3 thin films. This model was introduced in Refs.Zhang et al. (2009); Liu et al. (2010). Near the $\Gamma$ point, the effective low energy Hamiltonian is written as

[TABLE]

where the Dirac matrices are given by

[TABLE]

and $s=1,2,3$. In this convention the $\sigma_{1,2,3}$ and $\tau_{1,2,3}$ matrices act on the spin and orbital degrees of freedom respectively. Moreover, $\epsilon({\bf k})=C_{0}+C_{1}k_{z}^{2}+C_{2}k_{\perp}^{2}$ and ${\mathcal{M}}(k)=M_{0}+M_{1}k_{z}^{2}+M_{2}k_{\perp}^{2}$. The parameters $v_{F}$, $v_{3}$, $C_{i}$ and $M_{i}$ can be chosen carefully to reproduce the band structure near the $\Gamma$ point of the Sb2Te3 Liu et al. (2010). We should note that $M_{1}$ and $M_{2}$ coefficients are positive and $M_{0}$ is negative in the Bi2Se3 family of materials, e.g. Bi2Se3, Bi2Te3 and Sb2Te3.

In the presence of a magnetic field, the canonical momenta is modified into

[TABLE]

which satisfy the commutation relation, $[D_{x},D_{y}]=i{\hbar^{2}}/{\ell_{B}^{2}}$ where $\ell_{B}^{2}=\hbar/eB$. Let us introduce the ladder operators

[TABLE]

which obey $[a,a^{\dagger}]=1$. Then, the TI Hamiltonian in the presence of magnetic field becomes

[TABLE]

where $\widetilde{\mathcal{M}}(k_{z},a^{\dagger}a)=M_{0}+M_{1}k_{z}^{2}+2M_{2}/\ell_{B}^{2}(a^{\dagger}a+1/2)$ and $\widetilde{\epsilon}(k_{z},a^{\dagger}a)=C_{0}+C_{1}k_{z}^{2}+2C_{2}/\ell_{B}^{2}(a^{\dagger}a+1/2)$. Next, we bring the Hamiltonian into block diagonal form using the basis

[TABLE]

that is

[TABLE]

where $\ket{N}$ is the eigenstate of the Harmonic operator $a^{\dagger}a\ket{N}=N\ket{N}$ and $\widetilde{\mathcal{M}}_{N}^{\pm}=\widetilde{\epsilon}(k_{z},N)\pm\widetilde{\mathcal{M}}(k_{z},N)$. The zeroth Landau levels must be considered separately. They are given in the basis $\ket{\Psi_{0}}=(0,0,\alpha_{0}\ket{0},\beta_{0}\ket{0})^{T}$.

II.2 Surface effective model

In this part, we discuss the surface Dirac Hamiltonian  König et al. (2008); Zhou et al. (2008); Linder et al. (2009); Lu et al. (2010); Shan et al. (2010); Liu et al. (2010) and derive the Landau spectrum in the presence of magnetic field.

First, we need to find the zero modes which form the Dirac node on the surfaces. We model the top (bottom) surfaces by imposing an open boundary condition along the $z$-direction. This boundary condition breaks the translational symmetry and hence, $k_{z}$ is no longer a good quantum number and must be replaced by the gradient operator $-i\partial_{z}$ in real space. So, the zero modes are the solutions to the following Schröedinger equation,

[TABLE]

Since the above operator is diagonal in the spin basis, solutions can be written in the up/down spin states. Hence, we need to solve the two component equation

[TABLE]

In principle, the above equation can be solved for a system of thickness $d$ and the finite-size tunneling gap in the spectrum can be computed König et al. (2008); Zhou et al. (2008); Linder et al. (2009). Then, we do perturbation theory to the in-plane kinetic terms in the basis of zero modes and derive the surface Hamiltonian (normal to $z$ direction). We shall only quote the result. The surface Hamiltonian is found to be

[TABLE]

where $\eta_{i}$ are Pauli matrices in top (bottom) surface, $\Delta_{t}(k_{\perp})=\Delta_{0}+\Delta_{2}k_{\perp}^{2}$ denotes a tunneling amplitude between top and bottom surfaces ($k_{\perp}^{2}=k_{x}^{2}+k_{y}^{2}$). For a thick sample, the surface states decay exponentially into the bulk, where the characteristic length is

[TABLE]

The lowest order term, which gives the inter-plane tunneling, comes from the kinetic term as in

[TABLE]

Using the parameters of Liu et. al. Liu et al. (2010), M_{0}\approx-0.22\eV, M_{1}\approx 20\eVÅ2, M_{2}\approx 50\eVÅ2, v_{3}\approx 0.84\eVÅ, and v_{F}\approx 3.40\eVÅ, we obtain $\Delta_{0}\approx 6m$eV and $\Delta_{2}\approx 15$eVÅ2 for a sample of thickness $L=8nm$. We note that this calculation overestimates the finite-size energy gap $\Delta_{0}$, which is experimentally measured to be $\Delta_{0}^{\text{exp}}\approx 0.15m$eV.

In the presence of strong magnetic field, the zero energy modes on the surface are found by using the bulk zeroth Landau level Hamiltonian (14),

[TABLE]

in which $\widetilde{\cal M}=M_{0}+M_{2}/\ell_{B}^{2}-M_{1}\partial_{z}^{2}$. Note that compared with the zero-field Hamiltonian (18), the additional term $M_{2}/\ell_{B}^{2}$ makes the effective bulk gap smaller, and hence increases the penetration depth $\xi$ in (20). Given the parameter values mentioned earlier, we get $M_{2}/\ell_{B}^{2}\lesssim 16$meV for $B=20T$. This in turn leads to a negligible increase in the tunneling gap $\Delta_{0}$ of a sample with $L=8nm$ thickness.

Furthermore, the quadratic term $\Delta_{2}(k_{x}^{2}+k_{y}^{2})$ in the presence of magnetic field becomes $2\Delta_{2}/\ell_{B}^{2}(N+1/2)$. Hence, for the zeroth Landau level $N=0$, it increases the tunneling gap $\Delta_{0}$ by $\Delta_{2}/\ell_{B}^{2}\lesssim 5m$eV. This effect is much smaller than the estimation of Ref. Zhang et al. (2015), since $\Delta_{2}$ in our case is much smaller. Finally, the upper bound to the change in $\Delta_{0}$ is $5m$eV.

II.2.1 Surface Landau levels

Here, we study the Landau spectrum of the surface states. We consider the generic Hamiltonian (with already modified parameters),

[TABLE]

where $\omega_{0}=v_{F}\sqrt{2}/\ell_{B}$, is the energy scale of Landau levels and we also added a Zeeman term $\Delta_{z}$. The near zero Landau levels are given by

[TABLE]

where

[TABLE]

are zeroth Landau levels on top and bottom surfaces, respectively. The eigenstates are simply (anti-)bonding combinations of the surface zeroth Landau levels. It is important to note that the Zeeman field shifts both lowest Landau levels (LLLs) in the same direction and does not change the tunneling gap between LLLs.

For higher landau levels, we use the basis

[TABLE]

where the Hamiltonian in this subspace reads

[TABLE]

and $\alpha_{i}$’s are a set of Pauli matrices. The energy spectrum is given by

[TABLE]

It is evident from the above expression that the energy shifts associated with the Zeeman field and tunneling are smaller for higher LLs with larger $N$.

II.3 Effect of disorder on surface states

In this part, we study the effect of random impurity potential on surface states both in the presence and absence of a magnetic field. The competition between the magnetic field and random disorder potential can be characterized in terms of a dimensionless parameter $\Gamma/\hbar\omega_{0}$ which roughly speaking, compares the LL broadening and the LL spacing. Our results are summarized as $2$d phase diagrams in Figs. S12b and S14b where the corresponding zero-field systems are described by two massless Dirac cones and massive Dirac Hamiltonians due to the inter-surface tunneling gap, respectively. The important observation in either case is that the $\nu=1$ plateau is quite robust even in the strong disorder limit.

We model impurities (crystal defects, charged defects, etc.) by adding a random potential to the clean Hamiltonian of surface states (19),

[TABLE]

where

[TABLE]

and $v(\textbf{r})$ is the impurity potential profile which is given as a set of uncorrelated random numbers.

II.3.1 Landau-level broadening

Discrete Landau levels (38) are broadened due to scattering caused by the disorder potential. Following Ando (1984); Nomura et al. (2008), we use the long range disorder potential profile

[TABLE]

which consists of $N_{\text{imp}}$ impurities at random locations ${\bf R}_{j}$. To ensure the neutrality, we assume equal number of positive and negative potential energies $u_{j}=\pm u$. A measure of disorder strength is defined in terms of the LL broadening parameter,

[TABLE]

derived from self-consistent Born approximation Ando (1984). The quantity $\Gamma/\hbar\omega_{0}$ is then gives a relative ratio between LL spacing and the bandwidth. As the inset of Fig. S12a shows, when $\Gamma/\hbar\omega_{0}$ is small (e.g., the blue curve $\Gamma/\hbar\omega_{0}=0.4$), the DOS preserves its discrete form and we see the LLs are well separated. This corresponds to the presence of Hall plateaus at quantized values for the Hall conductance (e.g., blue curve in Fig. S12a). As we go to larger $\Gamma$, the LL mixing increases and in the extreme limit, the DOS becomes quite smooth (as in the yellow curve for $\Gamma/\hbar\omega_{0}=1.1$) and the Hall plateaus except for the $\nu=1$ are destroyed.

II.3.2 Finite-size tunneling

The top and bottom surfaces of a thin-film TI are always coupled via tunneling through the bulk (the tunneling amplitude however could be exponentially small, see Eq.(21)). This results in a small gap in the Dirac surface spectra. In principle, we expect to observe an insulating behavior near the charge neutrality point within the surface energy gap. However, as we see in this part random disorder smears the tunneling gap and drives the gapped system towards a metal. We investigate the gap smearing phenomena in the presence of strong magnetic field and our observation is summarized as follows: In the weak magnetic field limit (or strong disorder) $\Gamma/\hbar\omega_{0}\gg 1$ the finite-size tunneling gap is smeared into a critical metallic region, while in the strong magnetic field (or weak disorder) $\Gamma/\hbar\omega_{0}\ll 1$, we get an insulating phase between the two zeroths LLs which are separated by the finite-size tunneling gap.

We start by studying the zero-field limit analytically. We assume that $V(\textbf{r})$ are random numbers taken from a uniform (or Gaussian) distribution $[-W,W]$ and satisfy

[TABLE]

where overbar denotes the disorder average

Our goal is to see the band smearing in the presence of disorder. Because of disorder, the ensemble averaged propagator develops a self-energy part $\Sigma(i\omega_{n})$ as in

[TABLE]

It is this term $\Sigma(i\omega_{n})$ which could renormalize the original gap in $h(\textbf{k})$ and if it has an imaginary time it gives rise to scattering time

[TABLE]

The scattering time $\tau$ naturally leads to a mean-free path $l=v_{F}\tau$. In the limit $\tau<\infty$, the conduction becomes diffusive, otherwise the conduction is pseudo-ballistic. To approximate the self-energy, we use the self-consistent Born approximation

[TABLE]

which is equivalent to the propagator in the replica limit (non-crossing diagrams).

Let us plug in the surface Hamiltonian (19), and find the self-energy. The self-energy can be decomposed as $\Sigma=\Sigma_{0}\mathbb{I}+\Sigma_{1}\eta_{1}$, where $\Sigma_{0}=\Sigma_{0r}+i\Sigma_{0i}$ is complex-valued. The self-consistent equation for $\Sigma_{1}$ is given by

[TABLE]

in which we neglected $\Delta_{2}$ term. We also introduce the 3d bulk gap $\Lambda\sim M_{0}$ as a high energy cut-off. The crucial message is that $\Sigma_{1}<0$, which implies that the renormalized hybridization gap $\Delta_{0}+\Sigma_{1}$ is smaller than the original gap $\Delta_{0}$. Physically, this is due to the fact that the bands are being smeared by disorder Mong et al. (2012) and the system is driven towards a gapless metallic phase. In the above discussion, we use short range delta correlated disorder to derive this result analytically. However, it also holds for other types of scalar disorder potential.

Next, we study the effect of disorder on LLs associated with gapped Dirac surface states. Analytical calculations in this case are rather tedious and we resort to numerical investigations. We first compute the Thouless number to determine the fate of the insulating phase near the charge neutrality point. The Thouless number is a measure of longitudinal conductivity and is defined by

[TABLE]

where $\Delta\varepsilon$ is the energy shift induced by changing the boundary condition from periodic to anti-periodic. The mean level spacing is $\delta\varepsilon=1/(L^{2}\rho(\varepsilon))$ in terms of DOS $\rho(\varepsilon)$. The average energy shift is evaluated by $\braket{\Delta\varepsilon}=\exp(\overline{\ln\Delta\varepsilon})$. Figure S13a shows the results for various disorder strengths. In each panel, $g_{T}$ is plotted as a function of filling fraction $\nu$ for different system sizes $L\times L$. The peaks of $g_{T}$ indicate the center of LLs where the extended (delocalized) states reside. In fact, these peaks represent the critical metal (i.e., transition point) in the plateau-to-plateau transitions where the longitudinal conductivity is scale invariant (i.e., $g_{T}$ does not change with system size). We should note that the peaks are located at even filling fractions $\nu=2n$ in Fig. S13. The reason is here we set $\Delta_{z}=0$ as our focus is mostly the LLL physics which are just shifted by the Zeeman term (c.f. Eq.(30)). $\Delta_{z}=0$ implies that the higher LLs are two-fold degenerate as seen in Eq.(38). Between the peaks in Fig. S13, we observe that $g_{T}$ decreases as we increase the system size. This is a hallmark of an insulating behavior. In other words, $g_{T}$ curves in each panel of Fig. S13a shows a series of insulating regions separated by critical metals. These insulating regions in Fig. S13 are the integer quantum Hall (IQH) plateaus as shown in Fig. S14a.

Let us look more closely at the evolutions of $g_{T}$ peaks in the four panels of Fig. S13. For weak disorder ($\Gamma/\hbar\omega_{0}=0.1$ and $0.3$), the peaks are sharper and there are two peaks which separate the $\nu=0$ plateau from $\nu=\pm 1$ plateaus (see also blue curve of Fig. S14a). As the disorder strength is increased ($\Gamma/\hbar\omega_{0}=0.7$ and $1.0$), the peaks become wider. The insulating IQH plateaus between the peaks at higher filling fractions start to disappear while the $\nu=\pm 1$ plateaus still persist (see violet curve of Fig. S14a). This is an evidence for the robustness of $\sigma_{xy}=1$ plateaus. Moreover, the $\nu=0$ insulating phase seem to show a scale-invariant behavior at larger disorder strength. This means that the original $\nu=0$ plateau has turned into a critical metal.

We summarize our results in terms of a phase diagram in Fig. S14b. The destruction of higher LL plateaus which we also observe in our numerics is the well-known levitation phenomena where the extended states are effectively pushed to higher energies by disorder Khmelnitskii (1984); Laughlin (1984). We believe that the nature of the transition from $\nu=0$ plateau into a critical metal is similar to the disorder smearing effect explained for the zero-field regime.