# Moiré bands in twisted double-layer graphene

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 knownBernalHexagonal1 (); Ando (); CastroNeto (); BernalHexagonal2 (); BernalHexagonal3 (); BernalHexagonal4 (); FLGOptical (); FLGthermal () to be strongly dependent on stacking arrangement. In bulk graphite and relative orientations of the individual layer honeycomb lattices yield rhombohedral and Bernal crystals, but other twist angles also appear in many samplesrotatedGraphite (). Small twist angles are particularly abundant in epitaxial graphene layers grown on SiCrotatedEpitaxial (), but exfoliated bilayers can also appear with a twist, and arbitrary alignments between adjacent layers can be obtained by folding a single graphene layerfoldedBilayer (); Biro ().

Recent advances in FLG preparation methods have attracted theoretical attentionSantos (); Shallcross (); ShallcrossLong (); Mele (); FLGus (); localization () to the intriguing electronic properties of systems with arbitrary twist angles, usually focusing on the two-layer case. 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 structuresShallcrossLong (); FLGus (). 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 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 HamiltonianGrapheneReview () for a layer rotated by an angle with respect to a fixed coordinate system is

where is the Dirac velocity, is momentum magnitude measured from the layer’s Dirac point, and the spinor Hamiltonian acts on an individual layer’s sublattice degree-of-freedom. We choose the coordinate system depicted in Fig.1 for which the decoupled bilayer Hamiltonian is where projects onto layer .

We derive a continuum model for the tunneling term by assuming that the inter-layer tunneling amplitude between -orbitals is a smooth function of separation projected onto the graphene planes. A -band tight-binding calculation then yields (see Supplementary Information) a continuum limit in which tunneling is local. We find that the amplitude for an electron residing on sublattice in one layer to hop to sublattice on the other layer is

(1) |

where

(2) |

and . The three ’s in equation (1), depicted in Fig.1, are Dirac model momentum transfers that correspond to the three interlayer hopping processes. For and a vanishing twist angle the continuum tunneling matrix is , independent of position. By comparing with the experimentally knownKuzmenko () electronic structure of an AB stacked bilayer we set . The spectrum is independent of for . (see Supplementary Information). In the following we therefore set .

The continuum model hopping Hamiltonian captures the local stacking sequence of the misaligned layers. At , for example, only the , element of is non-zero because we have chosen our origin at a Bernal stacking point. The local lattice registration changes periodically in space with a pattern of AB points at which only is not zero, BA points at which only is non-zero and AA points at which only is non-zero. The periodicity is controlled by the reciprocal lattice vectors and . The magnitude of the moiré lattice vector is therefore , where is the carbon-carbon distance in grapheneMoire (). Because translational symmetry in the continuum model is broken only by the spatially periodic hopping Hamiltonian, we can apply Bloch’s theorem to obtain energy bands at any twist angle , independent of whether the underlying lattice is commensurate. We expect the continuum model to be accurate up to energies of 1eV and up to angles of . 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 order of . The dimension of the matrix which must be diagonalized numerically is roughly for small (measured in degrees), compared to the matrix dimension of microscopic tight-binding modelsShallcrossLong (); localization ().

Up to a scale factor the moiré bands depend on a single parameter . We have evaluated the moiré bands as a function of their Brillouin-zone momentum for many different twist angles; results for and are summarized in Fig.2. 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 inter-layer coupling effects appear only near the high energy van Hove singularities discussed by Andrei and co-workersvanHove_Andrei (). 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.2, we instead find that the Dirac-point velocity vanishes already at , 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 non-monotonic dependence on , vanishing repeatedly at the series of magic angles illustrated in Fig.3.

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 gives rise to the Hamiltonian

(3) |

where is in the moiré Brillouin-zone, and . This Hamiltonian acts on four two-component spinors; the first () is at a momentum near the Dirac point of one layer and the other three are at momenta near and in the other layer. The dependence of 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 (see inset of Fig.3). It is possible to demonstrate analytically (see Supplementary Information) that this Hamiltonian has two zero energy states at , and a Dirac velocity renormalized by

(4) |

The denominator in equation (4) captures the contribution of the ’s to the normalization of the wave function while the numerator captures their contribution to the velocity matrix elements. For small , equation (4) reduces to the expression , first obtained by Santos et al.Santos (). 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!

The distribution of the quasiparticle velocity between the two layers implies exotic transport characteristics for separately contacted layers. Consider a counterflow geometry in which the currents in the two layers flow anti-parallel to one another. The counter-flow velocity remains finite at the magic angle when the bands flatten and the density of states is enhanced. A Kubo formula calculation of the counterflow conductivity then yields (see Supplementary Information)

(5) |

where is the conductivity of an isolated single graphene layer. As is reduced from a large value towards , is reduced and the density-of-states 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 in expression (5) is replaced by . 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.

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 densities-of-states and 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
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 is this work are certain to be
important at magic twist angles in neutral systems and could give rise to counterflow superfluidityBECbilayers (); counterflowSC (),
flat-band magnetismflatBandFerromagnetism (), or other types of ordered states.

The authors acknowledge a helfpul conversation with Gene Mele. This work was supported by Welch Foundation grant F1473 and by the NSF-NRI SWAN program.

## References

- (1) McCann, E., Falko, V. I. Landau-Level Degeneracy and Quantum Hall Effect in a Graphite Bilayer. Phys. Rev. Lett. 96, 086805 (2006).
- (2) Koshino, M., Ando, T. Orbital diamagnetism in multilayer graphenes: Systematic study with the effective mass approximation. Phys. Rev. B 76, 085425 (2007).
- (3) Nilsson, J., Castro Neto, A. H., Guinea, F., Peres, N. M. R. Electronic properties of bilayer and multilayer graphene. Phys. Rev. B 78, 045405 (2008).
- (4) Min, H. MacDonald, A. H. Chiral decomposition in the electronic structure of graphene multilayers. Phys. Rev. B 77, 155416 (2008).
- (5) Koshino, M., McCann, E. Trigonal warping and Berrys phase N in ABC-stacked multilayer graphene Phys. Rev. B 80, 165409 (2009).
- (6) Zhang, F., Sahu, B., Min, H., MacDonald, A.H. Band structure of ABC-stacked graphene trilayers. Phys. Rev. B 82, 035409 (2010);
- (7) Wang, Y. et. al. Stacking-Dependent Optical Conductivity of Bilayer Graphene. ACS Nano, 4 (7), 40744080 (2010).
- (8) Ghosh, S. et al. Dimensional crossover of thermal transport in few-layer graphene. Nature Mat. 9, 555-558 (2010).
- (9) Rong, Z. Y. Kuiper, P. Electronic effects in scanning tunneling microscopy: Moiré pattern on a graphite surgace. Phys. Rev. B , 17427-17431 (1993).
- (10) Hass, J., et al. Why Multilayer Graphene on 4H-SiC Behaves Like a Single Sheet of Graphene. Phys. Rev. Lett. , 125504 (2008).
- (11) Schmidt, H., Luedtke, T., Barthold, P. Haug, R. J. Mobilities and scattering times in decoupled graphene monolayers. Phys. Rev. B 81, 121403(R) (2010).
- (12) Dobrik, G., Tapaszto, L., Nemes-Incze, P., Lambin, Ph., Biro, L. P. Crystallographically oriented high resolution lithography of graphene nanoribbons by STM lithography. Phys. Status Solidi B, 1-7 (2010).
- (13) Trambly de Laissardie‘re, G., Mayou, D., Magaud, L. Localization of Dirac Electrons in Rotated Graphene Bilayers. Nano Lett. 10, 804-808 (2010).
- (14) Lopes dos Santos, J. M. B., Peres, N. M. R., Castro Neto, A. H. Graphene Bilayer with a Twist: Electronic Structure. Phys. Rev. Lett. , 256802 (2007).
- (15) Shallcross, S., Sharma, S., Kandelaki, E., Pankratov, O. A. Electronic structure of turbostratic graphene. Phys. Rev. B , 165105 (2010).
- (16) Shallcross, S., Sharma, S., Pankratov, O. A. Quantum Interference at the Twist Boundary in Graphene. Phys. Rev. Lett. , 056803 (2008).
- (17) Mele, E. J. Commensuration and interlayer coherence in twisted bilayer graphene. Phys. Rev. B 81, 161405(R) (2010).
- (18) Bistritzer, R., MacDonald, A. H. Transport between twisted graphene layers. Phys. Rev. B 81, 245412 (2010).
- (19) Castro Neto, A. H., Guinea, F., Peres, N. M. R., Novoselov, K. S., Geim, A. K. The electronic properties of graphene. Rev. Mod. Phys. 81, 109-162 (2009).
- (20) Kuzmenko, A. B., Crassee, I., van der Marel, D., Blake, P., Novoselov, K. S. Determination of the gate-tunable bandgap and tight-binding parameters in bilayer graphene using infrared spectroscopy. Phys. Rev. B 80, 165406 (2009).
- (21) Amidror, I. The Theory of the Moire Phenomenon Vol. I (Springer, London 2009).
- (22) Li, G. et al. Observation of Van Hove singularities in twisted graphene layers. Nature Phys. , 109-113 (2010).
- (23) Miller, D. L. et al. Observing the Quantization of Zero Mass Carriers in Graphene. Science , 924-927 (2009).
- (24) Zhou, S. Y., Gweon, G.-H., Lanzara, A. Low energy excitations in graphite: The role of dimensionality and lattice defects. Ann. Phys.(N.Y.) , 1730-1746 (2006).
- (25) First, P. N. et al. Epitaxial Graphenes on Silicon Carbide. MRS Bulliten, 35, 296-305 (2010).
- (26) Eisenstein, J. P., MacDonald, A. H. BoseEinstein condensation of excitons in bilayer electron systems. Nature 432, 691 (2004).
- (27) Min, H., Bistritzer, R., Su, J. J., MacDonald, A. H. Room-temperature superfluidity in graphene bilayers. Phys. Rev. B 78, 121401(R) (2008).
- (28) Arita, R., Suwa, Y., Kuroki, K., Aoki, H. Gate-Induced Band Ferromagnetism in an Organic Polymer. Phys. Rev. Lett. 88, 127202 (2002).

Supplementary Information

## Appendix A Tunneling matrix

The tunneling matrix describes a process in which an electron with momentum residing on sublattice in one layer hops to a momentum state and sublattice in the other layer. A tight binding calculation shows that

(6) |

where is the unit cell area, linearly displaces one layer relative to the other, is the Fourier transform of the tunneling amplitude , is the vector connecting the two atoms within a unit cell, , and with being the Dirac momentum. In equation (6) we have placed the origin at a point where the sublattice of one layer lies above the sublattice of the other layer when , the sums are over reciprocal lattice vectors and the primed vectors and have been transformed by the rotation matrix . Note that the crystal momentum is conserved by the tunneling process because depends only on the difference between lattice positions. A closely related but slightly different expression appears in Ref.FLGus () 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.

The continuum model version for is obtained by measuring wavevectors in both layers relative to their Dirac points and assuming that the deviations are small compared to Brillouin-zone dimensions. vanishes rapidly with on a Brillouin-zone scale (this is partly because the inter-layer distance is larger than the honeycomb lattice constant). Therefore only values of for which contribute significantly to . The three possible values of in Eq.(6) are therefore where the latter two vectors connect a Dirac point with its equivalent first Brillioun zone counterparts. (See Fig.1.) Summing the three terms we find equation (1) in the main text.

## Appendix B Independence of the spectrum on

Here we 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 . Consider where is a momentum in the first moiré Brillouin zone. With each momentum on the k-space triangular sublattice (see Fig.1)

where , and we associate the phase

The phase associated with momentum on the other sublattice is as well. In terms of the new basis states the Hamiltonian is -independent.

## Appendix C Renormalized velocity in the 8-band model

For twist angles the 8-band Hamiltonian , given by Eq.(3) in the main text, adequately accounts for the low energy spectrum. The renormalized velocity follows from the spectrum of the twisted bilayer. We find the spectrum perturbatively in . The Hamiltonian is expressed as a sum of the term and the k-dependent term

and solved to leading order in .

Consider the term in the Hamiltonian. The wave functions of are expressed as each being a two component spinor. We now assume that has zero energy eigenstates and prove our assumption by explicitly finding these states. The zero energy eigenstates must satisfy

(7) |

Since

(8) |

the equation for the spinor is

(9) |

i.e. is one of the two zero energy states of the isolated layer. The two zero energy eigenstates of then follow from equations (7) and (9). Given that the normalization of is given by

(10) |

with being the identity matrix. For the second term in equation (10) we use the fact that is hermitian and the relations

(11) | |||||

to obtain

and .

To leading order in the energies are therefore given by the eigenvalues of

Using familiar Pauli matrix identities and equations (8,11) we obtain

Aside from a renormalized velocity

is therefore identical to the continuum model Hamiltonian of single layer graphene, , where the renormalized velocity is given by equation (4) in the main text.

## Appendix D Counterflow conductivity in the 8-band model

We find the counterflow conductivity of the bilayer system. This conductivity relates a counterflow current to the electric fields that induce it; the latter being oppositely orientated in the two layers. We restrict the calculations to twist angles for which the 8-band model is valid and to the semiclassical regime in which . We assume that only charge carriers in the linearly dispersing sector of the lowest conduction band contribute to the conductivity, i.e. that the Fermi momentum is much smaller than and that where is single particle lifetime.

Using the Kubo formula we find that

(12) |

where accounts for the spin and valley degeneracies,

(13) |

is the x-component of the counterflow velocity operator(we set the electric fields along the x-axis),

(14) |

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

(15) |

where is the density of states of the twisted bilayer. The vertex function

(16) |

where 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 given by equation (5) in the main text readily follows.