Magic Angle Hierarchy in Twisted Graphene Multilayers

TL;DR

This paper reveals a hierarchical relationship among magic angles in multilayer twisted graphene systems, showing how these angles relate to the bilayer case and enabling the tuning of multiple flat bands.

Contribution

It introduces a mathematical hierarchy linking multilayer magic angles to the bilayer case and demonstrates tunability of flat bands in multilayer graphene.

Findings

Magic angles in multilayer graphene relate to bilayer angles via specific mathematical sequences.

Multiple flat bands can be achieved simultaneously by tuning the twist angles.

For large layers, a continuum of magic angles exists for small twist angles.

Abstract

When two monolayers of graphene are stacked with a small relative twist angle, the resulting band structure exhibits a remarkably flat pair of bands at a sequence of 'magic angles' where correlation effects can induce a host of exotic phases. Here, we study a class of related models of $n$-layered graphene with alternating relative twist angle $\pm \theta$ which exhibit magic angle flat bands coexisting with several Dirac dispersing bands at the Moir\'e K point. Remarkably, we find that the Hamiltonian for the multilayer system can be mapped exactly to a set of decoupled bilayers at {\it different} angles, revealing a remarkable hierarchy mathematically relating all these magic angles to the TBG case. For the trilayer case ($n = 3$), we show that the sequence of magic angle is obtained by multiplying the bilayer magic angles by $\sqrt{2}$, whereas the quadrilayer case ($n = 4$) has two sequences of magic angles obtained by multiplying the bilayer magic angles by the golden ratio $\varphi = (\sqrt{5} + 1)/2 \approx 1.62$ and its inverse. We also show that for larger $n$, we can tune the angle to obtain several narrow (almost flat) bands simultaneously and that for $n \rightarrow \infty$, there is a continuum of magic angles for $\theta \lesssim 2^o$. Furthermore, we show that tuning several perfectly flat bands for a small number of layers is possible if the coupling between different layers is different. The setup proposed here can be readily achieved by repeatedly applying the "tear and stack" method without the need of any extra tuning of the twist angle and has the advantage that the first magic angle is always larger than the bilayer case.

Figures

Equations

$$H=\sum _{\ell =1}^n{c}_{\ell ,\mathbf{r}}^{†}(-i{v}_F{\mathbf{\sigma }}_{{\theta }_{\ell }}\cdot \nabla )c_{\ell ,\mathbf{r}}+\sum _{\ell =1}^{n-1}{c}_{\ell ,\mathbf{r}}^{†}{T}_{\mathbf{r}}^{\ell ,\ell +1}{c}_{\ell +1,\mathbf{r}}+\text{h.c.},$$

$$T^{ij}({\bm{r}})=\left(\begin{array}[]{cc}w^{ij}_{\textrm{AA}}U_{0}^{ij}({\bm{r}})&w^{ij}_{\textrm{AB}}U_{1}^{ij}({\bm{r}})\\ w^{ij}_{\textrm{AB}}U_{1}^{ij*}(-{\bm{r}})&w^{ij}_{\textrm{AA}}U_{0}^{ij}({\bm{r}})\end{array}\right),$$

$$U_m^{ij}(\mathbf{r})=\sum _{n=1}^3{e}^{im(n-1)\varphi }{e}^{-i{\mathbf{q}}_n^{ij}\cdot (\mathbf{r}-{\mathbf{D}}_{ij})}.$$

$${\mathbf{D}}_{ij}=\frac{{\mathbf{d}}_i+{\mathbf{d}}_j}{2}+i\cot ({\theta }_{ji}/2){\sigma }_y\frac{{\mathbf{d}}_i-{\mathbf{d}}_j}{2}.$$

$$\kappa =w_{\text{AA}}^{ij}/w_{\text{AB}}^{ij}.$$

$$\mathcal{H}=\left(\begin{array}[]{cc}\mathcal{M}&\mathcal{D}^{\dagger}\\ \mathcal{D}&\mathcal{M}\end{array}\right)_{\rm AB}\,,$$

$$\displaystyle\mathcal{D}=\left(\begin{array}[]{cc}-2i\bar{\partial}&WU_{1}({\bm{r}})\\ W^{T}U_{1}(-{\bm{r}})&-2i\bar{\partial}\end{array}\right),$$

$$\displaystyle\mathcal{M}=\kappa\left(\begin{array}[]{cc}0&WU_{0}({\bm{r}})\\ W^{T}U_{0}(-{\bm{r}})&0\end{array}\right)\,.$$

$$W=\left(\begin{array}[]{cccc}\alpha_{12}&0&0&\dots\\ \alpha_{23}&\alpha_{34}&0&\dots\\ 0&\alpha_{45}&\alpha_{56}&\dots\\ \dots&\dots&\dots&\dots\end{array}\right)\,.$$

$$V^{†}\mathcal{H}V=\{\begin{array}{ll}{\mathcal{H}}_{{\lambda }_1}^{(2)}\oplus {\mathcal{H}}_{{\lambda }_2}^{(2)}\cdots \oplus {\mathcal{H}}_{{\lambda }_{n_e}}^{(2)},& \text{n even},\\ {\mathcal{H}}_{{\lambda }_1}^{(2)}\oplus {\mathcal{H}}_{{\lambda }_2}^{(2)}\cdots \oplus {\mathcal{H}}_{{\lambda }_{n_e}}^{(2)}\oplus {\mathcal{H}}^{(1)},& \text{n odd},\end{array}$$

$$W=\alpha ({\delta }_{ij}+{\delta }_{i,j+1}).$$

$$\alpha^{(n)}_{k}=\alpha^{(2)}/\big{(}2\cos\frac{\pi k}{n+1}\big{)}\,.$$

$${\psi }_{\ell }^{(k)}(\text{r})=\sqrt{\frac{2n}{n+1}}(-1{)}^{\ell (k+1)}\sin (\frac{\pi k}{n+1}\ell ){\psi }_{{\lambda }_k}^{\text{TBG}}(\text{r}),$$

$${\alpha }_1^2={\alpha }_1^{(2)}{\alpha }_2^{(2)},2{\alpha }_1^2+{\alpha }_2^2=({\alpha }_1^{(2)}{)}^2+({\alpha }_2^{(2)}{)}^2,$$

$$T^{ij}=\int d\mathbf{r}d{\mathbf{r}}^{\prime }{c}_i^{†}(\mathbf{r})t^{ij}(\mathbf{r},{\mathbf{r}}^{\prime })c_j({\mathbf{r}}^{\prime }),$$

$$c_{\ell }(\mathbf{r})=\sum _{\mathbf{R},\alpha }{\varphi }_{\mathbf{R},\alpha }(R_{-{\theta }_{\ell }}(\mathbf{r}-{\mathbf{d}}_{\ell }))f_{\ell ,\mathbf{R},\alpha }=\sum _{\mathbf{R},\alpha }\varphi (R_{-{\theta }_{\ell }}(\mathbf{r}-{\mathbf{d}}_{\ell })-\mathbf{R}-{\tau }_{\alpha })f_{\ell ,\mathbf{R},\alpha }.$$

$$T^{ij}=\sum _{\mathbf{R},{\mathbf{R}}^{\prime },\alpha ,{\alpha }^{\prime }}{t}^{ij}(R_{{\theta }_i}(\mathbf{R}+{\tau }_{\alpha })+{\mathbf{d}}_i,R_{{\theta }_j}({\mathbf{R}}^{\prime }+{\tau }_{{\alpha }^{\prime }})+{\mathbf{d}}_j)f_{i,\mathbf{R},\alpha }^{†}{f}_{j,{\mathbf{R}}^{\prime },{\alpha }^{\prime }}.$$

$$t^{ij}(\mathbf{r},{\mathbf{r}}^{\prime })=\sum _{\mathbf{q}}{e}^{i\mathbf{q}\cdot (\mathbf{r}-{\mathbf{r}}^{\prime })}{t}_{\mathbf{q}}^{ij}.$$

$$f_{i,\mathbf{R},\alpha }=\sum _{\mathbf{p}\in {\mathrm{BZ}}_i}{e}^{i\mathbf{p}\cdot (\mathbf{R}+{\tau }_{\alpha })}{\psi }_{i,\mathbf{p},\alpha }.$$

$$T^{ij}=\sum _{\alpha ,{\alpha }^{\prime },\mathbf{G},{\mathbf{G}}^{\prime },\mathbf{p},{\mathbf{p}}^{\prime }}{\delta }_{R_{{\theta }_i}(\mathbf{G}+\mathbf{p}),R_{{\theta }_j}({\mathbf{G}}^{\prime }+{\mathbf{p}}^{\prime })}{t}^{ij}(\mathbf{G}+\mathbf{p})e^{i[\mathbf{G}\cdot {\tau }_{\alpha }-{\mathbf{G}}^{\prime }\cdot {\tau }_{{\alpha }^{\prime }}+R_{{\theta }_i}(\mathbf{G}+\mathbf{p})\cdot ({\mathbf{d}}_i-{\mathbf{d}}_j)]}{\psi }_{i,\mathbf{p},\alpha }^{†}{\psi }_{j,{\mathbf{p}}^{\prime },{\alpha }^{\prime }}.$$

$${\psi }_{i,\mathbf{p},\alpha }=\int d\mathbf{r}{e}^{-i(\mathbf{p}-K)\cdot \mathbf{r}}{\tilde{\chi }}_{i,K,\alpha }(\mathbf{r}).$$

$${\psi }_{i,\mathbf{p},\alpha }=\int d\mathbf{r}{e}^{-i(\mathbf{p}-K)\cdot R_{-{\theta }_i}(\mathbf{r}-{\mathbf{d}}_i)}{\chi }_{i,K,\alpha }(\mathbf{r}).$$

$$T^{ij}=\int d\mathbf{r}{c}_{i,\alpha }^{†}(\mathbf{r})[T^{ij}(\mathbf{r}){]}_{\alpha ,{\alpha }^{\prime }}{c}_{j,{\alpha }^{\prime }}(\mathbf{r}),$$

$$[T^{ij}(\mathbf{r}){]}_{\alpha ,{\alpha }^{\prime }}=w^{ij}\sum _{\mathbf{G}=0,{\mathbf{G}}_2,{\mathbf{G}}_3}{e}^{-i(R_{{\theta }_i}-R_{{\theta }_j})(\mathbf{G}+K)\cdot \mathbf{r}}{e}^{i(\mathbf{G}\cdot ({\tau }_{\alpha }-{\tau }_{{\alpha }^{\prime }})+R_{{\theta }_i}(\mathbf{G}+K)\cdot {\mathbf{d}}_i-R_{{\theta }_j}(\mathbf{G}+K)\cdot {\mathbf{d}}_j)}.$$

$$T^{ij}\to T^{ij}{e}^{-i(R_{{\theta }_i}K\cdot {\mathbf{d}}_i-R_{{\theta }_j}K\cdot {\mathbf{d}}_j)}=w^{ij}\sum _{\mathbf{G}=0,{\mathbf{G}}_2,{\mathbf{G}}_3}{e}^{-i(R_{{\theta }_i}-R_{{\theta }_j})(\mathbf{G}+K)\cdot \mathbf{r}}{e}^{i(\mathbf{G}\cdot ({\tau }_{\alpha }-{\tau }_{{\alpha }^{\prime }})+R_{{\theta }_i}\mathbf{G}\cdot {\mathbf{d}}_i-R_{{\theta }_j}\mathbf{G}\cdot {\mathbf{d}}_j)}$$

$$\displaystyle T^{ij}({\bm{r}})=\left(\begin{array}[]{cc}w^{ij}_{\rm AA}U_{0}^{ij}({\bm{r}})&w^{ij}_{\rm AB}U_{1}^{ij}({\bm{r}})\\ w^{ij}_{\rm AB}U_{1}^{ij*}(-{\bm{r}})&w^{ij}_{\rm AA}U_{0}^{ij}({\bm{r}})\end{array}\right)\,,$$

$$U_m^{ij}(\mathbf{r})=\sum _{n=1}^3{e}^{im(n-1)\varphi }{e}^{-i{\mathbf{q}}_n^{ij}\cdot (\mathbf{r}-{\mathbf{D}}_{ij})},m=0,1,{\mathbf{q}}_1^{ij}=2k_D\sin ({\theta }_{ji}/2)R_{{\varphi }_{ij}}(0,-1),{\mathbf{q}}_{2,3}^{ij}=R_{\pm \varphi }{\mathbf{q}}^{ij}$$

$$\varphi =2\pi /3,{\theta }_{ji}={\theta }_j-{\theta }_i,{\varphi }_{ij}=({\theta }_i+{\theta }_j)/2,{\mathbf{D}}_{ij}=\frac{{\mathbf{d}}_i+{\mathbf{d}}_j}{2}+i\cot ({\theta }_{ji}/2){\sigma }_y\frac{{\mathbf{d}}_i-{\mathbf{d}}_j}{2}$$

$$H=\sum _{\ell =1}^n{c}_{\ell ,\mathbf{r}}^{†}\cdot (-i{v}_F{\mathbf{\sigma }}_{{\theta }_{\ell }}\nabla )c_{\ell ,\mathbf{r}}+\sum _{\ell =1}^{n-1}{c}_{\ell ,\mathbf{r}}^{†}{T}^{\ell ,\ell +1}(\mathbf{r})c_{\ell +1,\mathbf{r}}+\text{h.c.},$$

$$\displaystyle\mathcal{H}=\left(\begin{array}[]{ccccc}-iv_{F}\bm{\sigma}_{\theta_{1}}\bm{\nabla}&T^{12}({\bm{r}})&0&\dots&0\\ T^{12{\dagger}}({\bm{r}})&-iv_{F}\bm{\sigma}_{\theta_{2}}\bm{\nabla}&T^{23}({\bm{r}})&\dots&0\\ 0&T^{23{\dagger}}({\bm{r}})&-iv_{F}\bm{\sigma}_{\theta_{3}}\bm{\nabla}&\dots&0\\ \dots&\dots&\dots&\ddots&T^{n-1,n}({\bm{r}})\\ 0&0&0&T^{n-1,n{\dagger}}({\bm{r}})&-iv_{F}\bm{\sigma}_{\theta_{n}}\bm{\nabla}\\ \end{array}\right)\,.$$

Full text

I Introduction

Recently, it was shown that two graphene layers twisted to a special (”magic”) angle exhibit a very interesting range of correlated phenomena including Mott insulating and superconducting phases Cao et al. (2018a, b); Yankowitz et al. (2019). This remarkable discovery has stimulated further extensive research into magic-angle superconductivity and correlated electron states in van der Waals heterostructures Po et al. (2018); Thomson et al. (2018); Zou et al. (2018); Guinea and Walet (2018); Carr et al. (2018); Su and Lin (2018); González and Stauber (2018); Wu et al. (2018a); Efimkin and MacDonald (2018); Yuan and Fu (2018); Xu and Balents (2018); Ochi et al. (2018); Wu et al. (2018b); Zhang et al. (2019); Wu et al. (2018c); Kang and Vafek (2018); Pizarro et al. (2018); Koshino et al. (2018); Kennes et al. (2018); Isobe et al. (2018); Rademaker and Mellado (2018); Qiao and He (2018); Chung et al. (2018); Fidrysiak et al. (2018); Peltonen et al. (2018); Tarnopolsky et al. (2019); Song et al. (2018); Hejazi et al. (2018); Po et al. (2018); Liu et al. (2019); Shen et al. (2019); Cao et al. (2019); Lee et al. (2019) and inspired a vast theoretical and experimental search to extend the family of systems which exhibit similar behavior Amorim and Castro (2018); Zhang et al. (2019); Zuo et al. (2018); Chen et al. (2019); Liu et al. (2019); Shen et al. (2019); Cao et al. (2019); Lee et al. (2019). Finding such systems achieves several goals. First, they expand the family of Moiré systems where correlated physics can be studied in a setting which have several advantages over traditional strongly-correlated systems (easier to fabricate, richer possibilities for tuning the band structures, etc). Second, finding systems which share similarities with TBG, but differ in some details – such as symmetries, bands topology and interaction strength, – can help provide a deeper understanding of the correlated physics in TBG itself. Furthermore, some of these systems may have practical advantages over TBG in terms of the ease of fabrication or tunability of the physical properties.

A distinguishing feature of the TBG physics is the appearance of remarkably flat bands at charge neutrality for magic twists $\theta_{*}$ Tarnopolsky et al. (2019); Bistritzer and MacDonald (2011); Cao et al. (2018a). The existence of such flat bands was predicted in Refs.  Dos Santos et al. (2007); Bistritzer and MacDonald (2011) using an effective continuum model for two graphene layers with twist-independent interlayer coupling. (see also Dos Santos et al. (2007); Bistritzer and MacDonald (2011); Lopes dos Santos et al. (2012); Shallcross et al. (2010) and Li et al. (2010); Trambly de Laissardiere et al. (2010); Lopes dos Santos et al. (2012)). In this model Bistritzer and MacDonald (2011), intra- and inter-sublattice hopping parameters were taken to be equal and band flattening happens at a certain sequence of magic angles for which the renormalized Fermi velocity vanishes at Dirac points. It was, however, recently realized that the effect of lattice relaxation in TBG leads to the expansion of the AB stacking regions relative to the AA regions in the Moiré pattern Carr et al. (2019). As a result, the intra-sublattice hopping parameter $w_{\text{AA}}$ is suppressed relative to the inter-sublattice hopping parameter $w_{\text{AB}}$ at small twists. Crucially, this results in band gap opening and further band flattening, down to the point when the bands can in principle become absolutely flat Tarnopolsky et al. (2019). Although it is understood that the existence of the flat bands is important for the correlated physics, it is still unclear which feature, - band-flattening, band isolation or band topology - is most decisive.

In this work, we report an infinite class of multilayer graphene systems which all manifest the remarkably flat bands and the corresponding magic angles patterns. Such systems, if realized experimentally, would provide a rich playground for correlated physics beyond TBG. It is worth noting that other multilayer systems studied in the literature such as ABC trilayer graphene stacked on hexagonal boron nitride or twisted double bilayer graphene do not exhibit magic angles or flat bands when realistic effects, e.g. trigonal warping terms, are included despite the recent experimental observation of correlated insulting states and superconductivity Liu et al. (2019); Cao et al. (2019); Shen et al. (2019); Lee et al. (2019).

We consider a model of alternating-twist multilayer graphene (ATMG) for which the relative twists between two neighboring layers have the same magnitude but alter in sign (see Fig.1). In general, for a system with $n$ graphene sheets, there will be several Moiré patterns – each one originating from a pair of adjacent layers. The overall periodicity of such patterns is determined by the relative inter-layer twists, while the origin of the pattern is controlled by the relative displacement. In the TBG case Bistritzer and MacDonald (2011); Tarnopolsky et al. (2019) the resulting Moiré physics is shift-independent, but in the case of a more complex multi-Moiré interference in twisted multilayers it can be relevant. In this work, we first focus on the case where the Moiré patterns are aligned so that the flat-band physics is most pronounced. At the end, we discuss how sensitive our results are to layer misalignment and show that the obtained flat bands are relatively stable to the inclusion of layer displacement.

Our main result is that the Hamiltonian for the ATMG with $n=2n_{e}$ ($n=2n_{e}+1$) layers can be mapped exactly to a sum of $n_{e}$ twisted bilayer models (plus a single layer model) at different twist angles. Using this mapping, we find that there are $n_{e}$ sequences of magic angles given by multiplying the TBG magic angles by $2\cos\frac{\pi k}{n+1}$ with $k=1,\dots,n_{e}$, which implies that the first magic angle for $n$-layered ATMG is larger than the first magic angle in TBG by a factor of $2\cos\frac{\pi}{n+1}$ which yields $\sqrt{2}$ for $n=3$ and approaches 2 as $n\rightarrow\infty$. This represents a practical advantage since samples with larger twist angles are generally more stable and easier to fabricate. In addition, the setup proposed here does not require independent tuning of the different relative twist angles since it can be achieved by repeatedly appying the ”tear and stack” method.

II Model

We begin by considering a general system of $n$ graphene layers with the $\ell$-th layer twisted counter-clockwise around a lattice cite by an angle $\theta_{\ell}$ and then displaced by a distance ${\bm{d}}_{\ell}$ relative to a fixed reference. Similar to the bilayer problem, the coupling between layers $i$ and $j$ is characterized by two parameters $w_{\rm AA}^{ij}$ and $w_{\rm AB}^{ij}$ which indicate intra- and intersublattice coupling, respectively. We take these parameters to be generally different between different layers and assume coupling takes place only between nearest neighboring layers. The resulting low-energy effective Hamiltonian reads (see Appendix A for details)

[TABLE]

where $v_{F}\approx 10^{6}$ m/s is the monolayer graphene Fermi velocity, ${\bm{\sigma}}_{\theta_{\ell}}=e^{\frac{i}{2}\theta_{\ell}\sigma_{z}}{\bm{\sigma}}e^{-\frac{i}{2}\theta_{\ell}\sigma_{z}}$ and the interlayer coupling matrix $T^{\ell,\ell+1}({\bm{r}})$ takes the form

[TABLE]

with the Moiré potentials $U_{m=0,1}^{ij}({\bm{r}})$ defined as

[TABLE]

Here, ${\bm{q}}^{ij}_{1}=2k_{D}\sin(\theta_{ji}/2)R_{\phi_{ij}}(0,-1)$, ${\bm{q}}^{ij}_{2,3}=R_{\pm\phi}{\bm{q}}^{ij}_{1}$, $\phi=2\pi/3$, $R_{\theta}=e^{-i\theta\sigma_{y}}$ denotes the counter-clockwise rotation operator with angle $\theta$, and $k_{D}={4\pi}/{3\sqrt{3}a}$ is the Dirac momentum of the monolayer graphene with lattice constant $a=1.42$ Å. We also introduced the auxiliary angle variables $\theta_{ji}=\theta_{j}-\theta_{i}$ and $\phi_{ij}=(\theta_{i}+\theta_{j})/2$. We can write the displacement vector ${\bm{D}}_{ij}$ of the Moiré pattern as

[TABLE]

In the bilayer case ($n=2$), the Hamiltonian (1) reduces to the TBG Hamiltonian Dos Santos et al. (2007); Bistritzer and MacDonald (2011); Tarnopolsky et al. (2019) up to the gauge transformation $c_{{\ell}}\rightarrow c_{{\ell}}e^{iR_{\theta_{{\ell}}}K\cdot{\bm{d}}_{{\ell}}}$. The advantage of the form we consider here is that it makes it clear how the layer displacements ${\bm{d}}_{i}$ enter the Hamiltonian by shifting the corresponding Moiré potentials. For $n$ layers, there are $n-1$ shift variables ${\bm{D}}_{\ell,\ell+1}$, $\ell=1,\dots,n-1$, one of which can be removed by redefining the origin, leaving $n-2$ variables which influence the spectrum. This is the reason why the shift vectors were unimportant in the bilayer case in contrast to the multilayer case considered here.

In general, the potential (2) will generate several overlapping Moiré patterns generated by the different angles and shift vectors between consecutive layers. For most of this paper, we will focus on the case of ”unshifted” ATMG which corresponds to the choice $\theta_{\ell}=(-1)^{\ell}\theta/2$ and ${\bm{d}}_{\ell}={\bm{d}}$ such that the nearest neighboring layers are aligned and have alternating relative twists of $\pm\theta$. In this case, $\phi_{\ell,\ell+1}=0$, $\theta_{\ell+1,\ell}=(-1)^{\ell+1}\theta$, ${\bm{D}}_{\ell,\ell+1}={\bm{d}}$ and there is a single Moiré pattern similar to the bilayer problem. We also assume that the ratio between $w_{\textrm{AA}}^{ij}$ and $w_{\textrm{AB}}^{ij}$ couplings is layer-independent and denote it as

[TABLE]

Assuming small twist angle $\theta$, we can neglect the phase factor in the Pauli matrices ${\bm{\sigma}}_{\theta_{\ell}}\to{\bm{\sigma}}$ and get rid of the angular dependence by introducing the dimensionless variables $\alpha_{ij}=w_{\textrm{AB}}^{ij}/(v_{F}k_{D}\theta)$ leading to111This assumption is unnecessary in the chiral limit ($\kappa=0$) since the phase can be removed by a gauge transformation.

[TABLE]

where AB indicates the matrix is in the sublattice space. The operators $\mathcal{D}$ and $\mathcal{M}$ are given by

[TABLE]

Here, we have rescaled the Hamiltonian so that all energies are measured in units of $v_{F}k_{D}\theta$. We also rescaled the coordinates so that they are measured in terms of the Moiré length scale ${\bm{r}}\rightarrow k_{D}\theta{\bm{r}}$ and introduced the derivatives $\partial$ and $\bar{\partial}$ relative to the dimensionless complex variable $z=x+iy$. The potentials $U_{m}({\bm{r}})$ are given by (3) with ${\bm{D}}_{ij}=0$ and ${\bm{q}}^{ij}_{n}=R_{(n-1)\phi}(0,-1)$. The operators $\mathcal{D}$ and $\mathcal{M}$ act on vectors which has the form $\psi=(\psi_{o},\psi_{e})^{T}$ where $\psi_{o/e}$ contain the wave functions for the odd/even layers given explicitly as $\psi_{o}=(\psi_{1},\psi_{3},\dots,\psi_{2n_{o}-1})^{T}$ and $\psi_{e}=(\psi_{2},\psi_{4},\dots,\psi_{2n_{e}})^{T}$ where $n_{e}=\lfloor n/2\rfloor$ and $n_{o}=\lceil n/2\rceil$ are numbers of even and odd layers. The matrix $W$ is $n_{o}\times n_{e}$ layer hopping matrix and given by

[TABLE]

III Results

III.1 Reduction to the bilayer problem

We now show that the Hamiltonian of the multilayer problem with $n=2{n_{e}}$ ($n={2n_{e}+1}$) layers can be mapped exactly to a direct sum of ${n_{e}}$ bilayer Hamiltonians (plus a single layer Hamiltonian). This is done by writing the singular value decomposition of $W$ as $W=A\Lambda B^{\dagger}$ where $A$ and $B$ are $n_{o}\times n_{o}$ and $n_{e}\times n_{e}$ unitary matrices, respectively, and $\Lambda$ is an $n_{o}\times n_{e}$ matrix with $\lambda_{k}$, $k=1,\dots,n_{e}$ on the diagonal and zeros everywhere else ($\lambda_{k}$ are square roots of the eigenvalues of $W^{T}W$). Applying the unitary transformation given by $V={\rm diag}(A,B)$ in the odd/even space to the Hamiltonian (6) yields

[TABLE]

where $\mathcal{H}^{(2)}_{\alpha}$ is the bilayer Hamiltonian with coupling parameters $\alpha=w_{\textrm{AB}}/(v_{F}k_{D}\theta)$ and $\kappa$ and $\mathcal{H}^{(1)}$ is the Hamiltonian for a single graphene layer.

A consequence of the preceding discussion is that the spectrum of the multilayer problem with coupling matrix $W$ is given by the union of the spectra of several bilayer problems whose coupling parameters are given by the eigenvalues of the matrix $\sqrt{W^{T}W}$ (in addition to a single layer graphene dispersion if the number of layers is odd). Moreover, the eigenstates of the multilayer problem are easily obtainable from the eigenstates of the single layer problem. This applies particularly for the case of flat bands where the eigenstates were shown to have a simple form Tarnopolsky et al. (2019).

III.2 The chiral limit

The chiral model for twisted bilayer graphene where the same-sublattice coupling set to zero $w_{\text{AA}}=0$ (or equivalently $\kappa=0$ in this work) was introduced in Ref. Tarnopolsky et al. (2019). It was shown that this model captures the essential phenomenology of magic angles where the different notions of flatness (vanishing Fermi velocity, minimum bandwidth, maximum band gap) all coincide due to the appearance of perfectly flat bands for special (magic) values of the dimensionless coupling $\alpha=w_{\text{AB}}/(v_{F}k_{D}\theta)$. This model is one of the simplest models exhibiting magic angle flat bands and its applicability to TBG is supported by the observation that lattice relaxation tends to reduce the size of AA regions relative to AB regions Carr et al. (2019), thus suppressing the value of $w_{\text{AA}}$ (intrasublattice coupling) relative to $w_{\text{AB}}$ (intersublattice coupling).

Let us first consider the standard setting where all interlayer couplings are the same $\alpha^{ij}=\alpha$. In this case, the layer hopping $n_{o}\times n_{e}$ matrix $W$ is given by

[TABLE]

The eigenvalues of $\sqrt{W^{T}W}$ can be easily computed for any number of layers $n$ and they are given by $\lambda_{k}=2\cos(\frac{\pi k}{n+1})\alpha$, $k=1,\dots,n_{e}$. Thus, the ATMG with $n$ layers has $n_{e}$ sequences of magic angles given by

[TABLE]

Here, $\alpha^{(2)}$ is the sequence of magic angles in the bilayer problem which was computed in Tarnopolsky et al. (2019) as $\alpha^{(2)}=0.586,2.221,3.75,5.276,6.795,\dots$. The magic angle sequence for any $n$ can then be easily computed as shown in Fig. 2 for $n$ up to 6.

It should be noted that the mapping to the bilayer problem can also be used to explicitly write the wave functions of the multilayer system in terms of their bilayer counterparts. For the eigenvalue $\lambda_{k}$, the corresponding eigenfunctions for all layers $\ell=1,\dots,n$ are given by

[TABLE]

normalized as $\sum_{\ell=1}^{n}|\psi_{\ell}^{(k)}(\textbf{r})|^{2}=n|\psi^{\textrm{TBG}}_{{\lambda_{k}}}(\textbf{r})|^{2}$ and $\psi^{\textrm{TBG}}_{\lambda_{k}}(\textbf{r})$ is the eigenfunction of the TBG Hamiltonian $\mathcal{H}^{(2)}_{\lambda_{k}}$.

For the trilayer case ($n=3$), the only eigenvalue of $\sqrt{W^{T}W}$ is $\lambda_{1}=\sqrt{2}\alpha$. Eq. 14 implies that the system is equivalent to the sum of a bilayer problem with coupling $\sqrt{2}\alpha$ and a single layer problem. As a result, we can immediately read off the magic angles, where a perfectly flat band appears, to be $\alpha^{(3)}=\alpha^{(2)}/\sqrt{2}=0.414,1.57,2.65,3.731,4.805,\dots$. This is verified in Fig. 3, where the band structure is computed numerically for the first two magic angles for the trilayer problem showing the existence of a perfectly flat band.

One particularly interesting feature here is that for the trilayer graphene, the first magic angle is larger by a factor of $\sqrt{2}$ compared to the TBG, which represents an experimental advantage. Apart from the scaling of the magic angles, the trilayer system differs from TBG in two main aspects. First, the flat band here coexists with a dispersing Dirac cone. If realized experimentally, this feature will distinguish the physics of the trilayer system from the TBG physics, and could help to elucidate whether band flatness or band isolation plays the bigger role in the correlated physics. Second, while the band structure for the flat band looks identical to the TBG band structure at the first magic angle $\theta\approx 1.08^{o}$, the actual scale for the Moiré pattern is determined by the actual angle $\theta\approx 1.53^{o}$ which determines the scale of the gaps and the interaction. We stress that the mapping does not only apply for the spectra but also for the wave functions. As a result, the physics of the trilayer model (including the interaction effects) will be identical to the physics of TBG with all distances scaled down by a factor of $\sqrt{2}$ and with an extra Dirac band from an individual graphene layer.

For the quadrilayer case ($n=4$), the matrix $\sqrt{W^{T}W}$ has two eigenvalues $\lambda_{1,2}=\alpha\varphi^{\pm 1}$, where $\varphi=(\sqrt{5}+1)/2$ is the golden ratio, yielding two sequences of magic angles $\alpha^{(4)}=\alpha^{(2)}/\varphi=0.362,1.373,2.318,3.261,4.2,5.138,5.075,\dots$ and $\alpha^{\prime(4)}=\alpha^{(2)}\varphi=0.948,3.594,6.069,8.537,\dots$. The quadrilayer ATMG maps to a sum of two TBGs and these two sequences correspond to points at which one of these two twisted bilayers hits a magic angle. The largest magic angle (the smallest $\alpha$) in this case is $\theta\approx 1.75^{\circ}$ – which is larger than the bilayer and trilayer cases.

We note that since the different magic angle sequences for a given $n$ have incommensurate periods, we can find some values of $\alpha$ which is close to several magic angles from different sequences simultaneously. This happens for example for $n=5$ for $\alpha\approx 2.2$ which is very close to the third magic angle in the first sequence (2.165) and the second magic angle in the second sequence (2.221). Another example happens when $n=6$ and $\alpha\approx 1.275$ which is very close to the second magic angle in the first sequence (1.23) and the first magic angle in the third sequence (1.32). In both cases, there are two pairs of very narrow bands coexisting at [math] as shown in Fig. 4.

It is instructive to consider the limit of large number of layers $n\rightarrow\infty$. In this case, the eigenvalues $\lambda_{k}=2\cos(\frac{\pi k}{n+1})\alpha$ of the matrix $\sqrt{W^{T}W}$ form a continuum from [math] to $2$, which implies that there is a continuum of magic angles: whenever $\alpha>\alpha^{(2)}_{1}/2\approx 0.293$, we are always arbitrarily close to a magic angle descending from the first magic angle of TBG where at least a single band is perfectly flat band. Similarly, there is a flat band deriving from the second magic angle for $\alpha>\alpha^{(2)}_{2}/2\approx 1.11$. In general, there will be exactly $k$ perfectly flat bands deriving from the first $m$ TBG magic angles whenever $\alpha^{(2)}_{m-1}/2\leq\alpha<\alpha^{(2)}_{m}/2$, where $\alpha^{(2)}_{m}$ is the $m$-th magic angle of TBG. This suggests an intriguing connection to possible flat-band-related phenomena in some samples of turbostratic graphites, if its layers are naturally assembled in small but very random alternating twists Shallcross et al. (2010).

When the number of layers $n$ is relatively small, it is still possible to achieve several flat bands at [math] simultaneously if we allow for different hopping parameters between different layers. So far, we have only considered the case where all the couplings $\alpha_{ij}=\alpha$ are equal which is naturally expected since all the graphene layers are identical. We now consider instead the possibility that the coupling between layers is non-uniform. For instance, the coupling to the outer layers (top and bottom) may differ slightly from that between inner layers. Another possibility is to artificially tune the couplings by including thin layers of a dielectric material between some of the layers or by depositing adatoms on the top or bottom surfaces to change the interlayer potential. Our purpose in this discussion is to show that this is an interesting theoretical possibility leaving the question of experimental realizability to future studies.

Let us now consider the simplest case with four layers $n=4$ such that the coupling to the outer layers $\alpha_{12}=\alpha_{34}=\alpha_{1}$ is different from the coupling between the middle layer $\alpha_{23}=\alpha_{2}$. In this case, we can achieve two perfectly flat bands simultaneously as follows: we require the two eigenvalues of the matrix $\sqrt{W^{T}W}$ to be equal to the first two magic angles. A simple way to achieve this is to require the determinant and the trace of this matrix to be equal to the product and sum of the first two bilayer magic angles $\alpha^{(2)}_{1,2}$ leading to the equations

[TABLE]

which can be easily solved for $\alpha_{1,2}$ yielding $\alpha_{1}=\sqrt{\alpha^{(2)}_{1}\alpha^{(2)}_{2}}=1.14$ and $\alpha_{2}=|\alpha^{(2)}_{1}-\alpha^{(2)}_{2}|=1.6{4}$. The band structure for this choice of parameters is shown in Fig. 5 showing two perfectly flat bands.

III.3 Switching on intrasublattice couplings $w_{\text{AA}}$

We now consider non-zero intra-sublattice coupling $w_{\text{AA}}$ which in our notation corresponds to non-zero $\kappa$ (cf. Eq. 5). The value of $\kappa$ generally depends on the angle; it is expected to be close to 1 for $\alpha\lesssim 0.25$ then starts decreasing as $\alpha$ is increased Carr et al. (2019). In the following, we will use the values of the relaxation parameter $\kappa$ computed for TBG Carr et al. (2019) as an estimate for the multilayer problem although a more involved ab-initio calculation is needed to refine this value and account for complex lattice relaxation effects in multilayer systems.

In TBG, it is known that the flat bands at the first magic angle are a lot more stable than higher magic angles when $\kappa$ is non-zero Tarnopolsky et al. (2019). Our mapping implies that the flat bands in the multilayer model inherit the stability of the corresponding bilayer flat bands. That is, multilayer magic angles which descend from the first bilayer magic angle are significantly more stable than those descendent from higher magic angles.

For the trilayer case ($n=3$), this means that the first magic angle $\alpha=0.414$ is stable for relatively large values of $\kappa$ including the realistic value $\kappa\approx 0.7-0.8$ estimated in Ref. Carr et al. (2019). The second magic angle $\alpha=1.57$ is however significantly less stable and the flat band gets destroyed for realistic values of $\kappa\approx 0.4$. For quadrilayer ($n=4$), the situation is different since the first two magic angles $\alpha=0.362,0.948$ are descendent from the first magic angle of TBG. In fact, the flat band $\alpha=0.948$ inherits the stability of the bilayer flat band while also having a significantly larger value of relaxation $\kappa\approx 0.5$ compared to the second magic angle in TBG which leads to extra stability. The flat bands for the first two magic angles for $n=3,4$ are shown in Fig. 3 for realistic values of the relaxation parameter $\kappa$ and we can see that they are stable in all cases except the second magic angle for $n=3$ as expected. The effect of $\kappa$ on the $n=4$ setup with different couplings leading to two perfectly flat bands can also be investigated and we find that it is relatively stable for realistic values of $\kappa$ around $0.5-0.6$ (cf. Fig. 5).

IV Experimental realization

So far we have focused on the setting where the layers are perfectly aligned, i.e. the displacement vectors ${\bm{D}}_{l,l+1}$ defined in (4) were assumed to be equal. For experimental realizations, it however very difficult to achieve perfect alignment on the atomic scale. Hence, it is crucial to investigate how our results are affected when we lift this assumption and consider unequal displacement vectors. In this case, the exact mapping to TBG no longer exists. However, the Hamiltonian (1) is still translationally invariant on the Moiré lattice (since the different Moiré potentials are only shifted relative to one another) and we can find the band structure within the Moiré Brillouin zone numerically.

For the trilayer case ($n=3$), we can set the displacement ${\bm{D}}_{12}$ to zero by shifting the origin leaving one relevant displacement ${\bm{D}}_{23}={\bm{D}}$. The bandwidth of the narrow band for $\kappa=0.8$ and for different values of the shift vector ${\bm{D}}$ is shown in Fig. 6. We can see that there is a range of ${\bm{D}}$ around [math] for which the bandwidth remains relatively small. For a reference, we can compare the bandwidth with the Coulomb energy scale given by ${E_{\textrm{Coulomb}}}=\frac{e^{2}\theta}{4\pi\epsilon\epsilon_{0}a}$. In our dimensionless unit, $\epsilon_{\textrm{Coulomb}}=\frac{E_{\textrm{Coulomb}}}{\hbar v_{F}k_{D}\theta}\approx\frac{1}{\epsilon_{0}}$. For $\epsilon_{0}\approx 4$, this scale exceeds the bandwidth for the whole range of ${\bm{D}}$ as shown in Fig. 6 implying strong interaction effects regardless of the layer displacement. We notice that, even when the band is not perfectly flat, it is associated with a very large peak in the density of states as shown in right panel of Fig. 6.

For the quadrilayer case ($n=4$), we can set ${\bm{D}}_{23}$ to zero leaving two different shift vectors ${\bm{D}}_{12}$ and ${\bm{D}}_{34}$ which affect the spectrum. However, if the quadrilayer system is made using the tear and stack method from the same TBG, then the displacements satisfy ${\bm{d}}_{4}-{\bm{d}}_{3}={\bm{d}}_{2}-{\bm{d}}_{1}$. For small enough angles, the displacement vectors are ${\bm{D}}_{12}=i\sigma_{y}({\bm{d}}_{1}-{\bm{d}}_{2})/\theta$ and ${\bm{D}}_{34}=i\sigma_{y}({\bm{d}}_{3}-{\bm{d}}_{4})/\theta={\bm{D}}_{12}$ 222Notice that the layer displacement ${\bm{d}}_{\ell}$ is defined up to translations of the original lattice whereas the effective shifts ${\bm{D}}_{\ell,\ell+1}$ are defined up to the Moiré lattice translations., hence there is only one shift vector ${\bm{D}}={\bm{D}}_{12}={\bm{D}}_{34}$ and we can again investigate how it influences the bandwidth as shown in Fig. 6. Similar to the trilayer case, we see that the bandwidth of the lowest band is small compared to the interaction scale for the whole range of displacements. We also see in this case that the narrow band is associated with a very large peak in the density of states.

It follows from the previous discussion that perfect alignment of the layers is not a requirement for the appearance of flat bands in ATMG since the bandwidth remains reasonably small for the whole range of displacements. It is worth noting that by using the ”tear and stack” trick repeatedly, we can ensure that the twist angles between consecutive layers are exactly equal and opposite without the need of any extra tuning for the angles. For trilayers, the procedure would start by tearing a part of a monolayer sample, twisting and stacking it, then tearing another piece of the base monolayer sample and stacking it on the top of the twisted bilayer without any extra twisting. This ensures that the top and bottom layers are aligned which implies that the two twist angles are opposite to a very good accuracy. The same can be done for quadrilayer by using the tear and stack method starting with a TBG sample, again without any additional twisting. Thus, we expect it to be possible to realize alternating twist angles of equal magnitude to a reasonable accuracy. Adding the fact that the magic angles in the multilayer setting are larger and thus more stable makes our current setting a very promising setup to observe magic angle physics beyond TBG.

V Conclusion

In conclusion, we have introduced a model of twisted multilayer graphene with alternating twist angle focusing on the limit of aligned layers. We have shown that this model for $n=2{n_{e}}$ ($n=2{n_{e}}+1$) layers maps exactly to a sum of ${n_{e}}$ twisted bilayer models (plus a single layer model) with different twist angles. Such mapping enabled us to determine the pattern of magic angles for arbitrary $n$ which is given by multiplying the bilayer magic angles by $2\cos\frac{\pi k}{n+1}$, $k=1,\dots,n_{e}$. Focusing on the trilayer and quadrilayer cases, we found that these models exhibit flat bands coexisting with other dispersing bands at zero energy and showed that such flat bands are relatively stable even when layer misalignment is taken into account. In addition, we found that for relatively large number of layers or when interlayer couplings are different, we can achieve several flat bands at zero energy simultaneously. Moreover, we show that there is a continuum of magic angles for $\theta\lesssim 2^{{\circ}}$ in the limit of very large number of layers $n\rightarrow\infty$. This might suggest an intriguing link to possible flat-band-related phenomena in some samples of turbostratic graphites, if its layers are naturally assembled in small but very random alternating twists. At the end, we discuss possible experimental realizations of the model and show that it can be achieved within current technology by applying the tear and stack method repeatedly. Compared to TBG, the multilayer setting has the advantage that the first magic angle is larger making it easier to realize experimentally.

Note: During the preparation of this manuscript, two related preprints appeared: Mora et al. (2019) which discussed twisted trilayer graphene with commensurate twist angles and Cea et al. (2019) which discusses several multilayer settings including a related one with equal rather than alternating twist angles.

Acknowledgements.

We thank Bertrand Halperin, Philip Kim and Pablo Jarillo-Herrero for fruitful discussions. A.J.K. was supported by the Swiss National Science Foundation’s grant P2ELP2_175278. G.T. was supported by the MURI grant W911NF-14-1-0003 from ARO and by DOE grant de-sc0007870 and DOE Grant No. DE-SC0019030. A.V. and E.K. were supported by a Simons Investigator award and by nsf-dmr 1411343.

Appendix A Derivation of the Hamiltonian

Our setting is a generalization of the twisted bilayer problem considered in Refs. Bistritzer and MacDonald (2011); Dos Santos et al. (2007). We consider $n$ graphene layers where the $\ell$-th layer is rotated counter-clockwise by an angle $\theta_{\ell}$ then displaced by distance ${\bm{d}}_{\ell}$ relative to a fixed reference and restrict ourselves to coupling only between next neighboring layers which can, in general, depend on the layer index $\ell$. Following Bistritzer and Macdonald Bistritzer and MacDonald (2011), the coupling between layers is given by the term

[TABLE]

where $c_{\ell}({\bm{r}})$ represents the annihilation operator for an electron in the $\ell$-th layer at position ${\bm{r}}$. The electron operator in the $\ell$-th layer can be expanded in terms of the graphene orbitals living on a honeycomb lattice which is twisted by angle $\theta_{\ell}$ and displaced by ${\bm{d}}_{\ell}$

[TABLE]

Here, $\phi_{{\bm{R}},\alpha}$ is the orbital centered at unit cell ${\bm{R}}=m_{1}{\bm{a}}_{1}+m_{2}{\bm{a}}_{2}$ in sublattice $\alpha$, $R_{\theta}$ is the rotation matrix $e^{-i\theta\sigma_{y}}$, and $\tau_{\alpha}$ is the real space position of the sublattice $\alpha$ given by $\tau_{A/B}=(\frac{\sqrt{3}}{2},\pm\frac{1}{2})a$. We now substitute in (19) and assume that the orbital $\phi_{{\bm{R}},\alpha}({\bm{r}})$ is strongly localized in space such that $\phi_{{\bm{R}},\alpha}({\bm{r}})\approx\delta({\bm{r}}-{\bm{R}}-\tau_{\alpha})$ yielding

[TABLE]

Next, we assume $t^{ij}({\bm{r}},{\bm{r}}^{\prime})$ only depends on $|{\bm{r}}-{\bm{r}}^{\prime}|$ and introduce the Fourier transform

[TABLE]

In addition, we expand the operator $f_{i,{\bm{R}},\alpha}$ in terms of the annihilation operator for the Bloch states

[TABLE]

Substituting in (21) and using $t_{{\bm{q}}}=t_{|{\bm{q}}|}$ yields

[TABLE]

Here, the momentum ${\bm{p}}({\bm{p}}^{\prime})$ is measured relative to the Brillouin zone center in the $i$-th ($j$-th layer).

Next, we consider the limit when both ${\bm{p}}$ and ${\bm{p}}^{\prime}$ are close to the $K$ point and the rotation angles $\theta_{i}$ and $\theta_{j}$ are small. In this case, we can restrict the sum in (24) to the largest terms which correspond to ${\bm{p}}\approx K$ and ${\bm{G}}=0,{\bm{G}}_{2},{\bm{G}}_{3}$ where ${\bm{G}}_{2,3}=\left(-\frac{2\pi}{\sqrt{3}a},\pm\frac{2\pi}{3a}\right)$. The Bloch state annihilation operator $\psi_{i,{\bm{p}},\alpha}$ in the vicinity of $K$ can be expressed in terms of the annihilation operator for real space slowly varying orbitals at the K valley as

[TABLE]

Here $\tilde{\chi}_{i,K,\alpha}({\bm{r}})$ denotes the annihilation operator for an electron in valley K, sublattice $\alpha$ at point ${\bm{r}}$ measured relative to the coordinate system of the $i$-th layer. The same operator in the reference coordinate system is given by $\chi_{i,K,\alpha}({\bm{r}})=\tilde{\chi}_{i,K,\alpha}(R_{-\theta_{i}}({\bm{r}}-{\bm{d}}_{i}))$. Substituting in (25) gives

[TABLE]

We now express the momenta close to the $K$ point in terms of the reference (unrotated Brillouin zone) ${\bm{p}}=R_{-\theta_{i}}{\bm{k}}+K$, ${\bm{p}}^{\prime}=R_{-\theta_{j}}{\bm{k}}^{\prime}+K$ to get

[TABLE]

To show that this reduces to the Bistritzer Macdonald case, we perform the gauge transformation $c_{i}\rightarrow c_{i}e^{iR_{\theta_{i}}K\cdot{\bm{d}}_{i}}$ which leads to

[TABLE]

which reduces to the expression of Bistritzer and Macdonald for $R_{\theta_{i}}=1$, ${\bm{d}}_{i}=0$. In the following, we will prefer to use the form (27) because the way displacement enters in the Hamiltonian is more transparent as we see below. We now assume that the value of $w^{ij}$ is different between the diagonal and off-diagonal terms to take into account the lattice relaxation effects Carr et al. (2019) leading to

[TABLE]

and $k_{D}=\frac{4\pi}{3\sqrt{3}a}$. The full Hamiltonian can then be written as

[TABLE]

where ${\bm{\sigma}}_{\theta_{\ell}}=e^{\frac{i}{2}\theta_{\ell}\sigma_{z}}{\bm{\sigma}}e^{-\frac{i}{2}\theta_{\ell}\sigma_{z}}$. The (first quantized) Hamiltonian can be written explicitly as

[TABLE]