Magnetic field independent shape of the zero-energy landau levels in a disordered T₃ model

The shape of the Landau levels (LLs) in a disordered two-dimensional electronic system is important to many physical measurements, such as specific heat [1, 2], magnetization [3], magnetocapacitance [4], and magnetoconductivity [5–9]. This subject has been theoretically and experimentally studied for many decades [10–22]. Based on the self-consistent Born approximation (SCBA), the shape of the lowest LL in the conventional two-dimensional electronic systems has been argued to be a semi-elliptical form, which is insensitive to the distribution of disorder [14]. Later, by mapping this problem onto that of a zero-dimensional ϕ⁴ model, a general nonperturbative approach was proposed to obtain an exact quasi-Gaussian shape for the lowest LL in the presence of the white-noise disorder [11]. After that, a supersymmetry method was extended to discuss the issue in the presence of several different kinds of short-range disorder distributions [12]. It is found that the density of states (DOS) of the lowest LL obeys the following scaling behavior [13]

$$\begin{eqnarray}&&\rho (E)=\sqrt{\displaystyle \frac{B}{\gamma }}f\left(\displaystyle \frac{E}{\sqrt{B\gamma }}\right),\end{eqnarray} \tag{ 1 }$$

where B and γ denote the magnetic field strength and disorder distribution variance, respectively. The function form f is closely related to the disorder distribution. Equation (1) also reflects that the width of the LLs Γ depends on both magnetic field strength and disorder strength as ${\rm{\Gamma }}\sim \sqrt{B\gamma }$. In our previous studies, we numerically computed the shape of the zero-energy LL in disordered graphene, whose low-energy quasiparticle is governed by the Dirac–Weyl equation with pseudospin-1/2 [18, 19]. Interestingly, our results are similar to the predictions of Wegner in the conventional two-dimensional electron systems.

Recently, there is a growing interest in studying the T₃ model [23–33]. Experimentally, the lattice can be realized in a trilayer structure of the face-centered cubic lattice, such as SrTiO₃/SrIrO₃/SrTiO₃ heterostructures [34, 35]. The low-energy quasiparticles of those systems can be well described by the pseudospin-1 Dirac–Weyl equation, predicting the existence of the flat band in their electronic structures. Such a special spectrum leads to distinctive features, such as super-Klein tunneling [29, 32], diverging dc conductivity [30], and anomalous Anderson localization [36]. A natural question arises as to how such a special spectrum affects the shape of the zero-energy LL with disorders. We notice that the disorder effect in the T₃ lattice has been discussed by Vidal et al in the large-flux limit (∼h/e) [37]. However, the study in the opposite small-flux limit (∼0.001h/e) is still lacking. In this case, the magnetic field is still strong because of the small lattice spacing.

In this work, using the Lanczos recursion method, we numerically determine the exact DOS of the T₃ lattice under a strong magnetic field with the short-range disorder. To illuminate the effect of the disorder on the shape of the LLs, we perform extensive simulations for two different types of disorders: the Gaussian random potential and the Anderson-type uniform random potential [38, 39]. Remarkably, our simulations show that the DOS of the zero-energy LL obeys a scaling law, which is similar to equation (1) but independent of the magnetic field. Concretely, the shape of the zero-energy LL under the Gaussian disorder can be well fitted by the Gaussian function, while that under the uniform distribution is approximately described by a semi-elliptical form. We also present an analytical study to compute the DOS of zero-energy LL based on the SCBA. As a result, we find that the magnetic field independent scaling behavior can be attributed to the magnetic field independent macroscopic degeneracy of the zero-energy LL. Moreover, to demonstrate the significance of the shape of the LLs, the magnetoconductivity is calculated using the linear-response formula [40]. It is shown that the shape of the LLs plays an important role in the transport properties, and the above scaling behavior leads to an interesting fact that the magnetoconductivity as a function of carrier density is independent of the disorder strength. In addition, the existence of the huge magnetic field independent macroscopic degeneracy of zero-energy LL results in blue shifts of the conduction and valence bands LLs, in contrast with what found in the disordered graphene [18, 19].

This paper is structured as follows. In section 2 the tight-binding Hamiltonian and low-energy effective Hamiltonian of the T₃ model is specified, the simulation method and parameters are given, and then the calculated results of magnetoconductivity using the linear-response formula are presented. In section 3, we discuss the DOS of the zero-energy LL from our simulations, where the scaling law is proved within numerical accuracy. Finally, we summarize our main outcomes in section 4.

2.1. Tight-binding model and simulation method

The geometry of the T₃ lattice resembles the honeycomb lattice and accordingly graphene, but with an additional atom that sits at the center of each hexagon, as shown in figure 1(a). Two of these sites are threefold coordinated, denoted by A and B respectively, while the third site is connected to six nearest neighbors (denoted by H) [24]. The nearest-neighbor hopping takes place between sites A, B, and H with strength $t/\sqrt{2}$, while the hopping between sites A and B is forbidden. Thus the T₃ lattice is bipartite and particle-hole symmetric [28]. The tight-binding Hamiltonian can be written as

$$\begin{eqnarray}&&H=\displaystyle \frac{t}{\sqrt{2}}\displaystyle \sum _{\langle i,j\rangle }{c}_{{Hi}}^{\dagger }({c}_{{Aj}}+{c}_{{Bj}})+{\rm{h}}.{\rm{c}}.,\end{eqnarray} \tag{ 2 }$$

where c_Ai, c_Bi and c_Hi denotes the annihilation operator for a particle of unit cell i on site A, B and H, respectively.

Zoom In Zoom Out Reset image size

The short-range disorder potential is modeled by random on-site energy

$$\begin{eqnarray}&&V=\displaystyle \sum _{i,X=A,B,C}{\epsilon }_{i}{c}_{{Xi}}^{\dagger }{c}_{{Xi}},\end{eqnarray} \tag{ 3 }$$

where 's are random variables and the disorder is assumed to be uncorrelated, that is, $\langle {\epsilon }_{i}{\epsilon }_{j}\rangle =\langle {\epsilon }_{i}\rangle \langle {\epsilon }_{j}\rangle$. In this work, two probability distributions of on-site energy are considered: the Gaussian distribution with standard derivation σ and the uniform distribution with an interval [−W/2, W/2]. Their variances γ equal to σ² and W²/12, respectively. The applied perpendicular magnetic field B = (0, 0, B) and the Landau gauge is chosen as A = (0, Bx, 0) in our simulations. In the tight-binding model, the magnetic field can be introduced via Peierls substitution: the hopping term t is multiplied by a phase factor ${{\rm{e}}}^{{\rm{i}}{\phi }_{{mn}}}$ involving the vector potential A

$$\begin{eqnarray}&&{\phi }_{{mn}}=\displaystyle \frac{2\pi {\rm{e}}}{h}{\int }_{m}^{n}{\bf{A}}{\rm{d}}{\bf{l}},\end{eqnarray} \tag{ 4 }$$

where the integration path is the shortest distance over the nearest neighbor bond from site m to site n.

The DOS is calculated with the Lanczos recursion method [18, 19]. The averaged DOS can be evaluated from the Green's function

$$\begin{eqnarray}&&\rho (E)=-\displaystyle \frac{1}{\pi }\overline{\langle {\psi }_{0}| \mathrm{Im}\displaystyle \frac{1}{E-H+i\eta }}| {\psi }_{0}\rangle .\end{eqnarray} \tag{ 5 }$$

The small artificial cutoff energy η is chosen as η = 10⁻³t, which is rather smaller than the broadening of the LLs investigated. The random phase state is used to reduce the time complexity [41]. The random phase state is defined through

$$\begin{eqnarray}&&| {\psi }_{0}\rangle =\displaystyle \frac{1}{\sqrt{N}}\displaystyle \sum _{n}{{\rm{e}}}^{{\rm{i}}{\theta }_{n}}| n\rangle ,\end{eqnarray} \tag{ 6 }$$

where N denotes the total number of sites and {θ_n} are a set of independent random variables uniformly distributed on [0, 2π]. Because of the large systems and weak disorder considered in this work, one or few different random phase states is sufficient to obtain a good estimation of the averaged DOS.

In our simulations, the shape of the finite cluster is diamond as shown in figure 1 (a) and the periodic boundary condition is imposed. The cluster comprises Ω/A = 4000² unit cells in this work, where Ω is the area of the system, $A=\sqrt{3}{a}^{2}/2$ denotes unit cell area, and a is the lattice constant as shown in figure 1(a). For simplicity, we define a unit magnetic field strength as ${B}_{0}=6h/e\sqrt{A{\rm{\Omega }}}$, i.e. the magnetic flux threading through per unit cell $\sim \tfrac{3}{2000}\tfrac{h}{e}$ in our simulations. If using the lattice constant of graphene a = 2.47 Å, we obtain B₀ ∼ 50T which is accessible experimentally. All calculations are carried out on a 2.00 GHz Xeon 3 CPU work station, and each simulation for a single realization of a random phase state under specific disorder realization requires around 20 h of computing time running on 32 cores.

2.2. Low-energy effective Hamiltonian and LLs

In the absence of the magnetic field, the energy spectrum of the T₃ lattice is plotted in figure 1(b), which contains three branches, including two electron–hole symmetric conduction and valence bands, and a unique flat band at the charge neutrality point. The dispersion of conduction and valence bands are identical to the bands in graphene and features two inequivalent Dirac cones.

In the presence of the magnetic field, here, we employ its effective Hamiltonian obtained by a long-wavelength approximation from the tight-binding Hamiltonian to describe the low-energy properties. For a given Dirac point, the Hamiltonian reads

$$\begin{eqnarray}&&{H}_{K}={\hslash }{v}_{F}{\bf{S}}\cdot ({\bf{k}}-e{\bf{A}}),\end{eqnarray} \tag{ 7 }$$

where ${v}_{F}=\tfrac{\sqrt{3}}{2}\tfrac{{at}}{{\hslash }}$ is the Fermi velocity, k = (k_x, k_y) is the momentum operator, and S denotes the pseudospin

$$\begin{eqnarray}{S}_{x}=\displaystyle \frac{1}{\sqrt{2}}\left(\begin{array}{ccc}0 & 1 & 0\\ 1 & 0 & 1\\ 0 & 1 & 0\end{array}\right),{S}_{y}=\displaystyle \frac{1}{\sqrt{2}}\left(\begin{array}{ccc}0 & -{\rm{i}} & 0\\ {\rm{i}} & 0 & -{\rm{i}}\\ 0 & {\rm{i}} & 0\end{array}\right).\end{eqnarray} \tag{ 8 }$$

Combining with

$$\begin{eqnarray}{S}_{z}=\left(\begin{array}{ccc}1 & 0 & 0\\ 0 & 0 & 0\\ 0 & 0 & -1\end{array}\right),\end{eqnarray} \tag{ 9 }$$

S_x, S_y form a 3-dimensional representation of SU(2), satisfying commutation relation [S_i, S_j] = i_ijkS_k, similar to the Pauli matrix. Although they have higher pseudospin, they can not form a Clifford algebra {S_i, S_j} $\ne$ 2δ_ij [29].

The corresponding LLs energies are expressed as

$$\begin{eqnarray}&&{\epsilon }_{{sn}}=s\displaystyle \frac{{\hslash }{v}_{F}}{l}\sqrt{2n-1},\end{eqnarray} \tag{ 10 }$$

where $l=\sqrt{{\hslash }/{eB}}$ is the magnetic length, s = ± denotes the conduction band and valence band respectively with the LL index n = 1, 2, ..., and s = 0 labels the LL at zero energy with index n = 0, 2, 3, ...[29, 42, 43]. The total number of the Landau indexes N_Φ for each band is determined by the semiclassical quantization rule [44]

$$\begin{eqnarray}&&2({N}_{{\rm{\Phi }}}+\displaystyle \frac{1}{2})\displaystyle \frac{2\pi {eB}}{{\hslash }}=\displaystyle \frac{{\left(2\pi \right)}^{2}}{A},\end{eqnarray} \tag{ 11 }$$

where the factor 2 stems from the valley degrees of freedom. Since N_Φ is large, we obtain N_Φ = πl²/A.

The corresponding eigenfunctions of the low-energy effective Hamiltonian in equation (7) are given by

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{{{snk}}_{y}}(x,y)={C}_{s}\left(\begin{array}{c}{a}_{{sn}}{\phi }_{n-2}(x-{k}_{y}{l}^{2})\\ s{\phi }_{n-1}(x-{k}_{y}{l}^{2})\\ {b}_{{sn}}{\phi }_{n}(x-{k}_{y}{l}^{2})\end{array}\right)\displaystyle \frac{{{\rm{e}}}^{-{{\rm{i}}{k}}_{y}y}}{\sqrt{{L}_{y}}},\end{eqnarray} \tag{ 12 }$$

where k_y is the wave number in the y direction, ϕ_n (n = 0, 1, 2, ...) are eigenfunctions of the one-dimensional harmonic oscillator, and for convenience we set ϕ_n = 0 if n is negative. Here, parameters C_s, a_sn, b_sn are given by

$$\begin{eqnarray}{a}_{{sn}}=\left\{\begin{array}{cc}\sqrt{\displaystyle \frac{n}{2n-1}} & \mathrm{for}\ {\text{}}s\,=\,0,{\text{}}n\,=\,0,2,3,...\\ \sqrt{\displaystyle \frac{n-1}{2n-1}} & \mathrm{for}\ {\text{}}s\,=\,\pm ,{\text{}}n\,=\,1,2,3,...,\end{array}\right.\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}{b}_{{sn}}=\left\{\begin{array}{cc}-\sqrt{\displaystyle \frac{n-1}{2n-1}} & \mathrm{for}\ {\text{}}s\,=\,0,{\text{}}n\,=\,0,2,3,...\\ \sqrt{\displaystyle \frac{n}{2n-1}} & \mathrm{for}\ {\text{}}s\,=\,\pm ,{\text{}}n\,=\,1,2,3,...,\end{array}\right.\end{eqnarray} \tag{ 14 }$$

and

$$\begin{eqnarray}{C}_{s}=\left\{\begin{array}{cc}1 & \mathrm{for}\ {\text{}}s\,=\,0\\ \displaystyle \frac{1}{\sqrt{2}} & \mathrm{for}\ {\text{}}s\,=\,\pm .\end{array}\right.\end{eqnarray} \tag{ 15 }$$

Notice that equations (13)–(14) are well defined only when s = ± with n = 1, 2, ...and s = 0 with n = 0, 2, 3, ..., otherwise, their values are 0. For example, a₀₁ = b₀₁ = 0.

2.3. Magnetoconductivity

To demonstrate significance of the shape of the LLs for physical properties, we calculate the longitudinal magnetoconductivity using the linear-response theory as an example [40]. The magnetoconductivity σ_yy^col is given by

$$\begin{eqnarray}&&{\sigma }_{{yy}}^{{\rm{col}}}=\displaystyle \frac{\beta {e}^{2}}{{\rm{\Omega }}}\displaystyle \sum _{\zeta ,\zeta ^{\prime} }{\mathsf{f}}({\varepsilon }_{\zeta })\left[1-{\mathsf{f}}({\varepsilon }_{\zeta ^{\prime} })\right]{W}_{\zeta \zeta ^{\prime} }{\left({x}_{\zeta }-{x}_{\zeta ^{\prime} }\right)}^{2},\end{eqnarray} \tag{ 16 }$$

where β = 1/k_BT, ${\mathsf{f}}(\varepsilon )$ is the Fermi–Dirac distribution and W_ζζ' is the transition rate in the presence of disorder. The sum is over all quantum number $| \zeta \rangle =| s,n,{k}_{y}\rangle$ and the expectation value of the x coordinate is ${x}_{\zeta }=\langle \zeta | x| \zeta \rangle ={k}_{y}{l}^{2}$. If the impurity potential is short-ranged, the transition rate is expressed as

$$\begin{eqnarray}&&{W}_{\zeta \zeta ^{\prime} }=\displaystyle \frac{2\pi }{{\rm{\Omega }}}\displaystyle \sum _{{\bf{q}}}\gamma | {F}_{\zeta \zeta ^{\prime} }{| }^{2}\delta ({\varepsilon }_{\zeta }-{\varepsilon }_{\zeta ^{\prime} }){\delta }_{{k}_{y},{k}_{y}^{{\prime} }+{q}_{y}}.\end{eqnarray} \tag{ 17 }$$

Here, the form factor F_ζζ' is defined as ${F}_{\zeta \zeta ^{\prime} }=\langle \zeta | {{\bf{e}}}^{{\rm{i}}{\bf{q}}\cdot {\bf{r}}}| \zeta ^{\prime} \rangle$.

The magnetoconductivity of the LLs with s = ±1 is calculated by Biswas et al [33]. In their work, the scattering between eigenstates is only allowed with identical Landau indexes n = n', because of the existence of the term δ(ε_ζ − ε_ζ') in equation (17). However, this restriction is not valid for the zero-energy LL, and the contribution from the scattering between eigenstates with different Landau indexes must be calculated carefully and integrally.

Using the identity

$$\begin{eqnarray}\begin{array}{rcl}\int {\rm{d}}x{\phi }_{n}(x){{\rm{e}}}^{{{\rm{i}}{q}}_{x}x}{\phi }_{n^{\prime} }(x-{q}_{y}{l}^{2}) & = & {{\rm{e}}}^{-\tfrac{1}{4}({q}^{2}{l}^{2}-2{{\rm{i}}{q}}_{x}{q}_{y}{l}^{2})}\\ & & \times \ \sqrt{\displaystyle \frac{n^{\prime} !}{n!}}{\left(-\displaystyle \frac{{{\rm{i}}{q}}_{x}l+{q}_{y}l}{\sqrt{2}}\right)}^{n-n^{\prime} }{L}_{n^{\prime} }^{(n-n^{\prime} )}\left(\displaystyle \frac{{q}^{2}{l}^{2}}{2}\right),\end{array}\end{eqnarray} \tag{ 18 }$$

we obtain the square of form factor for the zero-energy LL

$$\begin{eqnarray}\begin{array}{rcl}| {F}_{{nn}^{\prime} }^{s=0}{| }^{2} & = & {\delta }_{{k}_{y},{k}_{y}^{{\prime} }+{q}_{y}}{{\rm{e}}}^{-u}{u}^{n-n^{\prime} }\\ & & \times \ {\left[{a}_{0n-2}{a}_{0n^{\prime} -2}\sqrt{\displaystyle \frac{(n^{\prime} -2)!}{(n-2)!}}{L}_{n^{\prime} -2}^{(n-n^{\prime} )}(u)+{b}_{0n}{b}_{0n^{\prime} }\sqrt{\displaystyle \frac{n^{\prime} !}{n!}}{L}_{n^{\prime} }^{(n-n^{\prime} )}(u)\right]}^{2},\end{array}\end{eqnarray} \tag{ 19 }$$

where u = q²l²/2 and L_n^(α)'s are generalized Laguerre polynomials. Generally, the disorder scattering gives a finite relaxation time τ, which are ignored in equation (17). A simple and heuristic way to account for this disorder effect is to replace the delta function in equation (17) by a Lorentzian. The relaxation time τ can be estimated from the width of level broadening Γ, τ = ℏ/Γ. Therefore, we replace the delta function by the Lorentzian $\tfrac{1}{\pi }\tfrac{{\rm{\Gamma }}}{{{\rm{\Gamma }}}^{2}+{({\varepsilon }_{\zeta }-{\varepsilon }_{\zeta ^{\prime} })}^{2}}$ [45, 46]. For intra-level scattering, we can further replace the delta function by ${\delta }_{{\varepsilon }_{\zeta }}$, ε_ζ'/πΓ[33, 8]. As a result of the above, we have

$$\begin{eqnarray}&&{\sigma }_{{yy}}^{{\rm{col}}}=\displaystyle \frac{{e}^{2}}{h}\displaystyle \frac{\beta \gamma }{\pi {\rm{\Gamma }}{l}^{2}}\displaystyle \sum _{{nn}^{\prime} }{\mathsf{f}}({\varepsilon }_{\zeta })\left[1-{\mathsf{f}}({\varepsilon }_{\zeta ^{\prime} })\right]\times \int u| {F}_{{nn}^{\prime} }^{s=0}{| }^{2}{\rm{d}}u.\end{eqnarray} \tag{ 20 }$$

Using the orthogonality of Laguerre polynomials

$$\begin{eqnarray}&&{\int }_{0}^{\infty }{u}^{\alpha }{{\rm{e}}}^{-u}{L}_{n}^{(\alpha )}(u){L}_{m}^{(\alpha )}(u){\rm{d}}u={\delta }_{{nm}}\displaystyle \frac{{\rm{\Gamma }}(n+\alpha +1)}{n!},\end{eqnarray} \tag{ 21 }$$

where Γ(n) denotes the gamma function, and the recurrence relation

$$\begin{eqnarray}&&{{uL}}_{n}^{(\alpha )}(u)=-(n+1){L}_{n+1}^{(\alpha )}(u)+(2n+\alpha +1){L}_{n}^{(\alpha )}(u)-(n+\alpha ){L}_{n-1}^{(\alpha )}(u),\end{eqnarray} \tag{ 22 }$$

we derive

$$\begin{eqnarray}&&\int u| {F}_{{nn}^{\prime} }^{s=0}{| }^{2}{\rm{d}}u={a}_{0n-2}^{2}{a}_{0n^{\prime} -2}^{2}(n+n^{\prime} -3)+{b}_{0n}^{2}{b}_{0n^{\prime} }^{2}(n+n^{\prime} +1).\end{eqnarray} \tag{ 23 }$$

At low temperatures, we can make the approximation $\beta {\mathsf{f}}(\varepsilon )\left[1-{\mathsf{f}}(\varepsilon )\right]\approx \delta (\varepsilon -{\varepsilon }_{F})$. Near the zero energy, the averaged DOS is defined as $\rho ({\varepsilon }_{F})=2\tfrac{A}{{\rm{\Omega }}}{\sum }_{n,{k}_{y}}\delta ({\varepsilon }_{F}-{\epsilon }_{0n})=\delta ({\varepsilon }_{F}-{\epsilon }_{0n})$ in the clean limit, where 2 comes from the valley degeneracy. The energies ε_ζ's in equation (20) are unperturbed eigenvalues calculated in the clean system, which are highly degenerate, especially for the zero-energy LL. However, the disorder scattering will shift the unperturbed eigenvalues and lift the degeneracy. The distribution of the shifted eigenenergies is determined by the DOS. At low temperature, therefore, we can approximately replace the delta function by the numerically calculated DOS in our calculations [47, 33]. Thus the magnetoconductivity can be expressed in terms of the DOS as

$$\begin{eqnarray}&&{\sigma }_{{yy}}^{{\rm{col}}}=\displaystyle \frac{{e}^{2}}{h}\displaystyle \frac{2\gamma }{{\rm{\Gamma }}}\rho ({\varepsilon }_{F})\times {I}_{F}({N}_{{\rm{\Phi }}}),\end{eqnarray} \tag{ 24 }$$

where

$$\begin{eqnarray}&&{I}_{F}({N}_{{\rm{\Phi }}})=\displaystyle \sum _{{nn}^{\prime} }^{{N}_{{\rm{\Phi }}}}\int | {F}_{{nn}^{\prime} }^{s=0}{| }^{2}u{\rm{d}}u.\end{eqnarray} \tag{ 25 }$$

According to equations (13) and (14), the large N_Φ limit of the sum in equation (25) is equal to ${N}_{{\rm{\Phi }}}^{4}/4$. Within this approximation, the magnetoconductivity can be expressed as

$$\begin{eqnarray}&&{\sigma }_{{yy}}^{{\rm{col}}}\approx \displaystyle \frac{{e}^{2}}{h}\displaystyle \frac{\gamma }{2{\rm{\Gamma }}}{\left(\displaystyle \frac{\pi {l}^{2}}{A}\right)}^{4}\rho ({\varepsilon }_{F}).\end{eqnarray} \tag{ 26 }$$

Before investigating the simulation results, it is necessary to figure out how many states there are in each LL, in view of the finite system size. For the conduction and valence bands LLs, the degeneracy is determined by the total magnetic flux going through the system, denoted by ${N}_{k}=2\times \tfrac{{\rm{\Omega }}}{2\pi {l}^{2}}$, where the prefactor 2 comes from the valley degrees of freedom. Because of the existence of the flat band, the degeneracy of the zero-energy LL is much greater and equals to N_ΦN_k = Ω/A, which is independent of the magnetic field strength.

The DOS of the LLs under Gaussian disorder and the uniform disorder are plotted in figure 2(a). The DOS of the zero-energy LL is much higher and wider than the other LLs since the area of the peak is proportional to the corresponding LL degeneracy. Analogous to the case in graphene, the conduction and valence bands LLs in the T₃ lattice are distributed unevenly due to the unique linear dispersion relation. As depicted in figure 2(b), the conduction and valence bands LLs at both sides are pushed away from the zero energy, due to the repulsion from the zero-energy LL. This blue shifts phenomenon is in sharp contrast with the disordered graphene [18]. That is because of the extraordinarily huge degeneracy of the zero-energy LL in the T₃ lattice.

Zoom In Zoom Out Reset image size

Next, we focus on the dependence of the zero-energy LL on the disorder and magnetic field. Due to the magnetic field independent degeneracy, it is expected to exhibit a distinct different scaling behavior from equation (1) [13, 48]. Within the SCBA (see in appendix), we find that the zero-energy LL in the T₃ lattice is semi-elliptical: $\rho (E)=\tfrac{1}{\pi \sqrt{\gamma }}\sqrt{2-\tfrac{{E}^{2}}{\gamma }}$, which is independent of the magnetic field strength. Compared to equation (1), we propose the following magnetic field independent scaling law to describe the disorder broadened zero-energy LL:

$$\begin{eqnarray}&&\rho (E)=\displaystyle \frac{1}{\sqrt{\gamma }}f\left(\displaystyle \frac{E}{\sqrt{\gamma }}\right),\end{eqnarray} \tag{ 27 }$$

where the form of the function f is closely related to disorder distribution.

In the presence of Gaussian disorder, the scaled numerical calculated DOS ${\gamma }^{\tfrac{1}{2}}\rho$ of the zero-energy LL as a function of the scaled energy $E/{\gamma }^{\tfrac{1}{2}}$ is plotted in figure 3(a). The symbols represent the results for different disorder strengths and magnetic fields. They collapse onto the same curve, well supporting the scaling law in equation (27). Furthermore, the function f with Gaussian disorder can be well fitted by a Gaussian function, which is similar to that in two-dimensional gas and graphene [11, 12, 18]. As for the uniform disorder, the scaling law still holds, although the form of f has changed. In figure 3(b), the shape of the zero-energy LL can be approximately described by a semi-elliptical form derived from SCBA (see appendix for details). Notice that the deviation between our numerical results and the SCBA prediction (especially at E = 0) may come from some nonperturbative effects. In view of the existence of the flat band, the successful nonperturbative supersymmetry method [11, 12] cannot be applied to study the shape of the zero-energy LL directly in our case. This challenging problem is unsolved and needs further investigation.

Zoom In Zoom Out Reset image size

To further check the scaling law, we plot the LL width as a function of disorder strength $\sqrt{\gamma }$ in figure 4, since equation (27) indicates that there exists a linear relation between them. The slopes of the fitting curves are 1.96 ± 0.03 and 2.41 ± 0.07 for the Gaussian and the uniform disorder, respectively. Additionally, although the SCBA fails to describe the exact shape of LL, the linear relation between LL width and disorder from SCBA is in good agreement with the numerical results with the uniform disorder: ${\rm{\Gamma }}=\sqrt{6\gamma }\approx 2.45\sqrt{\gamma }$.

Zoom In Zoom Out Reset image size

Finally, as an example of physical measurements related to the DOS of LL, we calculate the longitudinal magnetoconductivity σ_yy^col. The carrier density is defined through

$$\begin{eqnarray}&&{n}_{e}({\varepsilon }_{F})={\int }_{0}^{{\varepsilon }_{F}}\rho (\epsilon ){\rm{d}}\epsilon ,\end{eqnarray} \tag{ 28 }$$

and we can derive it by numerical integration. Substituting the scaling law equation (27) in equation (28), the charge density is expressed as

$$\begin{eqnarray}&&{n}_{e}={\boldsymbol{ \mathcal F }}\left(\displaystyle \frac{{\varepsilon }_{F}}{\sqrt{\gamma }}\right),\quad {\boldsymbol{ \mathcal F }}(x)={\int }_{0}^{x}f{\rm{d}}x^{\prime} .\end{eqnarray} \tag{ 29 }$$

Using the equations (26) and (28), we calculate the longitudinal magnetoconductivity ${\sigma }_{{yy}}^{{\rm{col}}}$ as a function of the carrier density n_e as plotted in figure 5. Furthermore, by substituting equation (29) into equation (26), the magnetoconductivity can be written as

$$\begin{eqnarray}&&{\sigma }_{{yy}}^{{\rm{col}}}=\displaystyle \frac{{e}^{2}}{h}\displaystyle \frac{\sqrt{\gamma }}{2{\rm{\Gamma }}}{\left(\displaystyle \frac{\pi {\hslash }}{{eBA}}\right)}^{4}f({{\boldsymbol{ \mathcal F }}}^{-1}({n}_{e})).\end{eqnarray} \tag{ 30 }$$

Thus, the magnetoconductivity is proportional to B⁻⁴ in the T₃ model, and this is contrary to the case in graphene, where the longitudinal magnetoconductivity saturates to a value ∼e²/h under high magnetic field [18, 49]. Since the linear relation between LL width and $\sqrt{\gamma }$ is independent of disorder strength, equation (30) shows that the magnetoconductivity is independent of disorder strength for fixed carrier density. This conclusion is also well verified by our simulations, as shown in figure 5, where the magnetoconductivities with different disorder strengths collapse onto the same curves.

Zoom In Zoom Out Reset image size

In summary, we numerically determine the exact shape of the zero-energy LL in a disordered T₃ lattice under strong magnetic field. Because of the small lattice spacing, although the magnetic flux per unit cell is small, the magnetic field is large and strong enough to avoid the mixing of disorder broadened LLs. We confirm that the existence of the flat band in the T₃ model leads to magnetic field independent scaling behavior of the zero-energy LL, in sharp contrast with graphene. Such unconventional scaling behavior results in the interesting fact that the magnetoconductivity is proportional to 1/B⁴. Since the T₃ model may be realized in electronic materials and cold-atom systems [34, 35, 26], we believe our predictions can be tested by thermodynamics and transport measurements in experiments.

This work was partially supported by National Key Research & Development Program of China (No. 2016YFA0200604 and 2017YFA0204904), and by the National Natural Science Foundation of China (Nos. 11874337, 11634011, and 21873088). Computational resources are provided by the Chinese Academy of Sciences, Shanghai and USTC Supercomputer Centers.

To understand the unconventional scaling behavior above, we calculate the DOS of zero-energy LL within SCBA [14]. Since the energy spacing between the LLs with different band indexes s is large enough, the contribution to the self-energy from weak disorder scattering between different band indexes s is negligible. Because the translational invariance is recovered after impurity averaging, the self-energy becomes diagonal with respect to the wave factor k_y. The diagonality with respect to LL index n can be obtained as below. The first order of self-energy is given by

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{n,n^{\prime} }^{(1)}(s,{k}_{y})=\overline{\langle s,n,{k}_{y}| {{VG}}^{0}V| s,n^{\prime} ,{k}_{y}\rangle },\end{eqnarray} \tag{ A.1 }$$

where ${G}^{0}={\left(E-{H}_{K}+{0}^{+}\right)}^{-1}$. Here the disorder is short ranged $\overline{V({\bf{r}})V({\bf{r}}^{\prime} )}=\gamma \delta ({\bf{r}}-{\bf{r}}^{\prime} )$, where γ = W²/12 for the uniform disorder and γ = σ² for the Gaussian disorder. After some algebra, we obtain

$$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Sigma }}}_{n,n^{\prime} }^{(1)}(s,k)={\delta }_{{nn}^{\prime} }\gamma \displaystyle \frac{{N}_{k}A}{{\rm{\Omega }}}\displaystyle \sum _{{s}_{1}{n}_{1}}{C}_{s}^{2}{C}_{{s}_{1}}^{2}({a}_{{sn}}^{2}{a}_{{s}_{1}{n}_{1}}^{2}+{s}^{2}{s}_{1}^{2}+{b}_{{sn}}^{2}{b}_{{s}_{1}{n}_{1}}^{2}){G}_{{s}_{1}{n}_{1}}^{0},\end{array}\end{eqnarray} \tag{ A.2 }$$

where ${G}_{{sn}}^{0}={\left(E-{\varepsilon }_{{sn}}+{0}^{+}\right)}^{-1}$. Therefore the first order of self-energy is diagonal with respect to s, n, k_y. Using the recursion relation of the self-energy

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{sn},{sn}^{\prime} }^{(N)}({k}_{y})=\overline{\langle s,n,{k}_{y}| {{VG}}^{0}{{\rm{\Sigma }}}^{(N-1)}{G}^{0}V| s,n^{\prime} ,{k}_{y}\rangle },\end{eqnarray} \tag{ A.3 }$$

and the diagonality of G⁰ with respect to s, n, k_y, we can generalized the diagonality in equation (A.1) to arbitrary order of the self-energy. In the end, the self-energy within SCBA takes the form

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{{snk}}_{y},s^{\prime} n^{\prime} {k}_{y}^{{\prime} }}={\delta }_{{ss}^{\prime} }{\delta }_{{nn}^{\prime} }{\delta }_{{k}_{y}{k}_{y}^{{\prime} }}{{\rm{\Sigma }}}_{{sn}},\end{eqnarray} \tag{ A.4 }$$

thus the Green function is also diagonal

$$\begin{eqnarray}&&{G}_{{{snk}}_{y},s^{\prime} n^{\prime} {k}_{y}^{{\prime} }}={\delta }_{{ss}^{\prime} }{\delta }_{{nn}^{\prime} }{\delta }_{{k}_{y}{k}_{y}^{{\prime} }}{G}_{{sn}}.\end{eqnarray} \tag{ A.5 }$$

After some algebra, the impurity averaged Green's function and self-energy can be written as

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{sn}=\gamma \displaystyle \frac{{N}_{k}A}{{\rm{\Omega }}}\times \displaystyle \sum _{{s}^{{\rm{{\prime} }}}{n}^{{\rm{{\prime} }}}}{G}_{{s}^{{\rm{{\prime} }}}{n}^{{\rm{{\prime} }}}}{C}_{s}^{2}{C}_{{s}^{{\rm{{\prime} }}}}^{2}\left({a}_{sn}^{2}{a}_{{s}^{{\rm{{\prime} }}}{n}^{{\rm{{\prime} }}}}^{2}+{s}^{2}s{{\rm{{\prime} }}}^{2}+{b}_{sn}^{2}{b}_{{s}^{{\rm{{\prime} }}}{n}^{{\rm{{\prime} }}}}^{2}\right).\end{eqnarray} \tag{ A.6 }$$

Notice that Σ_+n is identical to Σ_−n explicitly because a_−n = a_+n, and b_−n = b_+n, since the electron–hole symmetry is retained.

Generally, equation (A.6) can be solved numerically, while the analytical treatment need some approximation. For the zero-energy LL in the large LL index limit, we obtain

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{0n}=\gamma \displaystyle \frac{{N}_{k}A}{{\rm{\Omega }}}\displaystyle \sum _{n^{\prime} }({G}_{+n^{\prime} }/2+{G}_{-n^{\prime} }/2+{G}_{0n^{\prime} }),\end{eqnarray} \tag{ A.7 }$$

where the valley degree of freedom has been taken into account in terms of the number N_k. It will be a proper approximation to keep G_0n only since the LL mixing can be neglected in a strong magnetic field and the contribution from the zero-energy LL is dominant [50]. Furthermore, in view of the macroscopic degeneracy of the zero-energy LL, the self-energy dependence of the zero-energy LL upon the Landau index n can be ignored, so that Σ_0n ≈ Σ₀ and G_0n ≈ G₀ and the sum over LL indexes equals N_Φ. All in all, we obtain

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{0}=\displaystyle \frac{\gamma }{2}\displaystyle \frac{1}{E-{{\rm{\Sigma }}}_{0}}.\end{eqnarray} \tag{ A.8 }$$

Notice that disorder strength is the only surviving physical parameter, thus the DOS is expected to be independent of magnetic field. The analytical solution of the zero-energy LL is

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{0}(E+{\rm{i}}0)=\displaystyle \frac{1}{2}(E-{\rm{i}}\sqrt{2\gamma -{E}^{2}}),\end{eqnarray} \tag{ A.9 }$$

and the DOS is expressed as

$$\begin{eqnarray}&&\rho (E)=\displaystyle \frac{1}{\pi \sqrt{\gamma }}\sqrt{2-\displaystyle \frac{{E}^{2}}{\gamma }}.\end{eqnarray} \tag{ A.10 }$$

In conclusion, a semi-elliptical DOS of the zero-energy LL is suggested within SCBA no matter what kind of disorder distribution is considered. On the other hand, equation (A.10) gives a linear relation between LL width and disorder ${\rm{\Gamma }}=\sqrt{6\gamma }\approx 2.45\sqrt{\gamma }$.