Quantum valley Hall effect, orbital magnetism, and anomalous Hall effect in twisted multilayer graphene systems

TL;DR

This paper investigates the topological electronic properties of twisted multilayer graphene systems, revealing how stacking and twisting influence valley Chern numbers, orbital magnetism, and potential anomalous Hall effects.

Contribution

It introduces a comprehensive analysis of the topological and magnetic properties of twisted multilayer graphene, highlighting the dependence of valley Chern numbers on stacking chirality and layer number.

Findings

Existence of two low-energy flat bands per valley in twisted multilayer graphene.

Valley Chern numbers depend on layer number and stacking chirality.

Potential for orbital ferromagnetism and anomalous Hall effect due to valley-contrasting orbital magnetizations.

Abstract

We study the electronic structures and topological properties of $(M+N)$-layer twisted graphene systems. We consider the generic situation that $N$-layer graphene is placed on top of the other $M$-layer graphene, and is twisted with respect to each other by an angle $\theta$. In such twisted multilayer graphene (TMG) systems, we find that there exists two low-energy flat bands for each valley emerging from the interface between the $M$ layers and the $N$ layers. These two low-energy bands in the TMG system possess valley Chern numbers that are dependent on both the number of layers and the stacking chiralities. In particular, when the stacking chiralities of the $M$ layers and $N$ layers are opposite, the total Chern number of the two low-energy bands for each valley equals to $\pm(M+N-2)$ (per spin). If the stacking chiralities of the $M$ layers and the $N$ layers are the same, then the total Chern number of the two low-energy bands for each valley is $\pm(M-N)$ (per spin). The valley Chern numbers of the low-energy bands are associated with large, valley-contrasting orbital magnetizations, suggesting the possible existence of orbital ferromagnetism and anomalous Hall effect once the valley degeneracy is lifted either externally by a weak magnetic field or internally by Coulomb interaction through spontaneous symmetry breaking.

Figures

Tables

Table 1. Table 1: The dependence of the Chern numbers on vertical electric potential for (2+2)-layer TMG

$V_{⟂}$	-60	-48	-36	-24	-12	0	12	24	36	48	60
$C_{1++}^K$	+1	+1	0	-2	-3	-	+3	+2	0	-1	-1
$C_{2++}^K$	+1	+1	+1	+2	+3	-	-3	-2	-1	-1	-1
$C_{1+-}^K$	0	0	+1	+1	0	0	0	+1	+1	0	0
$C_{2+-}^K$	+1	+1	+1	+1	+2	+2	+2	+1	+1	+1	+1

Table 2. Table 2: The dependence of the Chern numbers on vertical electric potential for (2+1)-layer TMG

$V_{⟂}$	-60	-48	-36	-24	-12	0	12	24	36	48	60
$C_{1+}^K$	0	-1	-1	-1	+2	+1	+2	+2	+2	+2	+2
$C_{2+}^K$	+2	+2	+2	+2	-1	0	-1	-1	-1	-1	-1

Equations

$$d_z(\mathbf{r})=d_0+2d_1\sum _{j=1}^3\cos ({\mathbf{g}}_j\cdot \mathbf{r}),$$

$$H_{\alpha ,{\alpha }^{\prime }}^K(M+N)=\left(\begin{array}{cc}{H}_{\alpha }^K(M)& \mathbb{U}\\ {\mathbb{U}}^{†}& H_{{\alpha }^{\prime }}^K(N)\end{array}\right),$$

$$H_{\alpha }^K(M)=\left(\begin{array}{ccccc}{h}_0(\mathbf{k})& h_{\alpha }& 0& 0& ...\\ h_{\alpha }^{†}& h_0(\mathbf{k})& h_{\alpha }& 0& ...\\ 0& h_{\alpha }^{†}& h_0(\mathbf{k})& h_{\alpha }& ...\\ & & & ...& \end{array}\right),$$

$$h_{+}=\left(\begin{array}{cc}0& 0\\ t_{\perp }& 0\end{array}\right),$$

$$\mathbb{U}=\left(\begin{array}{ccc}0& ...& 0\\ \text{mathmakebox}[\text{widthof}0][c]⋮& ...& 0\\ U(\mathbf{r})e^{-i\Delta \mathbf{K}\cdot \mathbf{r}}& ...& 0\end{array}\right),$$

$$U(\mathbf{r})=\left(\begin{array}{cc}{u}_{0}g(\mathbf{r})& u_0^{\prime }g(\mathbf{r}-{\mathbf{r}}_{AB})\\ u_0^{\prime }g(\mathbf{r}+{\mathbf{r}}_{AB})& u_{0}g(\mathbf{r})\end{array}\right),$$

$$h_{+}=\left(\begin{array}{cc}{t}_{2}f(\mathbf{k})& t_2{f}^{\ast }(\mathbf{k})\\ t_{\perp }-3t_3& t_{2}f(\mathbf{k})\end{array}\right),$$

$$C_{\alpha ,{\alpha }^{\prime }}^K=+(\alpha (M-1)-{\alpha }^{\prime }(N-1)),$$

$$C_{\alpha ,{\alpha }^{\prime }}^{K^{\prime }}=-(\alpha (M-1)-{\alpha }^{\prime }(N-1)),$$

$$t(\mathbf{d})=V_{\sigma }(\frac{\mathbf{d}\cdot \hat{\mathbf{z}}}{d}{)}^2+V_{\pi }(1-(\frac{\mathbf{d}\cdot \hat{\mathbf{z}}}{d}{)}^2)$$

$$M_z=\frac{e}{2\hslash (2\pi {)}^2}$$

$$\times (H_{\mathbf{k}}+{ϵ}_{n\mathbf{k}}-2\mu )|{\partial }_{\mathbf{k}}{u}_{n\mathbf{k}}⟩\},$$

$$H_{+,+}^K(2+2)=\left(\begin{array}{cccc}{h}_0(\hat{\mathbf{k}})& h_{+}& 0& 0\\ h_{+}^{†}& h_0(\hat{\mathbf{k}})& U(\mathbf{r})& 0\\ 0& U^{†}(\mathbf{r})& h_0(\hat{\mathbf{k}})& h_{+}\\ 0& 0& h_{+}^{†}& h_0(\hat{\mathbf{k}})\end{array}\right),$$

$$\displaystyle|\widetilde{\psi}_{Ms}\rangle=\frac{1}{2}\Big{(}\,|{\psi}_{Ms}\rangle+i|{\psi}_{(M+1)s}\rangle\,\Big{)}\;,$$

$$\displaystyle|\widetilde{\psi}_{(M+1)s}\rangle=\frac{1}{2}\Big{(}\,|{\psi}_{Ms}\rangle-i|{\psi}_{(M+1)s}\rangle\,\Big{)}\;,$$

$$|{\tilde{\psi }}_{ls}⟩=|{\psi }_{ls}⟩,\text{ if l=M,M+1}$$

$${\tilde{H}}_{+,+}^K(2+2)=\left(\begin{array}{cccc}{h}_0(\hat{\mathbf{k}})& \frac{h_{+}}{\sqrt{2}}& \frac{h_{+}}{\sqrt{2}}& 0\\ \frac{h_{+}^{†}}{\sqrt{2}}& h_0(\hat{\mathbf{k}}-\frac{e}{\hslash }\mathbf{A})& -3i{u}_0& -i\frac{h_{+}}{\sqrt{2}}\\ \frac{h_{+}^{†}}{\sqrt{2}}& 3i{u}_0& h_0(\hat{\mathbf{k}}+\frac{e}{\hslash }\mathbf{A})& i\frac{h_{+}}{\sqrt{2}}\\ 0& -i\frac{h_{+}^{†}}{\sqrt{2}}& i\frac{h_{+}^{†}}{\sqrt{2}}& h_0(\hat{\mathbf{k}})\end{array}\right)$$

$${\tilde{H}}_{+,+}^K(2+2)\approx \left(\begin{array}{cc}{h}_b& 0\\ 0& h_t\end{array}\right).$$

$$h_b=\left(\begin{array}{ccc}0& \hslash v_F({\hat{k}}_x+i{\hat{k}}_y)& 0\\ \hslash v_F({\hat{k}}_x-i{\hat{k}}_y)& 0& t_{\perp }\eta (k)\\ 0& t_{\perp }{\eta }^{\ast }(k)& 0\end{array}\right),$$

$$h_t=\left(\begin{array}{ccc}0& i{t}_{\perp }\eta (k)& 0\\ -i{t}_{\perp }{\eta }^{\ast }(k)& 0& \hslash v_F({\hat{k}}_x+i{\hat{k}}_y)\\ 0& \hslash v_F({\hat{k}}_x-i{\hat{k}}_y)& 0\end{array}\right),$$

$$H_{\text{eff}}^{\alpha }(M-1)$$

$$=$$

$$H_{K_s}^{\alpha }(\mathbf{k})=\left(\begin{array}{ccc}0& \frac{(k_x+i\alpha k_y{)}^{M-1}}{m_{M-1}}& 0\\ \frac{(k_x-i\alpha k_y{)}^{M-1}}{m_{M-1}}& 0& t_{\text{eff}}\\ 0& t_{\text{eff}}& 0\end{array}\right),$$

$$m_{M-1}=\frac{t_{\perp }^{M-2}}{(-1{)}^{M-2}(\hslash v_F{)}^{M-1}},$$

$$H_{K_s^{\prime }}^{{\alpha }^{\prime }}(\mathbf{k})=\left(\begin{array}{ccc}0& \frac{(k_x+i{\alpha }^{\prime }{k}_y{)}^{N-1}}{m_{N-1}}& t_{\text{eff}}^{\ast }\\ \frac{(k_x-i{\alpha }^{\prime }{k}_y{)}^{N-1}}{m_{N-1}}& 0& 0\\ t_{\text{eff}}& 0& 0\end{array}\right),$$

$${ϵ}_{1\mathbf{k}}(M-1)=-\sqrt{m_{M-1}^2{k}^{2M-2}+|t_{\text{eff}}{|}^2}$$

$${ϵ}_{2\mathbf{k}}=0$$

$${ϵ}_{3\mathbf{k}}(M-1)=\sqrt{m_{M-1}^2{k}^{2M-2}+|t_{\text{eff}}{|}^2}.$$

$$\displaystyle|\psi_{1\mathbf{k}}^{M-1}\rangle=\Big{(}\,\frac{(k_{x}+i\alpha k_{y})^{M-1}}{\sqrt{2}g_{M-1}(|\mathbf{k}|)}\,,\,-\frac{1}{\sqrt{2}}\,,\,\frac{\widetilde{t}_{\textrm{eff}}}{\sqrt{2}g_{M-1}(|\mathbf{k}|)}\,\Big{)}^{T}\;$$

$$\displaystyle|\psi_{2\mathbf{k}}^{M-1}\rangle=\Big{(}\,\frac{-\widetilde{t}_{\textrm{eff}}\,e^{i\alpha\theta_{\mathbf{k}}(M-1)}}{g_{M-1}(\mathbf{k})}\,,\,0\,,\,\frac{|\mathbf{k}|^{(M-1)}}{g_{M-1}(\mathbf{k})}\,\Big{)}^{T}\;$$

$$\displaystyle|\psi_{3\mathbf{k}}^{M-1}\rangle=\Big{(}\,\frac{(k_{x}+i\alpha k_{y})^{M-1}}{\sqrt{2}g_{M-1}(|\mathbf{k}|)}\,,\,\frac{1}{\sqrt{2}}\,,\,\frac{\widetilde{t}_{\textrm{eff}}}{\sqrt{2}g_{M-1}(|\mathbf{k}|)}\,\Big{)}^{T}$$

$$\displaystyle|\psi_{1\mathbf{k}}^{N-1}\rangle=\Big{(}\,-\frac{1}{\sqrt{2}}\,,\,\frac{(k_{x}-i\alpha^{\prime}k_{y})^{N-1}}{\sqrt{2}\,g_{N-1}(|\mathbf{k}|)}\,,\,\frac{\widetilde{t}_{\textrm{eff}}}{\sqrt{2}\,g_{N-1}(|\mathbf{k}|)}\,\Big{)}^{T}\;$$

$$\displaystyle|\psi_{2\mathbf{k}}^{N-1}\rangle=\Big{(}\,0\,,\,-\frac{\widetilde{t}_{\textrm{eff}}\,e^{-i\alpha^{\prime}\theta_{\mathbf{k}}(N-1)}}{g_{N-1}(\mathbf{k})}\,,\,\frac{|\mathbf{k}|^{(N-1)}}{g_{N-1}(\mathbf{k})}\,\Big{)}^{T}\;$$

$$\displaystyle|\psi_{3\mathbf{k}}^{N-1}\rangle=\Big{(}\,\frac{1}{\sqrt{2}}\,,\,\frac{(k_{x}-i\alpha^{\prime}k_{y})^{N-1}}{\sqrt{2}\,g_{N-1}(|\mathbf{k}|)}\,,\,\frac{\widetilde{t}_{\textrm{eff}}}{\sqrt{2}\,g_{N-1}(|\mathbf{k}|)}\,\Big{)}^{T}\;.$$

$$A_{1{\theta }_{\mathbf{k}}}^{M-1}=i⟨{\psi }_{1\mathbf{k}}^{M-1}|{\partial }_{{\theta }_{\mathbf{k}}}{\psi }_{1\mathbf{k}}^{M-1}⟩=-\frac{\alpha (M-1)|\mathbf{k}{|}^{2M-2}}{2{\tilde{t}}_{\text{eff}}^2+2|\mathbf{k}{|}^{2M-2}},$$

$$A_{1{\theta }_{\mathbf{k}}}^{N-1}=i⟨{\psi }_{1\mathbf{k}}^{N-1}|{\partial }_{{\theta }_{\mathbf{k}}}{\psi }_{1\mathbf{k}}^{N-1}⟩=\frac{{\alpha }^{\prime }(N-1)|\mathbf{k}{|}^{2N-2}}{2{\tilde{t}}_{\text{eff}}^2+2|\mathbf{k}{|}^{2N-2}},$$

Full text

Twisted bilayer graphene (TBG) has drawn significant attention recently due to the observations of the correlated insulating phases Cao et al. (2018a); Sharpe et al. (2019); Choi et al. (2019); Kerelsky et al. (2018); Codecido et al. (2019) and unconventional superconductivity Cao et al. (2018b); Codecido et al. (2019). At small twist angles, the low-energy states of TBG are characterized by four low-energy bands contributed by the two nearly decoupled monolayer valleys Lopes dos Santos et al. (2007); Bistritzer and MacDonald (2011). Around the “magic angles”, the bandwidths of the four low-energy bands become vanishingly small, and these nearly flat bands are believed to be responsible for most of those exotic properties observed in TBG. Numerous theories have been proposed to understand the electronic structures Po et al. (2018a); Yuan and Fu (2018); Koshino et al. (2018); Kang and Vafek (2018a); Song et al. (2018); Po et al. (2018b); Tarnopolsky et al. (2019); Pal et al. (2018); Liu et al. (2018a); Lian et al. (2018a); Zhang et al. (2019a), the correlated insulating phase Po et al. (2018a); Sboychakov et al. (2018); Isobe et al. (2018); Xu et al. (2018); Huang et al. (2018); Liu et al. (2018b); Rademaker and Mellado (2018); Venderbos and Fernandes (2018); Kang and Vafek (2018b); Xie and MacDonald (2018); Jian and Xu (2018); Bultinck et al. (2019), and the mechanism of superconductivity Xu and Balents (2018); Po et al. (2018a); Isobe et al. (2018); Wu et al. (2018a, b); Lian et al. (2018b); Huang et al. (2019); Liu et al. (2018b); Venderbos and Fernandes (2018); Kozii et al. (2018); Wu (2018); Roy and Jurivcić (2019).

On the other hand, interesting topological features have already emerged in the electronic structure of TBG. It has been shown that the four low-energy bands are topologically nontrivial in the sense that they are characterized by odd windings of Wilson loops Song et al. (2018); Ahn et al. (2018); Liu et al. (2018a), which is an example of the fragile topology Po et al. (2018b). The four flat bands have been further proposed to be equivalent to the zeroth pseudo Landau levels (LLs) with opposite Chern numbers and sublattice polarizations Liu et al. (2018a), which is the origin of the nontrivial band topology in the TBG system.

Moreover, recently unconventional ferromagnetic superconductivity and correlated insulating phase have been observed in twisted double bilayer graphene Liu et al. (2019); Cao et al. (2019). It implies that the low-energy flat bands, which are believed to be responsible for the correlated physics in TBG, may also exist in the twisted double bilayer graphene system. A recent theoretical study indeed revealed the presence of flat bands in twisted double bilayer graphene Chebrolu et al. (2019). Motivated by these works, in this paper we study the electronic structures and topological properties of twisted multilayer graphene (TMG). In particular, we consider the most generic situation that the $N$-layer chirally stacked graphene is placed on top of the other $M$-layer chirally stacked graphene, and they are twisted with respect to each other by a non-vanishing angle $\theta$, as schematically shown in Fig. 1(a) (for the case of $M\!=\!2,N\!=\!2$). In such a ($M+N$)-layer TMG system, we propose that there always exists two low-energy bands (for each valley), and that the bandwidths of the two low-energy bands become vanishingly small at the magic angles of twisted bilayer graphene (TBG) for arbitrary numbers of layers $M$ and $N$. The flat bands in the TMG system can be interpreted from the pseudo LL representation of TBG Liu et al. (2018a), and is protected by an approximate chiral symmetry in chiral graphene multilayers.

Moreover, we also find that there is a Chern-number hierarchy in the $(M+N)$-layer TMG system. In particular, when the stacking chiralities of the $M$ layers and $N$ layers are the same, the total Chern number of the two low-energy flat bands for each monolayer valley equals to $\pm(M-N)$ for each spin species 111Unless otherwise specified, in the rest of the paper the Chern number refers to the Chern number for each physical spin speicies.. On the other hand, if the stacking chiralities of the $M$ layers and the $N$ layers are opposite, then the total Chern number of the two low-energy bands for each valley is $\pm(M+N-2)$. The valley Chern numbers can be further tuned by an external electric field, leading to gate-tunable quantum valley Hall effect.

The nonzero valley Chern numbers of the low-energy flat bands are characterized by large and valley-contrasting orbital magnetizations. With the presence of an external magnetic field or the spontaneous symmetry breaking induced by the Coulomb interactions, the valley degeneracy is expected to be broken, and a valley-polarized (quantum) anomalous Hall state may be realized. The valley polarized state are associated with chiral current loops, which generate local magnetic fields peaked at the $AA$ region. The local magnetic fields generated by the chiral current loops may be a robust experimental signature for the nonzero valley Chern number and the valley polarized state in the TMG system. The flat bands at the universal magic angles, together with the Chern-number hierarchy and orbital magnetism, make the TMG systems a unique platform to study strongly correlated physics with nontrivial band topology, and may have significant implications on the observed ferromagnetic superconductivity and correlated insulating phase in twisted double bilayer graphene Liu et al. (2019); Cao et al. (2019).

I Electronic structures of the twisted multilayer graphene systems

I.1 The lattice structures

We consider the most generic case of chirally stacked twisted multilayer graphene, i.e., we place $N$ chiral graphene multilayers on top of $M$ chiral graphene multilayers, and twist them with respect to each other by an angle $\theta$. This is schematically shown in Fig. 1(a) for the case of $M\!=\!2,N=2\!$. Similar to the case of TBG, commensurate moiré supercells are formed when the twist angle $\theta(m)$ obeys the condition $\cos{\theta(m)}=(3m^{2}+3m+1/2)/(3m^{2}+3m+1)$ Lopes dos Santos et al. (2012), where $m$ is a positive integer. The lattice vectors of the moiré superlattice are expressed as $\mathbf{t}_{1}=(-\sqrt{3}L_{s}/2,L_{s}/2)$, and $\mathbf{t}_{2}=(0,L_{s})$, where $L_{s}=|\mathbf{t}_{1}|=a/(2\sin{(\theta/2)})$ is the size of the moiré supercell, and $a=2.46\,$Å is the lattice constant of graphene. In TBG it is well known that there are atomic corrugations, i.e., the variation of interlayer distances on the moiré length scale. In particular, in the $AB(BA)$ region of TBG, the interlayer distance $d_{AB}\!\approx\!3.35\,$Å  while in the $AA$-stacked region the interlayer distance $d_{AA}\!\approx\!3.6\,$Å Lee et al. (2008). Such atomic corrugations may be modeled as Koshino et al. (2018)

[TABLE]

where $\mathbf{g}_{1}$, $\mathbf{g}_{2}$, and $\mathbf{g}_{3}\!=\!\mathbf{g}_{1}+\mathbf{g}_{2}$ are the three reciprocal lattice vectors of the moire supercell. We take $d_{0}=3.433\,$Å and $d_{1}=0.0278\,$Å in order to reproduce the interlayer distances in $AA$- and $AB$-stacked bilayer graphene. In this paper, the atomic corrugations of the two twisted layers at the interface (between the $M$ layers and the $N$ layers) is also be modeled by Eq. (1). On the other hand, the interlayer distances within the untwisted $M$ layers and the untwisted $N$ layers are set to the interlayer distance of Bernal bilayer graphene $d_{AB}\!=\!3.35\,$Å. At a small twist angle $\theta$, the Brillouin zone (BZ) of the moiré supercell has been significantly reduced compared with those of the untwisted multilayers as shown in Fig. 1(b).

I.2 The effective Hamiltonian

The low-energy effective Hamiltonian of the twisted $(M+N)$-layer TMG of the $K$ valley is expressed as

[TABLE]

where $H^{K}_{\alpha}(M)$ and $H^{+}_{\alpha^{\prime}}(N)$ are the effective Hamiltonians for the $M$-layer and $N$-layer graphene with stacking chiralities $\alpha,\alpha^{\prime}=+/-$. In particular,

[TABLE]

where $h_{0}(\mathbf{k})\!=\!-\hbar v_{F}(\mathbf{k}-\mathbf{K}_{M})\cdot\mathbf{\sigma}$ stands for the low-energy effective Hamiltonian for monolayer graphene near the Dirac point $\mathbf{K}_{M}$, and $h_{\alpha}$ is the interlayer hopping, with

[TABLE]

and $h_{-}=h_{+}^{\dagger}$.

The off-diagonal term $\mathbb{U}$ represents the coupling between the twisted $M$ layers and $N$ layers. Here we assume that there is only the nearest neighbor interlayer coupling, i.e., the topmost layer of the $M$-layer graphene is only coupled with the bottom-most layer of the $N$-layer graphene, thus

[TABLE]

where the $2\!\times\!2$ matrix $U$ describes the tunneling between the Dirac states of the twisted bilayers Bistritzer and MacDonald (2011); Koshino et al. (2018)

[TABLE]

where $\mathbf{r}_{AB}\!=\!(\sqrt{3}L_{s}/3,0)$, $u_{0}^{\prime}$ and $u_{0}$ denote the intersublattice and intrasublattice interlayer tunneling amplitudes, with $u_{0}^{\prime}\!\approx\!0.098\,$eV, and $u_{0}\!\approx\!0.078$ eV Koshino et al. (2018). $u_{0}\!$ is smaller than $u_{0}^{\prime}$ due to the effects of atomic corrugations Koshino et al. (2018); Liu et al. (2018a). $\Delta\mathbf{K}=\mathbf{K}_{N}-\mathbf{K}_{M}=(0,4\pi/3L_{s})$ is the shift between the Dirac points of the $N$ layers and the $M$ layers. The phase factor $g(\mathbf{r})$ is defined as $g(\mathbf{r})=\sum_{j=1}^{3}e^{i\mathbf{q}_{j}\cdot\mathbf{r}}$, with $\mathbf{q}_{1}=(0,4\pi/3L_{s})$, $\mathbf{q}_{2}=(-2\pi/\sqrt{3}L_{s},-2\pi/3L_{s})$, and $\mathbf{q}_{3}=(2\pi/\sqrt{3}L_{s},-2\pi/3L_{s})$. It worth to note that Eq. (2) is the effective Hamiltonian for the $K$ valley. The Hamiltonian for the $K^{\prime}$ valley is readily obtained by applying a time-reversal operation to $H^{K}_{\alpha,\alpha^{\prime}}(M+N)$.

I.3 The emergence of two flat bands and the universal magic angles

We continue to study the electronic structures of the $(M+N)$-layer TMG systems using the effective Hamiltonian given by Eq. (2). The bandstructures for ($M\!=\!2$, $N\!=\!2$), ($M\!=\!2$, $N\!=\!4$), $(M\!=\!3,N\!=\!3)$, and $(M\!=\!3,N\!=\!4)$ at the first magic angle of TBG $\theta\!=\!1.05^{\circ}$ with the same stacking chiralities ($\alpha\!=\!\alpha^{\prime}=+$) are shown in Fig. 2(a)-(d) respectively 222The bandstructures with opposite stacking chiralities are very similar to those shown in Fig. 2.. Clearly there are two low-energy flat bands marked by the red lines that are separated from the other bands. The two low-energy bands are almost exactly flat at $\theta\!=\!1.05^{\circ}$ for all these TMG systems with different layers, indicating that the magic angle of TBG is universal for the TMG systems regardless the number of layers. It turns out that the two flat bands in TMG originates from the twisted bilayer at the interface, and they remain flat even after being coupled with the other graphene layers due to an (approximate) chiral symmetry of Eq. (2). More details about the origin of the flat bands in the TMG systems can be found in Appendix  A.

In realistic situations there are also further neighbor interlayer hoppings in graphene multilayers, which would break the chiral symmetry of the effective Hamiltonian in Eq. (2), and the flat bands shown in Fig. 2 would become more dispersive. In order to test the robustness of the flat bands, we have included all the second-neighbor and third-neighbor interlayer hoppings with intersite distances equal to $\sqrt{a^{2}/3+d_{AB}^{2}}$ and $\sqrt{a^{2}+d_{AB}^{2}}$ respectively ($d_{AB}\!\approx\!3.35\,$Åis the interlayer distance), and their amplitudes are denoted by $t_{2}$ and $t_{3}$. After including these terms, the interlayer hopping term with $+$ stacking chirality becomes

[TABLE]

where in $t_{2}\!=\!0.21$ eV, $t_{3}\!\approx\!0.05\,$eV are extracted from the Slater-Koster hopping parameters (see Eq. (9)). In order to be consistent with the choice of $t_{2}$ and $t_{3}$, we set $t_{\perp}\!=\!0.48\,$eV, which is also from the Slater-Koster formula (Eq. (9). The phase factor $f(\mathbf{k})\!=\!(e^{-i\sqrt{3}ak_{y}/3}+e^{i(k_{x}a/2+\sqrt{3}ak_{y}/6})+e^{i(-k_{x}a/2+\sqrt{3}ak_{y}/6})$. The interlayer hopping with $-$ stacking chirality $h_{-}\!=\!h_{+}^{\dagger}$. The bandstructures of (2+2)-layer TMG at $\theta\!=\!1.05^{\circ}$ with the new interlayer hopping term Eq. (7) are shown in Fig. 3, where (a) and (b) denote the cases with the same and opposite stacking chiralities respectively. Clearly the two low-energy bands marked by the red lines become more dispersive due to the presence of the further-neighbor interlayer hoppings, but the bandwidths are still small $\sim\!10-15$ meV, and the two low-energy bands are still separated from the high-energy bands.

II The Chern-number hierarchy and quantum valley Hall effect

II.1 The Chern-number hierarchy

The flat bands at the universal magic angle make the TMG systems a perfect platform to study the strongly correlated physics. In addition to the flat bands and the universal magic angles, the low-energy bands in the TMG systems also exhibit unusual topological properties with non-vanishing valley Chern numbers. To be specific, when the stacking chiralities of the $M$ layers and the $N$ layers are the same , the total Chern number of the two low-energy bands for each monolayer valley equals to $\pm(M-N)$. On the other hand, if the stacking chiralities of the $M$ layers and the $N$ layers are opposite, then the total Chern number of the two flat bands for each valley equals to $\pm(M+N-2)$. Such a Chern-number hierarchy is more concisely summerized in the following equation

[TABLE]

where $C_{\alpha,\alpha^{\prime}}^{K}$ ($C_{\alpha,\alpha^{\prime}}^{K^{\prime}}$) denotes the total Chern number of the two low-energy flat bands for the $K$ ($K^{\prime}$) valley, and the subscripts $\alpha,\alpha^{\prime}=\pm$ represent the stacking chiralities of the $M$ layers and $N$ layers. We would like to emphasize that the total Chern number of the two flat bands (per valley per spin) is a more robust quantity than the Chern number of each individual flat band. This is because the former is protected by the energy gaps between the two flat bands and the other high-energy bands, while the latter is crucially dependent on how the gap between the two flat bands is opened up.

In order to understand the Chern-number hierarchy of Eq. (8), we first divide the $(M+N)$-layer TMG system into three mutually decoupled subsystems: the TBG at the interface, the $(M-1)$ graphene monolayers below the interface TBG, and the $(N-1)$ graphene monolayers above the interface TBG, which are schematically shown in Fig. 5(a). We introduce a scaling parameter $0\!\leq\!\lambda\!\leq\!1$, and let the coupling strength between the three subsystems $t_{\perp}\!\to\!\lambda t_{\perp}$. We adiabatically turn on the coupling between the three subsystems by increasing $\lambda$ from 0 to 1, then inspect the evolution of the bandstructures of the TMG system.

In Fig. 5(b) we show the bandstructure of $(3+2)$-layer TMG (of the $K$ valley) at $\theta\!=\!1.05^{\circ}$ with the scaling parameter $\lambda\!=\!0$. When $\lambda\!=\!0$, the magic-angle TBG at the interface would give rise to two flat bands with total Chern number 0 as marked by the red lines in Fig. 5(a). The $(M-1)$ graphene monolayers below the TBG interface would contribute two low-energy bands with dispersions $\sim\!\pm|\mathbf{k}|^{M-1}$ around $K_{s}$ 333The Dirac points of the $M$ graphene layers $K_{M}$ and of the $N$ graphene monolayers $K_{N}$ are respectively mapped to $K_{s}$ and $K_{s}^{\prime}$ points of the moiré Brillouin zone, see Fig. 1(b) . Similarly, the $(N-1)$ layers above the TBG interface would contribute two low-energy bands with dispersions $\sim\!\pm|\mathbf{k}|^{N-1}$ around $K_{s}^{\prime}$. Since we have considered the case $M\!=\!3$ and $N\!=\!2$, there are quadratic band touching at $K_{s}$ and linear band touching at $K_{s}^{\prime}$ in Fig. 5(b). If $\lambda$ becomes nonzero, gaps would be opened up at $K_{s}$ and $K_{s}^{\prime}$ for the bands with $\pm|\mathbf{k}|^{M-1}$ and $\pm|\mathbf{k}|^{N-1}$ dispersions as shown in Fig. 5(c) for $\lambda\!=\!0.05$. As a result, the loop integral for the Berry’s connection of the conduction and valance bands around a loop enclosing the $K_{s}$ point would acquire the same Berry phase of $-\alpha(M-1)\pi$ ($\alpha$ is the stacking chirality of the $(M-1)$ layers), and contribute $-\alpha(M-1)/2$ to the total Chern number respectively, which adds up to $-\alpha(M-1)$ (see Appendix B). On the other hand, the conduction and valence bands around $K_{s}^{\prime}$ point contributed by the $(N-1)$ layers above the interface would acquire the same Berry phase of $\alpha^{\prime}(N-1)\pi$ ($\alpha^{\prime}$ is the stacking chirality of the $(N-1)$ layers), with the total Chern number of $\alpha^{\prime}(N-1)$ (see Appendix B). It is well known that the total Chern number of the bands from the $(M-1)$ and $(N-1)$ layers must cancel that of the two flat bands from the interface TBG, it follows that the total Chern number of the two flat bands for the $K$ valley equals to $\alpha(M-1)-\alpha^{\prime}(N-1)$. As $\lambda$ is further increased, the conduction and valence bands from the $(M-1)$ and $(N-1)$ layers are further pushed to high energies as shown in Fig. 5(d) for $\lambda\!=\!0.5$, and the Chern number of the two flat bands would remain unchanged. Thus Eq. (8) has been proved. We refer the readers to Appendix B for more details.

Eq. (8) has been numerically verified using the effective Hamiltonian of TMG shown in Eq. (2). In particular, in Fig. 4(a) we plot the Wilson-loop eigenvalues (denoted as $w(k)$) of the $(2+2)$ TMG ($M\!=\!2$, $N\!=\!2$) at the first magic angle with the same stacking chirality ($\alpha\!=\!\alpha^{\prime}\!=\!+$). The red diamonds and blue circles represent the Wilson loops of the $K$ and $K^{\prime}$ valleys respectively. As clearly shown in the figure, for each valley the total Chern number of the two flat bands vanishes. In Fig. 4(b) we plot the Wilson loops of the $(2+2)$ TMG at the first magic angle, but with opposite stacking chiralities ($\alpha\!=\!+,\alpha^{\prime}\!=\!-$). It is clearly seen that for the $K$ valley (blue circles) the two Wilson loops carry the same Chern number $+1$, giving rise to a total Chern number of $+2$ for the $K$ valley ($-2$ for the $K^{\prime}$ valley), which is consistent with Eq. (8). In Fig. 4(c) and (d) we plot the Wilson-loop eigenvalues for the the $(2+4)$ TMG ($M\!=\!2$, $N\!=\!4$) at the first magic angle. When the stacking chiralities are the same, the total Chern number of the two flat bands for the $K$ ($K^{\prime}$) valley equals to $-2$ ($+2$); while if the stacking chiralities are opposite, then the total Chern number of the two flat bands equals to $+4$ ($-4$) for the $K$ ($K^{\prime}$) valley. Again, this is in perfect agreement with Eq. (8). We have also numerically tested the other $(M+N)$-layer TMG with $M,N$ extending from 1 to 5, and they are all consistent with Eq. (8).

We have also considered the more realistic situation in which the interlayer hopping is given by Eq. (7) instead of Eq. (4). Since the chiral symmetry is broken in Eq. (7), the Chern-number hierachy given by Eq. (8) is no longer exact. However, since the Chern number of concern is the total Chern number of the two low-energy flat bands, it should remain unchanged as long as the two low-energy flat bands remain separated from the other high-energy bands. We have numerically calculated the total valley Chern numbers of the two low-energy bands in $(M+N)-$layer TMG systems ($M,N$ varies from 1 to 5) at $\theta\!=\!1.05^{\circ}$ using the more realistic interlayer hopping term Eq. (7), and find that the Chern-number hierarchy of Eq. (8) remains correct for the cases of $(M=1,N=2)$, $(M=1,N=3)$, $(M=1,N=4)$, $(M=1,N=5)$, $(M=2,N=2)$, $(M=2,N=3)$, and is partially correct (for one stacking configuration) for $(M\!=\!3,N\!=\!3)$, $(M\!=\!2,N\!=\!4)$, and $(M\!=\!3,N\!=\!5)$.

II.2 Gate tunable quantum valley Hall effect

The valley Chern numbers given by Eq. (8) can be further tuned by applying an vertical electric potential $V_{\perp}$. Taking the case of $(2+2)$-layer TMG and $(2+1)-$layer TMG as an example, we study the dependence of the valley Chern numbers of each of the two low-energy bands on the vertical electric potential $V_{\perp}$ at $\theta\!=\!1.05^{\circ}$. The valley Chern number of each band is denoted by $C_{n\alpha\alpha^{\prime}}^{K}$, with the band index $n=1,2$, the stacking chirality $\alpha,\alpha^{\prime}\!=\!\pm$, and the valley index $K$.

In Table  1 we show the dependence of $C_{1\alpha\alpha^{\prime}}^{K}$, $C_{2\alpha\alpha^{\prime}}^{K}$ on the vertical electric potential $V_{\perp}$ (in units of meV) for $(2+2)$-layer TMG. The valley Chern numbers are calculated using the effective Hamiltonian Eq. (2) with the more realistic interlayer hopping Eq. (7). When the stacking chiralities are the same (both $+$), the Chern number of each of the two low-energy bands becomes $\pm 3$ once a small $V_{\perp}\!\sim\!12\,$meV is applied, whereas the total Chern number of the two bands still sums to zero. As $|V_{\perp}|$ increases, the valley Chern number of the two flat bands are changed to $\pm 2$ at $|V_{\perp}|\!=\!24\,$meV, then become [math] and $\pm 1$ respectively at $|V_{\perp}|\!=\!36\,$meV, and both become $\pm 1$ when $|V_{\perp}|\gtrapprox 48\,$meV. On the other hands, when the stacking chiralities are opposite ($\alpha=+$ and $\alpha^{\prime}=-$), the total valley Chern number of the two bands equals to $-2$ at $|V_{\perp}|\!=\!0$, and remains unchanged for $|V_{\perp}|\lessapprox\!40\,$meV. Then the total valley Chern number of the two bands becomes $+1$ for $|V_{\perp}|\gtrapprox\!48\,$meV. Table 1 indicates that the topological phases of $(2+2)-$layer TMG is highly tunable by gate voltage, which is qualitatively in agreement with the results reported in Refs. Zhang et al., 2019b and Lee et al., 2019. Moreover, here we have also shown that it is sensitive to the stacking configurations.

In Table 2 we show $C_{n+}^{K}$ ($n=1,2$) vs. $V_{\perp}$ for $(2+1)$-layer TMG at $\theta\!=\!1.05^{\circ}$, where the subscript “$+$” means that the bottom bilayer has $+$ stacking chirality, and $n=1,2$ is the band index. When the vertical potential $V_{\perp}\!=\!0$, $C_{1+}^{K}\!=\!+1$ and $C_{2+}^{K}\!=\!0$, which is consistent with Eq. (8). Once a small positive or negative $V_{\perp}\!\sim\!10\,$meV is applied, $C_{1+}^{K}$ is changed to $+2$ and $C_{2+}^{K}$ becomes -1. However, when $|V_{\perp}|\!\gtrapprox\!20\,$meV, the valley Chern numbers become highly dependent on the sign of $V_{\perp}$ due to the asymmetric stacking configuration of $(2+1)$-layer TMG. Thus distinct topological phases could be realized by reversing the gate potentials.

III Verifications using an empirical tight-binding model

The flat bands at the first magic angle (Fig. 2) and the Chern-number hierarchy (Eq. (8)) in TMG have been verified using a realistic microscopic tight-binding (TB) model. To be specific, the hopping parameter between two $p_{z}$ orbitals at different carbon sites $i$ and $j$ in the multilayer system is expressed in the Slater-Koster form

[TABLE]

where $V_{\sigma}=V_{\sigma}^{0}\,e^{-(r-d_{c})/\delta_{0}}$, and $V_{\pi}=V_{\pi}^{0}\,e^{-(r-a_{0})/\delta_{0}}$. $\mathbf{d}=(d_{x},d_{y},d_{z})$ is the displacement vector between the two carbon sites. $a_{0}=a/\sqrt{3}=1.42\,$Å, $d_{c}=3.35\,$Å is the interlayer distance in AB-stacked bilayer graphene, and $\delta_{0}=0.184\,a$. $V_{\sigma}^{0}=0.48\,$eV and $V_{\pi}^{0}=-2.7\,$eV. The atomic corrugations at the interface between the $M$ layers and the $N$ layers in $(M+N)$ TMG are modeled by Eq. (1), and their effects can be taken into account by plugging Eq. (1) into the hopping parameter shown in Eq. (9).

The bandstructures calculated using the Slater-Koster TB model at $\theta\!\approx\!1.08^{\circ}$ for the $(2+2)$ TMG are shown in Fig. 5. To be specific, the bandstructure of (2+2) TMG with opposite and the same stacking chiralities are shown in Fig. 5(a) and (c) respectively. It is evident that there are four low-energy bands (contributed by the two valleys) which are separated from the other high-energy bands, and the bandwidths are on the order of 10-15 meV, which are greater than those from the continuum model (see Fig. 2(a). This is because in the continuum model only the nearest neighbor interlayer hopping is kept (see Eq. (13)), which imposes a chiral symmetry to the Hamiltonian of Eq. (2). As argued in Appendix A, the chiral symmetry of Eq. (2) would pin the zeroth pseudo LLs emerging from the twisted bilayer at the interface to zero energy, leading to almost vanishing bandwidth as shown in Fig. 2. However, such a chiral symmetry is broken in the realistic Slater-Koster TB model, and the bandwidth of the flat bands is expected to be enhanced due to the presence of further neighbor interlayer hoppings.

In Fig. 6(b) and (d) we show the Wilson loops of the four low-energy bands of (2+2) TMG at $\theta\!\approx\!1.08^{\circ}$ calculated using the Slater-Koster TB model. Fig. 6(b) ((d)) denotes the case with opposite (the same) stacking chiralities with the valley Chern number being $\pm 2$ ([math]). The topological equivalence between the band structures obtained by the more accurate tight binding model as shown in Fig. 6 and by the continuum model as shown in Fig. 4(a)-(b) provides a strong supporting fact for the Chern number hierarchy given by Eq. (8).

IV Valley-contrasting orbital magnetizations and valley-polarized states

The valley Chern numbers given by Eq. (8) implies opposite orbital magnetizations for the two monolayer valleys $K$ and $K^{\prime}$. In particular, according to the “modern theory” of orbital magnetizations Thonhauser et al. (2005); Xiao et al. (2005); Ceresoli et al. (2006); Shi et al. (2007), the orbital magnetization of the $(M+N)$-layer TMG for either $K$ or $K^{\prime}$ valley can be expressed as

[TABLE]

where $\epsilon_{n\mathbf{k}}$ and $|u_{n\mathbf{k}}\rangle$ are the eigenenergies and (the periodic part of) Bloch eigenstates of a $(M+N)$-layer TMG Hamiltonian (denoted by $H$) for either $K$ or $K^{\prime}$ valley, $\mu$ is the chemical potential, and $H_{\mathbf{k}}\!=\!e^{-i\mathbf{k}\cdot\mathbf{r}}He^{i\mathbf{k}\cdot\mathbf{r}}$. Since the $K$ and $K^{\prime}$ valleys are transformed to each other by a time-reversal operation, it is naturally expected that the valley-contrasting Chern numbers shown in Eq. (8) would lead to opposite orbital magnetization for the two valleys. The valley-contrasting orbital magnetizations further suggests that when the valley degeneracy is lifted by external magnetic fields or by Coulomb interactions, a valley-polarized state with non-vanishing or even quantized anomalous Hall conductivity will be generated due to the nonzero valley Chern numbers.

We have exploited this idea using the continuum model given by Eq. (2) with the interlayer hopping modeled by Eq. (7). To be explicit, we consider the case of $(2+2)$-layer TMG at the first magic angle $\theta\!=\!1.05^{\circ}$. We have calculated the orbital magnetizations ($M_{z}$) of the two low-energy bands for the $K$ valley using the Hamiltonian Eq. (2) with the interlayer hopping modeled by Eq. (7). The dependence of $M_{z}$ on the chemical potential $\mu$ is plotted in Fig. 7(a). The red diamonds in Fig. 7(a) represent the situation that the bottom bilayer and the top bilayer have opposite stacking chiralities with the valley Chern number $\pm 2$ (see Eq. (8)). In this case the magnitude of $M_{z}$ is large, which is on the order of $10\,\mu_{B}$ per moiré supercell when the two flat bands are completely filled. The large orbital magnetization is a manifestation of the band topology on the moiré pattern. In particular, the non-vanishing Chern number of the $K$ or $K^{\prime}$ valley implies that the ground state at a given filling would possess chiral current loops. The characteristic radius of the current loop is on the order of the moiré length scale $L_{s}$, which is associated with large orbital angular momentum $L_{z}\sim\mathbf{r}\times\mathbf{p}$, with $|\mathbf{p}|\!\sim\!\hbar K$ and $|\mathbf{r}|\!\sim\!L_{s}\!\gg\!a$. Therefore, the orbital magnetization generated by the current loops on the moiré scale is expected to be much greater than that on the microscopic lattice scale. We also note that $M_{z}$ increases almost linearly with $\mu$ when $\mu$ is in the gap above the two flat bands, which is a signature of the non-vanishing Chern number Ceresoli et al. (2006). On the other hand, the blue circles in Fig. 7(a) denote the case that the bottom bilayer and the top bilayer have the same stacking chirality with vanishing valley Chern number. In this case the orbital magnetization vanishes identically for any chemical potential due to the presence of a $\mathcal{P}\mathcal{T}$ symmetry ($\mathcal{P}$ is the 3D spatial inversion operation). It worth to note again that the orbital magnetization plotted in Fig. 7(a) is for the $K$ valley. The orbital magnetization for the $K^{\prime}$ valley is just opposite to that of the $K$ valley.

The large orbital magnetizations imply that the valley degeneracy of the system can be easily lifted by a weak external magnetic field or by spontaneous valley symmetry breaking from Coulomb interactions, leading to a valley-polarized (quantum) anomalous Hall state. A rough estimate reveals that a magnetic field of $2\,$T would give rise to an orbital (or valley) Zeeman splitting $\sim 2\,$meV (15% of the bandwidth), which would lead to a state with considerable valley polarization and anomalous Hall effect. Such a valley polarized anomalous Hall state is expected to possess chiral current loops, which are responsible for the large orbital magnetization shown in Fig. 7(a). Here we assume that the $K$ valley is 100% polarized either due to the presence of an external magnetic field, or due to spontaneous valley symmetry breaking from Coulomb interactions, and we calculate the local current density $\mathbf{J}(\mathbf{r})$ and the current-induced local magnetic field $B(\mathbf{r})$ with the two flat bands of the $K$ valley being completely filled. In Fig. 7(b) we show the distributions of $\mathbf{J}(\mathbf{r})$ and $B(\mathbf{r})$ within the moiré Wigner-Seitz cell for the $(2+2)$-layer TMG with opposite stacking chiralities, and with 100% valley polarization. The small filled circles in Fig. 7(b) represent the discretized real-space positions, with the color coding denoting the strength of the local magnetic field in units of Tesla. The black arrows represent the local current densities whose amplitudes are proportional to the lengths of the arrows. Clearly the valley polarized ground state possesses chiral current loops circulating around the $AA$ region. These circulating current loops further generate magnetic fields in the $AA$ region with the magnitude $\sim\!10^{-5}\,$T, which may be a strong experimental evidence for the non-vanishing valley Chern number in $(2+2)$-layer TMG with opposite stacking chiralties.

In Fig. 7(c) we plot the current densities $J_{x}$ (blue circles) and $J_{y}$ (red diamonds) along the real-space path $M_{r}-\Gamma_{r}-K_{r}$ in units of $\mu\textrm{A}/\mbox{\AA}^{2}$, where the $M_{r}$, $\Gamma_{r}$, and $K_{r}$ points are marked in Fig. 7(b). It is interesting to note that $J_{y}$ almost vanishes identically along the $M_{r}-\Gamma_{r}$ path due to the winding pattern of the current. The maximal magnitude of the current density $\sim 6\times 10^{-4}\mu\textrm{A}/\mbox{\AA}^{2}$. In Fig. 7(d) we show the local magnetic field plotted along the $M_{r}-\Gamma_{r}-K_{r}$ path. Clearly the magnetic field has a peak centered at $\Gamma_{r}$ (the $AA$ point) with the magnitude $\sim\!10^{-5}$ T. The details of the computing the current densities and local magnetic fields are presented in Appendix C.

V Summary

To summarize,we have studied the electronic structures and topological properties of the $(M+N)$-layer TMG system. We have proposed that, with the chiral approximation, there always exists two low-energy bands whose bandwidths become minimally small at the magic angle of TBG. We have further shown that the two flat bands in the TMG system are topologically nontrivial, and exhibits a Chern-number hierarchy. In particular, when the stacking chiralities of the $M$ layers and the $N$ layers are opposite, the total Chern number of the two low-energy bands for each valley equals to $\pm(M+N-2)$ (per spin). On the other hand, if the stacking chiralities of the $M$ layers and the $N$ layers are the same, then the total Chern number of the two low-energy bands for each valley is $\pm(M-N)$ (per spin). The non-vanishing valley Chern numbers are associated with large and valley-contrasting orbital magnetizations along $z$ direction, which implies that the valley degeneracy can be lifted by weak external magnetic fields or by Coulomb interactions, leading to a valley-polarized anomalous Hall state. Such a valley polarized state is associated with chiral current loops circulating around the $AA$ region, which generates local magnetic fields peaked at the $AA$ region. The local magnetic fields induced by the chiral current loops may be a robust experimental signature of the valley polarized state with non-vanishing Chern number. Our work is a crucial step forward in understanding the electronic properties of twisted multilayer graphene. The universal magic angles and the Chern-number hierarchy proposed in this work make the TMG system a perfect platform to study the interplay between electrons’ Coulomb correlations and nontrivial band topology.

Note added: During the preparation of our manuscript, we note two recently posted works Ref. Lee et al., 2019 and Ref. Koshino, 2019 . In the former, the electronic structures, superconductivity, and correlated insulating phase have been discussed for twisted double bilayer graphene with the same stacking chirality ($AB$-$AB$ stacking). In the latter, the bandstructures and valley Chern numbers of twisted double bilayer graphene with different stacking chiralities have been discussed.

Acknowledgements.

J.L. and X.D. acknowledge financial support from the Hong Kong Research Grants Council (Project No. GRF16300918). We thank Hongming Weng for invaluable discussions.

Appendix A The flat bands and magic angles in the TMG system

In this appendix we explain the origin of the flat bands and universal magic angles in the $(M+N)$-layer TMG system. After some gauge transformations, the constant wavevectors $\mathbf{K}_{M},\mathbf{K}_{N}$ in Eq. (3) can be removed. Taking the case of $(2+2)$-layer TMG with the same stacking chiralities as an example, the effective Hamiltonian is explicitly written as (after the gauge transformation)

[TABLE]

where $h_{0}(\hat{\mathbf{k}})=-\hbar v_{F}\hat{\mathbf{k}}\!\cdot\!\mathbf{\sigma}$, and $h_{\alpha}$ and $U$ are defined in Eq. (4) and Eq. (5) of the main text. Then we make the following unitary transformation to the basis functions of Eq. (11) (and more generally, to those of Eq. (2) of the main text)

[TABLE]

, where $|{\psi}_{ls}\rangle$ denotes the Bloch state at the $K$ point from the $l$th layer and the $s$ sublattice. Applying the unitary transformation Eq. (12) to Eq. (11) (letting $M\!=\!2$), then expanding the interlayer coupling term $U(\mathbf{r})$ to the linear order of $r/L_{s}$, Eq. (11) becomes

[TABLE]

where the pseudo vector potential $\mathbf{A}\!=\!(2\pi u_{0}^{\prime})/(L_{s}ev_{F})\,(\,y\,,-x\,)$ Liu et al. (2018a). Note that the diagonal blocks $h_{0}(\hat{\mathbf{k}}\pm\frac{e}{\hbar}\mathbf{A})\!=\!-\hbar v_{F}(\hat{\mathbf{k}}\pm e\mathbf{A}/\hbar)$ are equivalent to the Dirac fermions coupled with opposite pseudo magnetic fields, which would generate pseudo LLs of opposite Chern numbers $\pm 1$ Liu et al. (2018a). In particular, the two zeroth pseudo LLs have opposite sublattice polarizations, thus the intrasublattice coupling term $\pm 3iu_{0}$ in Eq. (13) vanishes in the subspace of the zeroth pseudo LLs Liu et al. (2018a). Therefore, if we re-write Eq. (13) in a mixed basis consisted of the Dirac states from the first and fourth layers, and the zeroth pseudo LLs from the second and third layers, then $\widetilde{H}_{+,+}^{K}(2+2)$ becomes

[TABLE]

Each element of Eq. (14) is a $3\times 3$ matrix. In particular,

[TABLE]

and

[TABLE]

where $\eta(k)$ denotes the coupling between the zeroth pseudo LL and the Dirac states of the first (fourth) layer, which can be expressed as some integral over the eigenfunctions of the zeroth pseudo LLs and the Bloch functions, and $k$ is index indicating the zeroth Landau-level degeneracy. Note that we have dropped the higher pseudo LLs in Eq. (14).

The Diagonalizations of Eq. (16) and Eq. (15) would always yield two decoupled zero modes at any $k$. These two zero modes originate from the two zeroth pseudo LLs contributed by the $M$th and $(M+1)$th twisted bilayer ($M\!=\!2$ for Eq. (14)), and they stay at zero energy even after being coupled with the other layers due to the chiral symmetry of the effective Hamiltonian of TMG: all the matrix elements in Eq. (2) and Eq. (11)-Eq. (14) are intersublattice. As a consequence, if we apply the gauge transformation such that all the basis functions at the $B$ sublattice changes sign, then the total Hamiltonian and eigenenergies would change sign as well. However, the eigenenergies are supposed to be invariant under such a gauge transformation, which thus enforces that both $E(\mathbf{k})$ and $-E(\mathbf{k})$ have to be the eigenenergies of the Hamiltonian. Therefore, a zero mode would stay at zero energy as long as the chiral symmetry is preserved. Similar argument applies to any $(M+N)$-layer TMG systems with either opposite or the same stacking chiralities. As long as the chiral symmetry is preserved, the zeroth pseudo LLs emerging from the interface between the $M$ layers and $N$ layers would be pinned to zero energy.

On the other hands, it is well known that at the magic angles of TBG the bandwidth of the two low-energy bands for each valley is minimal. From the perspective of the pseudo LLs Liu et al. (2018a), it means that at the magic angles, the states within the zeroth pseudo LLs are minimally coupled with each other (and to the higher pseudo LLs), thus they are almost exactly flat. As discussed above, by virtue of the chiral symmetry, these zeroth pseudo LLs that are maximally flat at the magic angles would remain flat even after being coupled with the other layers in the TMG systems. It follows that the magic angles in TBG should be universal in the TMG systems by virtue of the chiral symmetry in Eq. (2).

Appendix B Derivations of the Chern-number hierarchy

In this Appendix we mathematically prove the Chern-number hierarchy given by Eq. (8). As discussed in Sec. II.1, the $(M+N)$-layer TMG system can be decomposed into three subsystems: the TBG at the interface, the $(M-1)$ layers below the interface, and the $(N-1)$ layers above the interface. This is schematically shown in Fig. 5(a). In graphene multilayers with $+$ stacking chirality, the $B$ sublattice of the $n$th layer is strongly coupled with the $A$ sublattice of the $(n+1)th$ layer, forming pairs of bounding and anti-bounding states, leaving the $A$ sublattice of the first layer and the $B$ sublattice of the $N$th layer ($N$ is the number of layers) as two low-energy degrees freedom. Similarly, in graphene multlayers with $-$ stacking chirality, the $B$ sublattice of the 1st layer and the $A$ sublattice of the $N$th layer would be the low-energy degrees of freedom. The low-energy effective Hamiltonians of the $(M-1)$ layer graphene around the $K$ valley with $\alpha$ ($\alpha=\pm 1$) stacking chirality can be obtained by straightforward $(M-1)$th order perturbation theory Min and MacDonald (2008), which is expressed as

[TABLE]

where $t_{\perp}$ is the interlayer hopping within the $(M-1)$ layers. Eq. (17) is in the basis of $|1,A\rangle$ ($|1,B\rangle$) and $|M-1,B\rangle$ ($|M-1,A\rangle$) if the $(M-1)$ layers have $+$ ($-$) stacking chirality. On the other hand, there are two nearly flat bands at the magic angle contributed by the interface TBG. Around the $K_{s}$ and $K_{s}^{\prime}$ points these two flat bands are equivalent to zeroth pseudo LLs, and possess opposite sublattice polarizations as argued in Ref. Liu et al., 2018a. Let us first assume the coupling between the $(M-1)$ layers and the interface TBG is small, then the low-energy effective Hamiltonian around the $K_{s}$ point can be expressed as

[TABLE]

where

[TABLE]

and $t_{\textrm{eff}}$ is the low-energy effectively coupling between the states of the $(M-1)$ layers and the flat bands from the interface TBG. In particular, if the stacking chirality is $+$, $t_{\textrm{eff}}$ represents the coupling between the Bloch states from the $B$ sublattice of the $(M-1)th$ layer and one of the two flat bands with $A$ sublattice polarization contributed by the interface TBG. If the stacking chirality is $-$, then $t_{\textrm{eff}}$ represents the coupling between the Bloch states from the $A$ sublattice of the $(M-1)th$ layer and one of the two flat bands with $B$ sublattice polarization contributed by the interface TBG. Note that although $\mathbf{k}$ in Eq. (17) denotes the $\mathbf{k}$ point in the original primitive-cell BZ, while $\mathbf{k}$ in Eq. (18) represents the $\mathbf{k}$ point in the moiré supercell BZ. This is because the form of the low-energy effective Hamiltonian (Eq. (17)) is unchanged after the BZ folding.

Similarly one can write down the low-energy effective Hamiltonian around the $K_{s}^{\prime}$ point, which consists of the two low-energy states from the $(N-1)$ layers above the TBG interface and one of the two flat bands from the TBG interface. More explicitly, it can be expressed as

[TABLE]

where $\alpha^{\prime}=\pm$ represents the stacking chirality of the $N$ layers, and $m_{N-1}$ is given by Eq. (19).

Both Eq. (18) and Eq. (20) can be solved analytically. The eigenenergies of Eq. (18) are expressed as:

[TABLE]

The eigenenergies of Eq. (20) have exactly the same analytic expression but with $(M-1)$ replaced by $(N-1)$. The eigenstates of $H_{K_{s}}^{\alpha}(\mathbf{k})$ are expressed as

[TABLE]

where $g_{M-1}(|\mathbf{k}|)\!=\!\sqrt{\widetilde{t}_{\textrm{eff}}^{2}+|\mathbf{k}|^{2M-2}}$, $\widetilde{t}_{\textrm{eff}}\!=\!t_{\textrm{eff}}/m_{M-1}$, and $\theta_{\mathbf{k}}\!=\!\arctan(k_{y}/k_{x})$. The eigenstates of $H_{K_{s}^{\prime}}^{\alpha^{\prime}}(\mathbf{k})$ are expressed as

[TABLE]

Given that $(k_{x}\pm\alpha k_{y})^{M-1}$ can be rewritten as $|\mathbf{k}|^{M-1}e^{i\pm\alpha\theta_{\mathbf{k}}(M-1)}$, it is straightforward to calculate the Berry connections of the valence states and the conduction states. For the states around the $K_{s}$ point (Eq. (22), it turns out that

[TABLE]

and $A_{3\theta_{\mathbf{k}}}^{M-1}\!=\!A_{1\theta_{\mathbf{k}}}^{M-1}$. For the states around the $K_{s}^{\prime}$ point, we have

[TABLE]

and $A_{3\theta_{\mathbf{k}}}^{N-1}\!=\!A_{1\theta_{\mathbf{k}}}^{N-1}$. It is interesting to note that the valence and conduction states $|\psi_{1\mathbf{k}}\rangle$ and $|\psi_{3\mathbf{k}}\rangle$ have the same Berry connections, therefore they have the same Berry phase by virtue of the chiral symmetry of Eq. (18) and Eq. (20). Taking the limit $\widetilde{t}_{\textrm{eff}}\!\to\!0$, from Eq. (24-25) it follows that $A_{1\theta_{\mathbf{k}}}^{M-1}\!=\!A_{3\theta_{\mathbf{k}}}^{M-1}\!\to\!-\alpha(M-1)/2$, and $A_{1\theta_{\mathbf{k}}}^{N-1}\!=\!A_{3\theta_{\mathbf{k}}}^{N-1}\!\to\!\alpha^{\prime}(N-1)/2$. Therefore, around $K_{s}$ ($K_{s}^{\prime}$) point, the sum of the Berry fluxes of the conduction and the valence bands equals to $-\alpha(M-1)$ ($\alpha^{\prime}(N-1)$). Then the total Chern number of the conduction and the valence bands equal to $-\alpha(M-1)+\alpha^{\prime}(N-1)$. The Chern number of the two flat bands ($|\psi_{2\mathbf{k}}\rangle$) must cancel the total Chern number of the valence and conduction bands ($|\psi_{1\mathbf{k}}\rangle$ and $|\psi_{3\mathbf{k}}\rangle$), it follows that the total Chern number of the two flat bands equals to $\alpha(M-1)-\alpha^{\prime}(N-1)$ for the $K$ valley. The total Chern number of the two flat bands of the $K^{\prime}$ valley is just opposite to that of the $K$ valley, thus Eq. (8) has been proved.

Appendix C Calculating the local charge current using the continuum model

The current density operator at a position $\mathbf{r}$ is expressed as

[TABLE]

where $\hat{\rho}(\mathbf{r})$ is the local density operator at $\mathbf{r}$, and the velocity operator $\hat{\mathbf{v}}$ satisfies the Schördinger equation $-i\hbar\hat{\mathbf{v}}\!=[H(\mathbf{r}),\mathbf{r}]$. The Hamiltonian at $\mathbf{r}$ is given by Eq. (2), with $\mathbf{k}\!=\!-i\partial_{\mathbf{r}}$. Then it is straightforward to calculate the velocity operator by performing the commutator of $H(\mathbf{r})$ and $\mathbf{r}$. The expectation value of $\hat{\mathbf{J}}(\mathbf{r})$ then equals to

[TABLE]

where $\hat{\rho}\!=\!\sum_{n\mathbf{k}}|\psi_{n\mathbf{k}}\rangle\langle\psi_{n\mathbf{k}}|\theta(\mu-\epsilon_{n\mathbf{k}})$ is the density matrix at zero temperature with the chemical potential $\mu$, with $|\psi_{n\mathbf{k}}\rangle$ and $\epsilon_{n\mathbf{k}}$ being the $n$th eigenstate and eigenenergy of the Hamiltonian Eq. (2) at the $\mathbf{k}$ point in the moiré BZ. $|s\mathbf{G},\mathbf{k}\rangle$ is the plane-wave function, where $\mathbf{G}$ represents a reciprocal lattice vector of the moiré cell, and $s$ is the index for the layer and sublattice degrees of freedom. To be more explicit, in the plane-wave basis, $\hat{\rho}$, $\rho(\mathbf{r})$ and $\hat{v}$ are expressed as

[TABLE]

where the $C_{s\mathbf{G},n}(\mathbf{k})$ is the plane-wave coefficient of the eigenstate $|\psi_{n\mathbf{k}}\rangle$, i.e., $|\psi_{n\mathbf{k}}\rangle=\sum_{s\mathbf{G}}C_{s\mathbf{G},n}(\mathbf{k})|s\mathbf{G},\mathbf{k}\rangle$. $V$ is the total volume of the system, and $H_{\mathbf{k}}=e^{-i\mathbf{k}\cdot\mathbf{r}}He^{i\mathbf{k}\cdot\mathbf{r}}$, where the Hamiltonian $H$ is given by Eq. (2). Given the current density distribution $\mathbf{J}(\mathbf{r})$, the magnetic field $B(\mathbf{r})$ can be solved using the Ampère’s law, $\nabla\times\mathbf{B}(\mathbf{r})=\mu_{0}\mathbf{J}(\mathbf{r})$, where $\mu_{0}$ is the magnetic permeability of the vacuum.