Quantum oscillations and three-dimensional quantum Hall effect in ZrTe₅

We calculate the DOS of 3D ZrTe₅ under a magnetic field. With the help of the retarded Green’s function $G^{R}(\varepsilon)=(\varepsilon-H+i\eta)^{-1}$ Y.X.Wang2021 , the DOS is defined as

	$\displaystyle D(\varepsilon)$	$\displaystyle=-\frac{g}{\pi L_{z}}\sum_{k_{z}}\text{Tr}[\text{Im}G^{R}(\varepsilon,k_{z})]$
		$\displaystyle=\frac{g}{\pi L_{z}}\sum_{k_{z}}\sum_{ns\lambda}\frac{\eta}{[\varepsilon-\varepsilon_{ns\lambda}(k_{z})]^{2}+\eta^{2}},$		(6)

where $L_{z}$ is system size in the $z$ direction, $g=\frac{1}{2\pi l_{B}^{2}}$ is the LL degeneracy in the $x$-$y$ plane and can be denoted as the uniform DOS, and $\eta$ represents the linewidth broadening phenomenologically that is induced by impurity scatterings. Consider the short-range pointlike impurities, $U=u_{0}\sum_{i}\delta(\bm{r}-\bm{R}_{i})$, and the impurity concentration $n_{i}$, then in the Born approximation, the scattering time $\tau=\frac{1}{2\eta}=\frac{1}{2\pi\gamma D_{F}}$ P.Hosur ; A.A.Burkov2014 ; J.Klier ; H.W.Wang , with $\gamma=\frac{1}{2}n_{i}u_{0}^{2}$ characterizing the strength of the impurity potential, and $D_{F}$ the DOS at the Fermi energy. Note that when the band inversions are absent, $\xi=\xi_{z}=0$, the divergence in the DOS is one over the square-root type S.Kaushik ; J.Klier ; Y.X.Wang2019 . Now the LL energies show complicated dependencies on $k_{z}$, leading to the intricated behavior of the DOS, as shown below.

On the other hand, the carrier density $n_{0}$ is related to the DOS by Ashby

$\displaystyle n_{0}(\mu,B)=\int_{0}^{\infty}d\varepsilon D(\varepsilon)f(\varepsilon)+\int_{-\infty}^{0}d\varepsilon D(\varepsilon)[1-f(\varepsilon)],$

(7)

where $f(\varepsilon)=\big{[}\text{exp}(\frac{\varepsilon-\mu}{k_{B}T})+1\big{]}^{-1}$ is the Fermi distribution function, with $\mu$ being the chemical potential, $k_{B}$ the Boltzmann constant, and $T$ the temperature. In the following, we focus on zero temperature.

We also calculate the conductivity and resistivity of the system. The conductivity tensors $\sigma_{\alpha\beta}$ can be derived from the linear-response Kubo-Streda formula L.Smrcka ; G.D.Mahan ,

	$\displaystyle\sigma_{\alpha\beta}=$	$\displaystyle\frac{1}{2\pi V}\sum_{\bm{k}}\int_{-\infty}^{\infty}d\varepsilon f(\varepsilon)\Big{[}\text{Tr}\big{(}J_{\alpha}\frac{dG^{R}}{d\varepsilon}J_{\beta}(G^{A}-G^{R})$
		$\displaystyle-J_{\alpha}(G^{A}-G^{R})J_{\beta}\frac{dG^{A}}{d\varepsilon}\big{)}\Big{]},$		(8)

where $V$ is the system volume, $J_{\alpha}=e\frac{\partial H}{\partial k_{\alpha}}$ is the current density operator, and $G^{A}(\varepsilon)=(\varepsilon-H-i\eta)^{-1}$ is the advanced Green’s function. In the transport theory S.Datta , the relaxation time denotes the time that an electron experiences from the undisturbed state to the equilibrium state through impurity scatterings, while the inverse LL broadening is related to the lifetime that an electron stays in the state. Here we assume that both quantities are equal to the scattering time $\tau$ S.Galeski2021 . Such an assumption has the advantage that the conductivities can be directly determined from the band structures.

After a straightforward calculation, we obtain the longitudinal conductivity $\sigma_{xx}$ and transverse Hall conductivity $\sigma_{xy}$, which are expressed as a summation over the LL index (see Appendix A),

$\displaystyle\sigma_{xx}=$

$\displaystyle\sigma_{0}\frac{4g\eta^{2}}{L_{z}}\sum^{\prime}\frac{M_{ns\lambda;n+1,s^{\prime}\lambda^{\prime}}^{2}}{[(\mu-\varepsilon_{ns\lambda})^{2}+\eta^{2}][(\mu-\varepsilon_{n+1,s^{\prime}\lambda^{\prime}})^{2}+\eta^{2}]},$

(9)

and

	$\displaystyle\sigma_{xy}=$	$\displaystyle\sigma_{0}\frac{4g\pi}{L_{z}}\sum^{\prime}\frac{(\varepsilon_{ns\lambda}-\varepsilon_{n+1,s^{\prime}\lambda^{\prime}})^{2}-\eta^{2}}{[(\varepsilon_{ns\lambda}-\varepsilon_{n+1,s^{\prime}\lambda^{\prime}})^{2}+\eta^{2}]^{2}}M_{ns\lambda;n+1,s^{\prime}\lambda^{\prime}}^{2}$
		$\displaystyle\times[\theta(\mu-\varepsilon_{ns\lambda})\theta(\varepsilon_{n+1,s^{\prime}\lambda^{\prime}}-\mu)$
		$\displaystyle-\theta(\mu-\varepsilon_{n+1,s^{\prime}\lambda^{\prime}})\theta(\varepsilon_{ns\lambda}-\mu)].$		(10)

where $\sigma_{0}=\frac{e^{2}}{h}$ is the unit of the quantum conductivity, $\theta(x)$ is the step function, the summation sign is

$\displaystyle\sum^{\prime}=\sum_{k_{z}}\sum_{n\geq 0,s\lambda}\sum_{s^{\prime}\lambda^{\prime}},$

(11)

and the matrix element $M_{ns\lambda;n+1,s^{\prime}\lambda^{\prime}}$ is

	$\displaystyle M_{ns\lambda;n+1,s^{\prime}\lambda^{\prime}}=v(c_{ns\lambda}^{1}c_{n+1,s^{\prime}\lambda^{\prime}}^{2}-c_{ns\lambda}^{4}c_{n+1,s^{\prime}\lambda^{\prime}}^{3})$
	$\displaystyle+\frac{\sqrt{2(n+1)}\xi}{l_{B}}(-c_{ns\lambda}^{1}c_{n+1,s^{\prime}\lambda^{\prime}}^{1}+c_{ns\lambda}^{4}c_{n+1,s^{\prime}\lambda^{\prime}}^{4})$
	$\displaystyle+\frac{\sqrt{2n}\xi}{l_{B}}(c_{ns\lambda}^{2}c_{n+1,s^{\prime}\lambda^{\prime}}^{2}-c_{ns\lambda}^{3}c_{n+1,s^{\prime}\lambda^{\prime}}^{3}).$		(12)

Note that the nonvanishing matrix element $\langle ns\lambda|J_{x/y}|n^{\prime}s^{\prime}\lambda^{\prime}\rangle$ determines the common selection rules $n\rightarrow n\pm 1$ Y.Jiang ; L.You .

In the transport experiments, since the measured quantities are the resistivities $\rho_{\alpha\beta}$, we need to convert conductivities into resistivities. Using the relation ${\bm{\sigma}}\cdot{\bm{\rho}}=I$, we obtain

$\displaystyle\rho_{xx}=\frac{\sigma_{xx}}{\sigma_{xx}^{2}+\sigma_{xy}^{2}},\quad\rho_{xy}=\frac{\sigma_{xy}}{\sigma_{xx}^{2}+\sigma_{xy}^{2}}.$

(13)

In this section, we study the magnetotransport in ZrTe₅ under the condition of fixed low carrier density, which is set as $n_{0}=5.2\times 10^{16}$ cm^-3.

The calculated DOS and chemical potential are plotted in Fig. 1(a). We see that when $B^{-1}<0.58$ T^-1, the system lies in the quantum limit, with an anomalous peak in the DOS. When $B^{-1}>0.58$ T^-1, the system lies in the QO regime. In this regime, the DOS exhibits the fourfold peak structure and the chemical potential oscillates with descending amplitudes. They both show a period $\Delta_{1/B}=0.95$ T^-1. When $B^{-1}$ is asymptotically large (not shown here), the system will enter the semiclassical regime H.W.Wang , in which the energy spacing between two adjacent LLs becomes smaller than the linewidth broadening $\eta$.

The oscillation period can be understood with the Onsager’s relation L.Onsager ; D.Shoeberg . It tells us that the period is inversely proportional to the extremal cross-sectional Fermi surface area $S_{e}$ in the plane perpendicular to the magnetic field:

$\displaystyle\Delta_{1/B}^{O}=\frac{2\pi e}{S_{e}}.$

(14)

When the magnetic field is absent, the chemical potential for the fixed $n_{0}$ is $\mu_{0}=23.41$ meV (see Appendix B), as indicated by the red dotted line in Fig. 1(a). The corresponding extremal Fermi surface is plotted in the inset of Fig. 1(a). The Fermi surface is almost an ellipse with the area $S_{e}=0.01005$ nm^-2. Then we obtain $\Delta_{1/B}^{O}=0.954$ T^-1, agreeing well with the period extracted from Fig. 1(a).

To illustrate the chemical potential oscillations, we show how $\mu$ crosses the $n=3$ LLs. The LL dispersions are plotted in Figs. 1(b) and 1(c), with the magnitudes of $B^{-1}$ corresponding to the blue lines $b$ and $c$ in Fig. 1(a), respectively. Actually with increasing $B^{-1}$, $\mu$ is determined by the competition between the local DOS gaining due to the 1D LL dispersion and the uniform DOS dropping induced by the magnetic field (i.e., the LL degeneracy $g$). In Fig. 1(b), when $B^{-1}=2.5$ T^-1, we have $\mu=24.39$ meV, which meets the $\zeta$ point of the $(3+-)$ LL [see Fig. 1(b), inset]. Since the local DOS gaining surpasses the uniform DOS dropping, to conserve the carrier density, $\mu$ will decrease and thus $\mu=24.39$ meV behaves as a peak. In Fig. 1(c), when $B^{-1}=3.078$ T^-1, the chemical potential decreases to $\mu=22.74$ meV, which lies above the $\Gamma$ point of the $(3++)$ LL [see Fig. 1(c), inset]. With further increasing $B^{-1}$, the local DOS gaining is inferior to the uniform DOS dropping and thus $\mu=22.74$ meV behaves as a valley. After that, $\mu$ goes into the next oscillation. Importantly, in the local region that includes four saddle points of the $n$-th LLs, the decreasing chemical potential with $B^{-1}$ can defer its crossing over these saddle points. Such a deferring effect favors distinguishing the saddle points with weak energy difference.

In the DOS, since the peaks are related to the saddle points of the bands, here the anomalous peak in the quantum limit is caused by the $\Gamma$ point of the zeroth LL, and the fourfold peaks in the $n\geq 1$ LLs are attributed to the fact that the $\lambda=\pm 1$ branches are splitted and each branch owns two saddle points, the $\Gamma$ point and the $\zeta$ point.

We discuss the saddle point sequence that the chemical potential sweeps with increasing $B^{-1}$. If

$\displaystyle\varepsilon_{n+-}(\Gamma)>\varepsilon_{n++}(\zeta),$

(15)

the saddle point is swept in the sequence $\zeta$-$\zeta$-$\Gamma$-$\Gamma$. By contrast, if

$\displaystyle\varepsilon_{n+-}(\Gamma)<\varepsilon_{n++}(\zeta),$

(16)

the sequence will be $\zeta$-$\Gamma$-$\zeta$-$\Gamma$. For example, in the inset of Fig. 1(b), we have $\varepsilon_{3+-}(\Gamma)>\varepsilon_{3++}(\zeta)$, thus the saddle points are swept in the former sequence, which is also valid for other LLs. In a previous work J.Wang , the carrier density of the ZrTe₅ crystal is reported to be above $10^{17}$ cm^-3 and thus the quantum limit is reached at the critical magnetic field $B_{c}>10$ T. Such a strong magnetic field results in the well-separated $\lambda=\pm 1$ branches and the latter sequence $\zeta$-$\Gamma$-$\zeta$-$\Gamma$.

If we include the spin Zeeman effect into the system, the Hamiltonian is

$\displaystyle H_{Z}=-\frac{1}{2}g\mu_{B}B\sigma_{z},$

(17)

where $g$ denotes the Landé $g$-factor and $\mu_{B}$ is the Bohr magneton. Taking the value $g\simeq 10$ that is reported in ZrTe₅ R.Y.Chen ; Z.G.Chen , the total Zeeman splittings will get strengthened, leading to the further separation of the two $\lambda=\pm 1$ branches. Thus, for the lower LLs, the saddle point sequence may become $\zeta$-$\Gamma$-$\zeta$-$\Gamma$, and for the higher LLs, the sequence remains unchanged. This means that with increasing LL index $n$, there exists a transition of the saddle point sequence, which depends heavily on the carrier density and the $g$-factor of the LLs.

Next we study the conductivities. The results are plotted as a function of the inverse magnetic field $B^{-1}$ in Figs. 2(a) and 2(b). We analyze the longitudinal conductivity $\sigma_{xx}$ in Fig. 2(a). In the limiting clean case $\eta=0$, which corresponds to the infinite scattering time, $\sigma_{xx}$ vanishes. The impurity-induced scatterings can give rise to the drift current along the electric field direction, leading to the increasing of $\sigma_{xx}$ with weak $\eta$. When $\eta=0.01$ meV, $\sigma_{xx}$ exhibits clear oscillations and fourfold peaks that are similar to the DOS. When $\eta=0.1$ meV, the fourfold peaks are smoothened and become broad. Further increasing $\eta$ to be as strong as $\eta\geq 5$ meV, we observe that the oscillations in $\sigma_{xx}$ are smeared out, implying that the LL structures are broken. Moreover, in the quantum limit, $\sigma_{xx}$ will get further increased, while in the QO regime, $\sigma_{xx}$ decreases. Such an effect of impurity scatterings on $\sigma_{xx}$ is consistent with our self-consistent Born approximation study of the diagonal disorder in a Weyl semimetal system Y.X.Wang2020 .

We also analyze the Hall conductivity $\sigma_{xy}$ in Fig. 2(b). When $\eta=0$, $\sigma_{xy}$ is proportional to the inverse magnetic field $B^{-1}$, which retrieves the classical Hall conductivity expression V.Konye ,

$\displaystyle\sigma_{xy}=\frac{n_{0}e}{B}.$

(18)

With the extracted slope $k=8.31\times 10^{-2}$ m$\Omega^{-1}$ cm^-1 T, the carrier density is obtained as $n_{0}=\frac{k}{e}=5.19\times 10^{16}$ cm^-3, which is consistent with the chosen carrier density. For the weak impurity scattering $\eta<0.1$ meV, $\eta$ can be ignored in Eq. (10), thus $\sigma_{xy}$ shows certain robustness to weak $\eta$. When $\eta$ increases, the robustness disappears and $\sigma_{xy}$ will get suppressed from the higher LLs to lower ones Y.X.Wang2019 . At strong $\eta\geq 5$ meV, $\sigma_{xy}$ even becomes negative at large $B^{-1}$. This is because when the chemical potential lies between the $n$-th and $(n+1)$-th LLs, the dominate contributions to $\sigma_{xy}$ come from the LL transition $(n,1,\lambda)\rightarrow(n+1,1,\lambda)$, which is similar to two-dimensional (2D) Dirac fermion in graphene V.P.Gusynin . Since the energy spacing between the higher LLs is smaller, the associated $\sigma_{xy}$ will get suppressed at first, which will occur successively for the lower LLs. In the case when the energy spacing is less than the strong $\eta$, $\sigma_{xy}$ becomes negative.

Then we study the resistivity. To make direct comparisons with experiments, the resistivities are plotted as a function of the magnetic field $B$ in Figs. 2(c) and 2(d). With weak $\eta=0.01$ meV, the conductivities satisfy $\sigma_{xx}\ll\sigma_{xy}$, and correspondingly, the longitudinal resistivity and Hall resistivity are given as

$\displaystyle\rho_{xx}=\sigma_{xx}\sigma_{xy}^{-2},\quad\rho_{xy}=\sigma_{xy}^{-1},$

(19)

respectively. Thus in Fig. 2(c), $\rho_{xx}$ exhibits the SdH oscillations on an increasing background, which keeps the features of $\sigma_{xx}$. The fourfold peak structure of $\rho_{xx}$ is more clearly seen in the inset of Fig. 2(c). On the other hand, $\rho_{xy}$ shows a linear dependence on $B$. The linear $\rho_{xy}$ was derived by Abrikosov for Dirac electrons occupying the lowest zeroth LL Abrikosov , where the large and nonsaturated magnetoresistivity was understood with fixed carrier density. In experiment, such a linear $\rho_{xy}$ was observed in the Weyl semimetal NbP C.Shekhar and the Dirac semimetal Cd₃As₂ L.P.He ; T.Liang2015 even at a very low magnetic field.

At strong $\eta=5$ meV, $\sigma_{xx}$ becomes comparable to $\sigma_{xy}$. When the chemical potential crosses the $n$th LLs, the fourfold peaks in $\sigma_{xx}$ merge into the single peak structure. Now in Fig. 2(d), for $\rho_{xx}$, the SdH oscillations are blurred but the kinks are seen. For $\rho_{xy}$, it exhibits a ramp in the quantum limit that gradually increases with $B$. Note that $\rho_{xy}$ shows a plateau in a broad range $0.604$ T$<B<1.186$ T. In fact, such a plateau originates from the plateaus in $\sigma_{xx}$ and $\sigma_{xy}$ that are induced by impurity scatterings, as indicated by the arrows in Figs. 2(a) and 2(b). These results strongly suggest that 3D QHE cannot be explained with fixed carrier density.

In this section, we study the magnetotransport of ZrTe₅ under the condition of fixed chemical potential, which is set as $\mu=30$ meV.

The calculated DOS and carrier density are plotted in Fig. 3(a). The quantum limit is reached at the critical magnetic field $B_{c}=1.88$ T. Now as the $\Gamma$ point of the zeroth LL is always below the chemical potential,

$\displaystyle\varepsilon_{0+}(\Gamma)=M-\frac{\xi}{l_{B}^{2}}<\mu,$

(20)

the anomalous peak that appears in Fig. 1(a) is absent. In the QO regime, we can see that the fourfold peaks of the DOS are overlapping and thus are difficult to be distinguished. On the other hand, the carrier density oscillates with a period $\Delta_{1/B}=0.55$ T, implying that the charging energies exist in the system. The period can also be understood with the Onsager’s relation. Without a magnetic field, the extremal Fermi surface in the $k_{x}$-$k_{y}$ plane is plotted in the inset of Fig. 3(a), with its area $S_{e}=0.01724$ nm^-2. Then, according to Eq. (14), we have $\Delta_{1/B}^{O}=0.556$ T^-1, which shows good consistency.

The competition between the local DOS and the uniform DOS that was used to analyze the chemical potential oscillations in the previous section can be extended to illustrate the carrier density oscillations here. It is worth pointing out that the oscillations in Figs. 1(a) and 3(a) are out of phase. For example, in the insets of Figs. 1(b) and 3(b), although both the chemical potentials $\mu$ meet the $\zeta$ point of the $(3+-)$ LL, in Fig. 1(a), $\mu$ behaves as a peak, whereas in Fig. 3(a), the corresponding carrier density behaves as a valley.

Next we turn to the conductivity and resistivity. Figures 4(a) and 4(b) show that the effects of the increasing impurity scatterings $\eta$ on the conductivities $\sigma_{xx}$ and $\sigma_{xy}$ are the same as those in Figs. 2(a) and 2(b), respectively. Figure 4(c) shows that with weak linewidth $\eta=0.01$ meV, in the longitudinal resistivity $\rho_{xx}$, the SdH oscillations are evident, but the fourfold peaks can only be distinguished for the $n=1$ LLs and are overlapping for other LLs, with details seen in the inset of Fig. 4(c). This means that the saddle points of the inverted LLs cannot be well distinguished with fixed chemical potential L.You , even when the spin Zeeman effect is included. In the Hall resistivity $\rho_{xy}$, the quasi-quantized plateaus are clearly observable.

Actually, a 3D Hall system under a magnetic field can be regarded as the 2D slice stacking A.A.Burkov2011 ; Y.X.Wang2019 . Then the Hall conductivity (in unit of $\sigma_{0}$) is a summation of the Chern numbers of the occupied LLs (see Appendix C),

$\displaystyle\sigma_{xy}=\sigma_{0}\sum_{ns\lambda}\int_{-\pi}^{\pi}\frac{dk_{z}}{2\pi}C_{ns\lambda}^{k_{z}}\theta[\varepsilon_{ns\lambda}(k_{z})<\mu],$

(21)

in which the Chern number of the $ns\lambda$ LL at wave vector $k_{z}$ is $C_{ns\lambda}^{k_{z}}=1$ for Dirac fermions. In the quantum limit, the Fermi wave vector $k_{F0}$ is obtained as

$\displaystyle k_{F0}\simeq\sqrt{\frac{\mu+M}{\xi_{z}}-\frac{\xi}{\xi_{z}l_{B}^{2}}}\simeq\sqrt{\frac{\mu+M}{\xi_{z}}}=0.43\text{ nm}^{-1}.$

(22)

in which the $B$-dependent term $\frac{\xi}{\xi_{z}l_{B}^{2}}$ is much weaker than $\frac{\mu+M}{\xi_{z}}$ and can be neglected. Then we have $\sigma_{xy}=\sigma_{0}\frac{k_{F0}}{\pi}$, leading to the quasi-quantized plateau in $\rho_{xy}$ that is estimated to be $\rho_{xy}=\sigma_{xy}^{-1}=18.9$ m$\Omega$ cm, as shown in Fig. 4(c). For the $n\geq 1$ LL, if occupied, the Fermi wavevector $k_{Fn}$ is obtained as

$\displaystyle k_{Fn}\simeq\Big{(}\frac{1}{\xi_{z}}\sqrt{\mu^{2}-\frac{2nv^{2}}{l_{B}^{2}}}+\frac{M}{\xi_{z}}\Big{)}^{\frac{1}{2}}.$

(23)

Here, as the $B$-dependent term cannot be neglected, $k_{Fn}$ decreases with $B$. Thus in Fig. 4(c), the $n\geq 1$ plateau of $\rho_{xy}$, which is given as $\rho_{xy}=\pi\sigma_{0}^{-1}(k_{F0}+2\sum_{n}k_{Fn})^{-1}$, with the factor $2$ accounting for the two LL branches, increases slowly with $B$. For the kink of $\rho_{xy}$ around $B=1.8$ T in Figs. 4(c) and 4(d), it occurs when the $\lambda=-1$ branch is occupied and the $\lambda=1$ branch is unoccupied, with the magnitude $\rho_{xy}=\pi\sigma_{0}^{-1}(k_{F0}+k_{F1})^{-1}$.

In Fig. 4(d) when $\eta$ increases to 1 meV, we observe that in $\rho_{xx}$, the fourfold peaks no longer exist, but the SdH oscillations are still discernible, while in $\rho_{xy}$, the Hall plateaus are preserved, with the neighboring plateau transitions becoming smooth. More importantly, the nonvanishing dips in $\rho_{xx}$ correspond to the plateaus in $\rho_{xy}$, and the peaks in $\rho_{xx}$ correspond to the plateau transition regions in $\rho_{xy}$. Although the peak heights in $\rho_{xx}$ and the plateau magnitudes in $\rho_{xy}$ depend heavily on the model parameters, the numerical results are qualitatively consistent with the behaviors of $\rho_{xx}$ and $\rho_{xy}$ observed experimentally in ZrTe₅ F.Tang and HfTe₅ P.Wang , which thus capture the main features of 3D QHE.

In a 2D electron system, the magnetotransport is often described by fixed chemical potential, resulting in the well-defined Hall plateaus V.P.Gusynin ; D.N.Sheng . Here also under the condition of fixed chemical potential, we can understand 3D QHE as a close analog to 2D QHE, in which the quasi-quantized plateaus in $\rho_{xy}$ are determined by the dependence of the Fermi wavevectors on the magnetic field. Moreover, the impurity scatterings play indispensable roles in driving 3D QHE, since they can enhance the peak heights in $\rho_{xx}$ as well as smoothen the plateau transitions in $\rho_{xy}$. Although our results based on a non-interacting model are quite different from the interaction-induced CDW scenario F.Tang ; P.Wang ; F.Qin , but are supported by a recent study S.Galeski2021 , in which the CDW states are demonstrated to be absent in ZrTe₅.

Note that in a Dirac semimetal $(\xi=\xi_{z}=0)$ with fixed chemical potential, the Fermi wavevectors are solved as

$\displaystyle k_{F0}=\frac{1}{v_{z}}\sqrt{\mu^{2}-M^{2}},$

(24)

for the $n=0$ LL, and

$\displaystyle k_{Fn}=\frac{1}{v_{z}}\sqrt{\mu^{2}-M^{2}-\frac{2nv^{2}}{l_{B}^{2}}},$

(25)

for the $n\geq 1$ LL. We see that the Fermi wavevectors in Eqs. (24) and (25) show similar dependence on the magnetic field as those in Eqs. (22) and (23). Based on these results, we suggest that the proposed mechanism for 3D QHE is justified for a 3D Dirac fermion system, no matter the ground state is a strong TI ($\xi>0$ and $\xi_{z}>0$), a weak TI ($\xi>0$ and $\xi_{z}<0$), or even a Dirac semimetal.

The conductivity tensors can be derived from the Kubo-Streda formula in Eq. (8). At zero temperature, the longitudinal conductivity $\sigma_{xx}$ is written as

	$\displaystyle\sigma_{xx}$	$\displaystyle=-\frac{1}{\pi V}\sum_{\bm{k}}\int d\varepsilon\frac{df(\varepsilon)}{d\varepsilon}\text{Tr}\Big{[}J_{x}\text{Im}G^{R}(\varepsilon)J_{x}\text{Im}G^{R}(\varepsilon)\Big{]}$
		$\displaystyle=\frac{1}{\pi V}\sum_{\bm{k}}\text{Tr}\Big{[}J_{x}\text{Im}G^{R}(\mu)J_{x}\text{Im}G^{R}(\mu)\Big{]},$		(A1)

in which the identity $\frac{df(\varepsilon)}{d\varepsilon}=-\delta(\varepsilon-\mu)$ is used in the second line. Note that the expression of $\sigma_{xx}$ in Eq. (A1) is the same as those used in previous studies H.W.Wang ; Y.X.Wang2020 . Consider the LL degeneracy $g=\frac{1}{2\pi l_{B}^{2}}$, and insert the energies and wavefunctions of the LLs, then we have

	$\displaystyle\sigma_{xx}$	$\displaystyle=\frac{g\eta^{2}}{\pi L_{z}}\sum_{k_{z}}\sum_{ns\lambda}\sum_{n^{\prime}s^{\prime}\lambda^{\prime}}|\langle ns\lambda|J_{x}|n^{\prime}s^{\prime}\lambda^{\prime}\rangle|^{2}$
		$\displaystyle\times[(\mu-\varepsilon_{ns\lambda})^{2}+\eta^{2}]^{-1}[(\mu-\varepsilon_{n^{\prime}s^{\prime}\lambda^{\prime}})^{2}+\eta^{2}]^{-1}.$		(A2)

Similarly, the Hall conductivity $\sigma_{xy}$ is

	$\displaystyle\sigma_{xy}=$	$\displaystyle\frac{g}{2\pi L_{z}}\sum_{k_{z}}\sum_{ns\lambda}\sum_{n^{\prime}s^{\prime}\lambda^{\prime}}\int d\varepsilon f(\varepsilon)\Big{[}\langle ns\lambda|J_{x}\frac{dG^{R}}{d\varepsilon}|n^{\prime}s^{\prime}\lambda^{\prime}\rangle\langle n^{\prime}s^{\prime}\lambda^{\prime}|J_{y}(G^{A}-G^{R})|ns\lambda\rangle-\langle ns\lambda|J_{x}(G^{A}-G^{R})|n^{\prime}s^{\prime}\lambda^{\prime}\rangle$
		$\displaystyle\times\langle n^{\prime}s^{\prime}\lambda^{\prime}|J_{y}\frac{dG^{A}}{d\varepsilon}|ns\lambda\rangle\Big{]}$
	$\displaystyle=$	$\displaystyle\frac{g}{L_{z}}\sum_{k_{z}}\sum_{ns\lambda}\sum_{n^{\prime}s^{\prime}\lambda^{\prime}}\text{Im}\Big{[}\frac{\langle ns\lambda|J_{x}|n^{\prime}s^{\prime}\lambda^{\prime}\rangle\langle n^{\prime}s^{\prime}\lambda^{\prime}|J_{y}|ns\lambda\rangle}{(\varepsilon_{ns\lambda}-\varepsilon_{n^{\prime}s^{\prime}\lambda^{\prime}}+i\eta)^{2}}-\frac{\langle ns\lambda|J_{y}|n^{\prime}s^{\prime}\lambda^{\prime}\rangle\langle n^{\prime}s^{\prime}\lambda^{\prime}|J_{x}|ns\lambda\rangle}{(\varepsilon_{ns\lambda}-\varepsilon_{n^{\prime}s^{\prime}\lambda^{\prime}}-i\eta)^{2}}\Big{]}\theta(\mu-\varepsilon_{ns\lambda})\theta(\varepsilon_{n^{\prime}s^{\prime}\lambda^{\prime}}-\mu).$		(A3)

In previous literatures, the expressions of $\sigma_{xy}$ are much more complicated H.W.Wang ; B.Fu ; C.Wang and include a normal part and an anomalous part, where the former is determined by the states around the Fermi level and the latter reflects the thermodynamic property of the system. Here Eq. (A3) is quite concise and thus are facilitate to the numerical calculations.

After calculating the current density matrix elements $\langle ns\lambda|J_{x/y}|n^{\prime}s^{\prime}\lambda^{\prime}\rangle$, we obtain

$\displaystyle\sigma_{xx}=$

$\displaystyle\sigma_{0}\frac{2g\eta^{2}}{L_{z}}\sum^{\prime}\frac{M_{ns\lambda;n+1,s^{\prime}\lambda^{\prime}}^{2}}{[(\mu-\varepsilon_{ns\lambda})^{2}+\eta^{2}][(\mu-\varepsilon_{n+1,s^{\prime}\lambda^{\prime}})^{2}+\eta^{2}]}+\sigma_{0}\frac{2g\eta^{2}}{L_{z}}\sum^{\prime\prime}\frac{M_{ns\lambda;n-1,s^{\prime}\lambda^{\prime}}^{2}}{[(\mu-\varepsilon_{ns\lambda})^{2}+\eta^{2}][(\mu-\varepsilon_{n-1,s^{\prime}\lambda^{\prime}})^{2}+\eta^{2}]},$

(A4)

and

	$\displaystyle\sigma_{xy}=$	$\displaystyle\sigma_{0}\frac{4g\pi}{L_{z}}\sum^{\prime}\frac{(\varepsilon_{ns\lambda}-\varepsilon_{n+1,s^{\prime}\lambda^{\prime}})^{2}-\eta^{2}}{[(\varepsilon_{ns\lambda}-\varepsilon_{n+1,s^{\prime}\lambda^{\prime}})^{2}+\eta^{2}]^{2}}M_{ns\lambda;n+1,s^{\prime}\lambda^{\prime}}^{2}\theta(\mu-\varepsilon_{ns\lambda})\theta(\varepsilon_{n+1,s^{\prime}\lambda^{\prime}}-\mu)$
		$\displaystyle-\sigma_{0}\frac{4g\pi}{L_{z}}\sum^{\prime\prime}\frac{(\varepsilon_{ns\lambda}-\varepsilon_{n-1,s^{\prime}\lambda^{\prime}})^{2}-\eta^{2}}{[(\varepsilon_{ns\lambda}-\varepsilon_{n-1,s^{\prime}\lambda^{\prime}})^{2}+\eta^{2}]^{2}}M_{ns\lambda;n-1,s^{\prime}\lambda^{\prime}}^{2}\theta(\mu-\varepsilon_{ns\lambda})\theta(\varepsilon_{n-1,s^{\prime}\lambda^{\prime}}-\mu),$		(A5)

where the summation signs are

$\displaystyle\sum^{\prime}=\sum_{k_{z}}\sum_{n\geq 0,s\lambda}\sum_{s^{\prime}\lambda^{\prime}},\quad\sum^{\prime\prime}=\sum_{k_{z}}\sum_{n\geq 1,s\lambda}\sum_{s^{\prime}\lambda^{\prime}},$

(A6)

and the matrix element $M_{ns\lambda;n+1,s^{\prime}\lambda^{\prime}}$ is given in Eq. (12).

We note that the components of the conductivity satisfy the following relations:

	$\displaystyle\sigma_{xx}(ns\lambda\rightarrow n+1,s^{\prime}\lambda^{\prime})=\sigma_{xx}(n+1,s^{\prime}\lambda^{\prime}\rightarrow ns\lambda),$		(A7)
	$\displaystyle\sigma_{xy}(ns\lambda\rightarrow n+1,s^{\prime}\lambda^{\prime})=-\sigma_{xy}(n+1,s^{\prime}\lambda^{\prime}\rightarrow ns\lambda),$		(A8)

implying that the contributions to $\sigma_{xx}$ from the LL transition $ns\lambda\rightarrow(n+1,s^{\prime}\lambda^{\prime})$ and from $(n+1,s^{\prime}\lambda^{\prime})\rightarrow ns\lambda$ are equal, while those to $\sigma_{xy}$ are opposite. According to these properties, we can change the index $n\rightarrow n+1$, $s\leftrightarrow s^{\prime}$, and $\lambda\leftrightarrow\lambda^{\prime}$ in the second summation of Eqs. (A4) and (A5), and obtain the final expressions in Eqs. (9) and (10) in the main text.

When the magnetic field is absent, the energies for the Hamiltonian $H(\bm{k})$ in Eq. (1) can be obtained directly, which are given as

	$\displaystyle\varepsilon_{s\bm{k}}=s\varepsilon_{\bm{k}}=$	$\displaystyle s\Big{[}v^{2}(k_{x}^{2}+k_{y}^{2})+v_{z}^{2}k_{z}^{2}+\big{(}M-\xi(k_{x}^{2}+k_{y}^{2})$
		$\displaystyle-\xi_{z}k_{z}^{2}\big{)}^{2}\Big{]}^{\frac{1}{2}},$		(A9)

where $s=\pm 1$ denotes the conduction/valence band. Note that both the conduction and valence bands own twofold degeneracy.

The DOS is calculated as,

$\displaystyle D(\varepsilon>0)=\frac{g_{b}\eta}{8\pi^{4}}\iiint dk_{x}dk_{y}dk_{z}\frac{1}{(\varepsilon-\varepsilon_{\bm{k}})^{2}+\eta^{2}},$

(A10)

where $g_{b}=2$ denotes the twofold degeneracy. To do the integrations, we make the following substitutions Y.X.Wang2020 : $x=vk_{x}$, $y=vk_{y}$, $z=v_{z}k_{z}$, and $\xi^{\prime}=\frac{\xi}{v^{2}}$, $\xi_{z}^{\prime}=\frac{\xi_{z}}{v_{z}^{2}}$, and change the Cartesian into the Cylindrical system: $x=\rho\text{cos}\theta$, $y=\rho\text{sin}\theta$. After completing the integration over the azimuth angle $\theta$, we have

	$\displaystyle D(\varepsilon)=\frac{g_{b}\eta}{4\pi^{3}v^{2}v_{z}}\int_{-\infty}^{\infty}dz\int_{0}^{\infty}d\rho\rho$
	$\displaystyle\times\Big{[}\Big{(}\varepsilon-\sqrt{\rho^{2}+z^{2}+(M-\xi^{\prime}\rho^{2}-\xi_{z}^{\prime}z^{2})^{2}}\Big{)}^{2}+\eta^{2}\Big{]}^{-1}.$		(A11)

On the other hand, the carrier density is defined in Eq. (7), which, at zero temperature, is written as

$\displaystyle n_{0}$

$\displaystyle=\int_{0}^{\mu}d\varepsilon D(\varepsilon).$

(A12)

Evidently, the integrands in Eqs. (A11) and (A12) show complicated dependencies on $\rho$ and $z$, thus the integration cannot be solved analytically. But we can resort to the numerical calculations. In Figs. 5(a) and 5(b), we plot the calculated DOS as a function of the energy $\varepsilon$ and the carrier density $n_{0}$ as a function of the chemical potential $\mu$, respectively. We see that the DOS increase with $\varepsilon$ and $n_{0}$ increases with $\mu$, both in a monotonous tendency. In Fig. 5(b), the red dashed lines denote that the chosen carrier density $n_{0}=5.2\times 10^{16}$ cm^-3 in Sec. IVA corresponds to the chemical potential $\mu=23.41$ meV.

Here we give a detailed derivation of Eq. (21). The 2D Chern number of the $n$-th band is defined as an integration over the BZ,

$\displaystyle C_{n}$

$\displaystyle=-\frac{1}{2\pi i}\int_{\text{BZ}}d^{2}k\Omega_{n},$

(A13)

where $\Omega_{n}(\bm{k})=\nabla_{\bm{k}}\times A_{n}(\bm{k})\big{|}_{z}$ is the Berry curvature and $A_{n}(\bm{k})=\langle u_{n}|\nabla_{\bm{k}}|u_{n}\rangle$ is the Berry connection. Then we have

	$\displaystyle C_{n}$	$\displaystyle=-\frac{1}{\pi}\text{Im}\sum_{n^{\prime}}\int_{\text{BZ}}d^{2}k\langle\frac{\partial u_{n}}{\partial k_{x}}|u_{n^{\prime}}\rangle\langle u_{n^{\prime}}|\frac{\partial u_{n}}{\partial k_{y}}\rangle$
		$\displaystyle=-\frac{1}{\pi}\text{Im}\sum_{n^{\prime}}\int_{\text{BZ}}d^{2}k\frac{\langle u_{n}|v_{x}|u_{n^{\prime}}\rangle\langle u_{n^{\prime}}|v_{y}|u_{n}\rangle}{(\varepsilon_{n}-\varepsilon_{n^{\prime}})^{2}},$		(A14)

where $\varepsilon_{n}$ and $u_{n}$ are the energy and wave function of the $n$-th band. Note that in the second line of Eq. (A14), the identity $\langle\frac{\partial u_{n}}{\partial k_{\alpha}}|u_{n^{\prime}}\rangle=\frac{\langle u_{n}|v_{\alpha}|u_{n^{\prime}}\rangle}{\varepsilon_{n}-\varepsilon_{n^{\prime}}}$ is used. Comparing Eq. (A14) with Eq. (A3), we obtain Eq. (21) in the main text.