# Moiré bands in twisted double-layer graphene

See allHide authors and affiliations

Contributed by Allan H. MacDonald, June 7, 2011 (sent for review December 8, 2010)

## Abstract

A moiré pattern is formed when two copies of a periodic pattern are overlaid with a relative twist. We address the electronic structure of a twisted two-layer graphene system, showing that in its continuum Dirac model the moiré pattern periodicity leads to moiré Bloch bands. The two layers become more strongly coupled and the Dirac velocity crosses zero several times as the twist angle is reduced. For a discrete set of magic angles the velocity vanishes, the lowest moiré band flattens, and the Dirac-point density-of-states and the counterflow conductivity are strongly enhanced.

Low-energy electronic properties of few layer graphene (FLG) systems are known (1–8) to be strongly dependent on stacking arrangement. In bulk graphite 0° and 60° relative orientations of the individual layer honeycomb lattices yield rhombohedral and Bernal crystals, but other twist angles also appear in many samples (9). Small twist angles are particularly abundant in epitaxial graphene layers grown on SiC (10, 11), but exfoliated bilayers can also appear with a twist, and arbitrary alignments between adjacent layers can be obtained by folding a single graphene layer (12, 13).

Recent advances in FLG preparation methods have attracted theoretical attention (14–20) to the intriguing electronic properties of systems with arbitrary twist angles, usually focusing on the two-layer case. The geometry of the bilayer system is characterized by a twist angle *θ* and by a translation vector **d**. Commensurability is determined only by *θ*. Sliding one layer with respect to the other in a commensurate structure modifies the unit cell but leaves the bilayer crystalline. In this work we find it convenient to regard the AB stacking as the aligned configuration. The positions of the carbon atoms in the two misaligned layers labeled by **R** and **R**^{′} are then related by **R**^{′} = *M*(*θ*)(**R** - *τ*) + **d**, where *M* is a 2-D rotation matrix within the graphene plane, and *τ* is a vector connecting the two atoms in the unit cell.

The problem is mathematically interesting because a bilayer forms a two-dimensional crystal only at a discrete set of commensurate rotation angles; for generic twist angles Bloch’s theorem does not apply microscopically and direct electronic structure calculations are not feasible. For twist angles larger than a few degrees the two layers are electronically isolated to a remarkable degree, except at a small set of angles which yield low-order commensurate structures (16, 19). As the twist angles become smaller, interlayer coupling strengthens, and the quasiparticle velocity at the Dirac point begins to decrease.

Here we focus on the strongly coupled small twist angle regime. We derive a low-energy effective Hamiltonian valid for any value of **d** and for *θ* ≲ 10° irrespective of whether or not the bilayer structure is periodic. We show that it is meaningful to describe the electronic structure using Bloch bands even for incommensurate twist angles and study the dependence of these bands on *θ*.

## Model

We construct a low-energy continuum model Hamiltonian that consists of three terms: two single-layer Dirac–Hamiltonian terms that account for the isolated graphene sheets and a tunneling term that describes hopping between layers. The Dirac–Hamiltonian (21) for a layer rotated by an angle *θ* with respect to a fixed coordinate system is where *v* is the Dirac velocity, **k** is momentum measured from the layer’s Dirac point, *θ*_{k} is the momentum orientation relative to the x axis, and the spinor Hamiltonian acts on the individual layer’s *A* and *B* sublattice degrees-of-freedom. We choose the coordinate system depicted in Fig. 1 for which the decoupled bilayer Hamiltonian is , where projects onto layer *i*.

We derive a continuum model for the tunneling term by assuming that the interlayer tunneling amplitude between *π*-orbitals is a smooth function *t*(*r*) of spatial separation projected onto the graphene planes. The matrix element [1]of the tunneling Hamiltonian *H*_{T} describes a process in which an electron with momentum **p**^{′} = *M***p** residing on sublattice β in one layer hops to a momentum state **k** and sublattice α in the other layer. In a *π*-band tight-binding model the projection of the wave functions of the two layers to a given sublattice are [2]and [3]Here *τ*_{A} = 0, *τ*_{B} = *τ*, and **R** is summed over the triangular Bravais lattice. Substituting Eqs. **2** and **3** in Eq. **1** and invoking the two-center approximation, [4]for the interlayer hopping amplitude in which *t* depends on the difference between the positions of the two carbon atoms we find that [5]Here Ω is the unit cell area, *t*_{q} is the Fourier transform of the tunneling amplitude *t*(**r**), the vectors **G**_{1} and **G**_{2} are summed over reciprocal lattice vectors, and . The bar notation over momenta in Eq. **5** indicates that momentum is measured relative to the center of the Brillouin zone and not relative to the Dirac point. Note that crystal momentum is conserved by the tunneling process because *t* depends only on the difference between lattice positions.*

The continuum model for *H*_{T} is obtained by measuring wave vectors in both layers relative to their Dirac points and assuming that the deviations are small compared to Brillouin-zone dimensions. The model’s utility rests centrally on the observation that, although *t*_{q} is not precisely known, it should nevertheless fall to zero very rapidly with *q* on the reciprocal lattice vector scale. This behavior follows from the property that the graphene layer separation *d*_{⊥} exceeds the separation between carbon atoms within the layers by more than a factor of 2. Because the two-center integral *t*(*r*) varies with the three-dimensional separation the strong small *r* hopping processes vary with *r* on the scale of *d*_{⊥}. For this reason *t*_{q} begins to decline rapidly for *qd*_{⊥} > 1. Fig. 2 plots *t*_{q} values obtained numerically from the *π*-band tight-binding models proposed in refs. 19, 22, and 23. The largest *t*_{q} values that enter the tunneling near the Dirac point have *q* = *k*_{D}, the Brillouin-zone corner (Dirac) wave vector magnitude, and correspond to the three reciprocal vectors **0**, , and where the latter two vectors connect a Dirac point with its equivalent first Brillouin-zone counterparts (See Fig. 1). When only these terms are retained, [6]where *w* = *t*_{kD}/Ω is the hopping energy, [7]and *ϕ* = 2*π*/3. The three **q**_{j}’s in Eq. **6** are Dirac model momentum transfers that correspond to the three interlayer hopping processes.

For **d** = 0 and a vanishing twist angle the continuum tunneling matrix is 3*wδ*_{αA}*δ*_{βB}, independent of position. By comparing with the experimentally known (24) electronic structure of an AB stacked bilayer we set *w* ≈ 110 meV for exfoliated samples, however experiments suggest (25) that *w* may be smaller in some epitaxial graphene samples. As we show below the spectrum is independent of **d** for *θ* ≠ 0. In the following we therefore set **d** = 0.

## Results

### Moiré Bloch Bands.

In the continuum model hopping is local and periodic, allowing Bloch’s theorem to be applied at any rotation angle irrespective of whether or not the bilayer is crystalline. We solve numerically for the moiré bands using the plane wave expansion illustrated in Fig. 1. Convergence is attained by truncating momentum space at lattice vectors of the order of *w*/*ℏv*. The dimension of the matrix, which must be diagonalized numerically, is roughly ∼10 *θ*^{-2} for small *θ* (measured in degrees), compared to the ∼10^{4} *θ*^{-2} matrix dimension of microscopic tight-binding models (14, 16).

Up to a scale factor the moiré bands depend on a single parameter *α* = *w*/*vk*_{θ}. We have evaluated the moiré bands as a function of their Brillouin-zone momentum **k** for many different twist angles; results for *w* = 110 meV, and *θ* = 5°, 1.05°, and 0.5° are summarized in Fig. 3. For large twist angles the low-energy spectrum is virtually identical to that of an isolated graphene sheet, except that the velocity is slightly renormalized. Large interlayer coupling effects appear only near the high energy van Hove singularities discussed by Andrei and coworkers (26). As the twist angle is reduced, the number of bands in a given energy window increases and the band at the Dirac point narrows. This narrowing has previously been expected to develop monotonically with decreasing *θ*. As illustrated in Fig. 3, we instead find that the Dirac-point velocity vanishes already at *θ* ≈ 1.05°, and that the vanishing velocity is accompanied by a very flat moiré band which contributes a sharp peak to the Dirac-point density-of-states (DOS). At smaller twists the Dirac-point velocity has a nonmonotonic dependence on *θ*, vanishing repeatedly at the series of magic angles illustrated in Fig. 4.

Partial insight into the origin of these behaviors can be achieved by examining the simplest limit in which the momentum-space lattice is truncated at the first honeycomb shell. Including the sublattice degree of freedom, this truncation gives rise to the Hamiltonian [8]where **k** is in the moiré Brillouin-zone and **k**_{j} = **k** + **q**_{j}. This Hamiltonian acts on four two-component spinors Ψ = (*ψ*_{0},*ψ*_{1},*ψ*_{2},*ψ*_{3}). The first (*ψ*_{0}) is at a momentum near the Dirac point of one layer and the other three *ψ*_{j} are at momenta near **q**_{j} and in the other layer. The dependence of *h*(*θ*) on angle is parametrically small and can be neglected. We have numerically verified that this approximation reproduces the velocity with reasonable accuracy down to the first magic angle (Fig. 4, *Inset*).

The renormalized velocity follows from the spectrum of the twisted bilayer. The Hamiltonian is expressed as a sum of the **k** = 0 term and the *k*-dependent term and solved to leading order in **k**.

Consider the **k** = **0** term in the Hamiltonian. We assume that has zero energy eigenstates and prove our assumption by explicitly finding these states. The zero energy eigenstates must satisfy [9]Because [10]the equation for the *ψ*_{0} spinor is *h*_{0}*ψ*_{0} = 0, i.e., *ψ*_{0} is one of the two zero energy states and of the isolated layer. The two zero energy eigenstates of then follow from Eq. **9**. Given that , the wave functions should be normalized by |Ψ|^{2} = 1 + 6*α*^{2}. The effective Hamiltonian matrix to leading order in *k* is therefore Aside from a renormalized velocity [11]the Hamiltonian is identical to the continuum model Hamiltonian of single-layer graphene. The denominator in Eq. **11** captures the contribution of the Ψ_{j}’s to the normalization of the wave function whereas the numerator captures their contribution to the velocity matrix elements. For small α, Eq. **11** reduces to the expression *v*^{⋆}/*v* = 1 - 9*α*^{2}, first obtained by Lopes dos Santos et al. (15). The velocity vanishes at the first magic angle because it is in the process of changing sign. The eigenstates at the Dirac point are a coherent combination of components in the two layers that have velocities of opposite sign.

### Counterflow Conductivity.

The distribution of the quasiparticle velocity between the two layers implies exotic transport characteristics for separately contacted layers. Consider a counterflow geometry in which currents in the two layers flow antiparallel to one another. We focus on twist angles *θ* ≳ 2° for which the eight-band model is valid and to the semiclassical regime in which *ϵ*_{F}*τ* > 1 and find the counterflow conductivity *σ*_{CF}. We assume that the Fermi momentum is much smaller than *k*_{θ} and that 1/*τ*_{0} < *ℏvk*_{θ}, where *τ*_{0} is single particle lifetime. Using the Kubo formula we find that [12]where [13]is the x component of the counterflow velocity operator (we set the electric fields along the x axis), is the retarded Green function with μ labeling the two Dirac bands, and is the energy dispersion of the twisted bilayer at small momenta. For an electron-doped system the valence band can be neglected and [14]where *ν*^{⋆} is the DOS of the twisted bilayer. The vertex function [15]where *v*_{CF} = *v*(1 + 3*α*^{2})/(1 + 6*α*^{2}) follows directly from the previous section if we notice the sign differences between the counterflow velocity operator and the normal velocity operator. The counterflow conductivity is then [16]where *σ*_{0} ∼ *e*^{2}*ϵ*_{F}*τ*/*π* is the conductivity of an isolated single graphene layer. As *θ* is reduced from a large value toward 1°, *v*^{⋆} is reduced and the DOS is correspondingly increased. The counterflow conductivity is enhanced because of an increased density of carriers, which is not accompanied by a decrease in counterflow velocity. For a conventional measurement in which the current in the bilayer is unidirectional *v*_{CF} in Eq. **16** is replaced by *v*^{⋆}. The increase in the DOS is then exactly compensated by the reduction of the renormalized velocity and the single-layer value of the conductivity is regained.

### Dependence of the Spectrum on **d**.

We now show that the spectrum of misaligned bilayers is independent of linear translations of one layer with respect to the other using a unitary transformation that makes the Hamiltonian independent of **d**. Consider *H*_{Q} where **Q** is a momentum in the first moiré Brillouin zone. With each momentum on the *k*-space triangular Bravais lattice (see Fig. 1) where and , we associate the phase The phase associated with momentum on the other sublattice is *ϕ*_{k} as well. In terms of the new basis states the Hamiltonian *H*_{Q} is **d**-independent.

Physically, the lack of dependence on **d** can be understood by noticing that varying **d** just shifts the moiré pattern in space. The bilayer spectrum does depend on **d** at *θ* = 0, and at other commensurate angles. We expect that dependence on **d** will be observable only at short period (large *θ*) commensurate angles.

## Discussion

Twisted double-layer graphene is, for most values of *θ*, a quasi-periodic structure that has no unit cell. Nevertheless, we find that for *θ* ≲ 10° it is meaningful to describe the electronic structure of the system in terms of Bloch bands. The hidden periodic structure is shown to be related to the moiré pattern of the overlaid layers (27).

The leading corrections to the periodic moiré band Hamiltonian involve hopping amplitudes with the smallest momenta **g** that satisfy the crystal momentum conservation condition in Eq. **5** and are larger than *k*_{D}. As we showed in ref. 19, real space commensuration between the two rotated hexagonal lattice is concomitant to momentum space commensuration of the Dirac points in the extended-zone scheme (see figures 2 and 3 in ref. 19). The commensurate vector **g** can therefore be found using the formula for the moiré periodicity if the lattice vector of graphene (where *a* is the carbon-carbon distance in a single-layer graphene) is replaced by the reciprocal lattice vector *G*. It follows that *g*(*θ*) ≈ *G*/*θ*. For example *g*(10°) = 24/*a* and *g*(2°) = 120/*a*. As Fig. 2 demonstrates, the hopping amplitudes for these large wave vectors are indeed negligible compared to the value of *t*_{kD}. We therefore expect the continuum model to be very accurate up to energies of approximately 1 eV and up to angles of approximately 10°.

The Bloch band model has a simple and appealing physical interpretation. The hopping Hamiltonian is local in space. At each position, its 4 × 4 matrix, describes sublattice-dependent interlayer hopping, which depends on the local coordination between the atoms in the two layers. In Fig. 5 we have plotted the moiré pattern of atomic positions and the smaller of the two positive eigenvalues of the hopping Hamiltonian as a function of position on the same length scale. At each position, the local interlayer tunneling Hamiltonian, is that of a system in which the local coordination is maintained through all space. At AB and BA points, for example, the tunneling Hamiltonian is that of AB and BA systems, for which tunneling does not produce a gap so that the smallest positive eigenvalue vanishes. On the other hand the gap reaches its maxima (6*w*) at AA points in the moiré pattern.

In summary we have formulated a continuum model description of the electronic structure of rotated graphene bilayers. The resulting moiré band structure can be evaluated at arbitrary twist angles, not only at commensurate values. We find that the velocity at the Dirac point oscillates with twist angle, vanishing at a series of magic angles which give rise to large DOS and to large counterflow conductivities. Many properties of the moiré bands are still not understood. For example, although we are able to explain the largest magic angle analytically, the pattern of magic angles at smaller values of *θ* has so far been revealed only numerically. Additionally the flattening of the entire lowest moiré band at *θ* ≈ 1.05° remains a puzzle. Interesting new issues arise when our theory is extended to graphene stacks containing three or more layers. Finally, we remark that electron-electron interactions neglected in this work are certain to be important at magic twist angles in neutral systems and could give rise to counterflow superfluidity (28, 29), flat-band magnetism (30), or other types of ordered states.

## Acknowledgments

We acknowledge a helpful conversation with Gene Mele. This work was supported by Welch Foundation Grant F1473 and by the National Science Foundation-Nanoelectronics Research Initiative South West Academy of Nanoelectronics program.

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. E-mail: macd{at}physics.utexas.edu.

Author contributions: R.B. and A.H.M. designed research; R.B. performed research; and R.B. and A.H.M. wrote the paper.

The authors declare no conflict of interest.

↵

^{*}A closely related but slightly different expression appears in ref. 19 in which we chose the origin at a honeycomb lattice point. The present convention is more convenient for the discussion of small rotations relative to the Bernal arrangement.

## References

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Schmidt H,
- Luedtke T,
- Barthold P,
- Haug RJ

- ↵
- Dobrik G,
- Tapaszto L,
- Nemes-Incze P,
- Lambin Ph,
- Biro LP

- ↵
- ↵
- ↵
- ↵
- ↵
- Mele EJ

- ↵
- ↵
- Morell ES,
- Correa JD,
- Vargas P,
- Pacheco M,
- Barticevic Z

- ↵
- ↵
- Pereira VM,
- Castro Neto AH,
- Peres NMR

- ↵
- ↵
- ↵
- ↵
- Li G,
- et al.

- ↵
- Amidror I

- ↵
- ↵
- Min H,
- Bistritzer R,
- Su JJ,
- MacDonald AH

- ↵

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Physics