Twisted Multilayer Graphene Model

• Twisted multilayer graphene is a structure where rotated graphene layers form moiré superlattices, yielding flat bands and exotic topological features.
• The model employs atomistic tight-binding, continuum, and chiral-limit approaches to capture interlayer coupling, band decomposition, and symmetry effects.
• Its insights underpin research into correlated insulators, superconductivity, and quantum Hall phenomena, guiding both experimental and computational studies.

Twisted multilayer graphene (TMG) model denotes the family of atomistic, continuum, chiral-limit, and effective many-body descriptions for graphene stacks in which two or more ordered graphene blocks are relatively rotated, producing moiré superlattices and, in more general settings, moiré-of-moiré structures. In current usage, the term covers generic $(M+N)$-layer twisted systems, alternating-twist multilayers, mirror-symmetric double-twisted stacks, and small-angle models derived from flat-band projections. These descriptions are used to organize flat bands, Dirac and higher-order band crossings, valley Chern hierarchies, moiré phonons, magnetic quantization, and correlated phases within a common framework (Do et al., 2021, Zhang et al., 2020).

1. Geometry, moiré scales, and stack classifications

The geometric starting point is monolayer graphene with lattice constant $a \approx 2.46\,\text{\AA}$, Dirac momentum $K=4\pi/(3a)$, and a relative twist $\theta$ between adjacent layers or graphene blocks. For a single twisted interface, the moiré period is

$L_m=\frac{a}{2\sin(\theta/2)},$

and the characteristic momentum transfer is $k_\theta = 2K\sin(\theta/2)$. In commensurate constructions, moiré lattice vectors and reciprocal vectors can be written explicitly in terms of the rotated monolayer lattices, while for small $\theta$ the continuum limit replaces atomistic periodicity by a long-wavelength moiré Brillouin zone (Do et al., 2021).

Several stack classes recur across the literature. Generic twisted $(M+N)$-layer systems twist one ordered $M$-layer block against one ordered $N$-layer block. Alternating-twist multilayers assign layer angles $a \approx 2.46\,\text{\AA}$0, producing a single aligned moiré throughout the stack. Helical trilayers use consecutive twists of the same sign, as in $a \approx 2.46\,\text{\AA}$1. Double-twisted multilayer graphene introduces two moiré interfaces with a mirror plane, yielding stacks such as $a \approx 2.46\,\text{\AA}$2 or $a \approx 2.46\,\text{\AA}$3. Twisted double bilayer graphene (TDBG) is the AB–AB case in which two Bernal bilayers are relatively twisted across the middle interface, with interlayer spacing $a \approx 2.46\,\text{\AA}$4 (1901.10485, Ding et al., 2022).

Beyond the single-angle setting, flat-band constructions also exist when neighboring twist angles are unequal but close to coprime ratios $a \approx 2.46\,\text{\AA}$5. In that regime one obtains a fast moiré lattice of scale $a \approx 2.46\,\text{\AA}$6 and a slow “moire-of-moire” modulation of scale $a \approx 2.46\,\text{\AA}$7, so a local-displacement description becomes natural. This broadens the notion of a magic-angle point into an extended magic phase in angle-disordered multilayers and narrow turbostratic graphite samples (Foo et al., 2023).

2. Hamiltonian families and decomposition principles

At the single-particle level, TMG is not represented by one universal Hamiltonian but by a hierarchy of models adapted to scale and observable. Atomistic tight-binding descriptions retain graphene $a \approx 2.46\,\text{\AA}$8 orbitals, intralayer hoppings, and position-dependent interlayer Slater–Koster transfer integrals. Continuum descriptions generalize Bistritzer–MacDonald by coupling rotated Dirac cones with three moiré harmonics. Chiral-limit models set $a \approx 2.46\,\text{\AA}$9 and keep only inter-sublattice moiré tunneling, giving exact analytic control of flat bands. Projected low-energy models further reduce the problem to Wannier, four-band, or Hubbard descriptions when only the flat-band manifold is retained (Do et al., 2021, 1901.10485, Xu et al., 2018).

Model family	Core ingredients	Typical use
Atomistic tight binding	$K=4\pi/(3a)$0 orbitals, Slater–Koster hoppings, real-space lattice geometry	large-angle spectra, magnetic quantization, incommensurate structures
Continuum moiré model	rotated Dirac blocks plus three-harmonic $K=4\pi/(3a)$1 tunneling	small-angle minibands and valley-resolved topology
Chiral-limit continuum	$K=4\pi/(3a)$2, off-diagonal $K=4\pi/(3a)$3	exact flat bands, ideal geometry, rigorous spectral results
Projected effective model	flat-band Wannier or triangular-lattice Hubbard reduction	superconductivity, Mott physics, low-energy responses

A generic continuum Hamiltonian for a twisted $K=4\pi/(3a)$4-layer stack can be written in block form as

$K=4\pi/(3a)$5

where $K=4\pi/(3a)$6 and $K=4\pi/(3a)$7 are the chiral multilayer blocks of the untwisted stacks and $K=4\pi/(3a)$8 is the interface moiré coupling (Liu et al., 2019). For Bernal-aligned interfaces, the hybrid $K=4\pi/(3a)$9-tight-binding framework refines the minimal continuum tunneling by importing the full Slonczewski–Weiss–McClure parametrization, including $\theta$0, $\theta$1, $\theta$2, and, in multilayers, $\theta$3 and $\theta$4; this largely preserves flat-band structure but restores electron–hole asymmetry (Garcia-Ruiz et al., 2021).

A distinctive feature of TMG theory is the existence of exact decomposition rules. In alternating-twist $\theta$5-layer graphene with uniform couplings, singular-value decomposition of the layer-space ladder matrix yields decoupled bilayer channels with effective couplings

$\theta$6

plus a decoupled monolayer when $\theta$7 is odd. This produces a magic-angle hierarchy directly inherited from twisted bilayer graphene (1901.10485). For arbitrary stacking, chiral decomposition reorganizes the low-energy content into two moiré-induced flat bands per valley together with pseudospin doublets whose chirality $\theta$8 is fixed by the lengths of rhombohedral-like stacking segments away from the twisted interface (Zhang et al., 2020). In mirror-symmetric double-twisted systems, an exact mirror symmetry decomposition splits the Hamiltonian into even- and odd-parity sectors, each equivalent to either a single-twist multilayer or an untwisted multilayer. In $\theta$9-DTMLG this produces one TMBG-like and one TBG-like sector; in $L_m=\frac{a}{2\sin(\theta/2)},$0-DTMLG the even sector acquires a $L_m=\frac{a}{2\sin(\theta/2)},$1 enhancement of moiré tunneling and an enlarged magic angle (Ding et al., 2022).

3. Flat bands, higher-order crossings, and topological structure

The canonical small-angle control parameter is

$L_m=\frac{a}{2\sin(\theta/2)},$2

or its equivalent in a given normalization. In the chiral limit, the first magic coupling coincides with the chiral twisted-bilayer value, numerically $L_m=\frac{a}{2\sin(\theta/2)},$3, and the corresponding multilayer bands become exactly flat at charge neutrality. In the hierarchy of two twisted sheets of $L_m=\frac{a}{2\sin(\theta/2)},$4-layer Bernal graphene, the flat-band wavefunctions are holomorphic in momentum space, the band geometry is ideal, and the flat-band Chern number has magnitude $L_m=\frac{a}{2\sin(\theta/2)},$5 per valley, with sign set by convention (Wang et al., 2021, Yang, 2023).

A major refinement of the chiral theory is the rigorous trichotomy proved for the first two bands of twisted $L_m=\frac{a}{2\sin(\theta/2)},$6-layer Bernal sheets. Away from special couplings, the bands meet at the moiré Dirac points through an order-$L_m=\frac{a}{2\sin(\theta/2)},$7 tangential crossing. At the discrete magic set $L_m=\frac{a}{2\sin(\theta/2)},$8, the first bands are perfectly flat. At a second discrete set $L_m=\frac{a}{2\sin(\theta/2)},$9, the spectrum instead develops Dirac cones, while the remaining low-energy branch retains order-$k_\theta = 2K\sin(\theta/2)$0 tangency. This overturns the earlier expectation that chiral multilayer models interpolate only between higher-order crossings and flat bands (Li et al., 4 Aug 2025).

Topological classification in more generic $k_\theta = 2K\sin(\theta/2)$1-layer systems is often expressed as a valley Chern hierarchy. If the stacking chiralities of the bottom and top stacks are labeled $k_\theta = 2K\sin(\theta/2)$2, the total Chern number of the two low-energy interface bands in valley $k_\theta = 2K\sin(\theta/2)$3 is

$k_\theta = 2K\sin(\theta/2)$4

Opposite stacking chiralities therefore give $k_\theta = 2K\sin(\theta/2)$5, while identical stacking chiralities give $k_\theta = 2K\sin(\theta/2)$6. The sign reverses in the opposite valley (Liu et al., 2019). In arbitrary-stacking TMMG, these Chern numbers are often obscured by energetic entanglement between flat bands and pseudospin doublets, but vertical displacement fields can gap the crossings and make the flat-band valley Chern numbers well defined (Zhang et al., 2020).

Optical driving adds a further topological control parameter. In Floquet-engineered TMG, circularly polarized light produces a valley-contrasting mass

$k_\theta = 2K\sin(\theta/2)$7

and the central Floquet flat bands obey a Chern hierarchy

$k_\theta = 2K\sin(\theta/2)$8

where $k_\theta = 2K\sin(\theta/2)$9 labels conduction or valence, $\theta$0 the moiré valley, and $\theta$1 same- or opposite-chirality stacking. This makes controlled valley-polarized high-Chern flat bands possible in systems such as TDBG (Lu et al., 2020).

4. Magnetic quantization, Hall response, and collective lattice modes

A lattice-resolved treatment of magnetic field begins with the generalized Peierls substitution,

$\theta$2

In monolayer and Bernal bilayer graphene this reproduces the standard Dirac and chiral bilayer Landau levels. Applied to twisted bilayer graphene (TBG) and AB–AB TDBG, it yields a unified account of Landau levels, inter-Landau-level selection rules, and quantum Hall conductivity across twist angles. In large-angle TBG the Hall sequence is $\theta$3, reflecting doubled degeneracy relative to monolayer graphene; in AB–AB TDBG the sequence is $\theta$4 with a possible $\theta$5 jump when the $\theta$6 and $\theta$7 levels overlap; in magic-angle TBG the flat-band regime gives $\theta$8 and small $\theta$9 plateaus near zero energy because the $(M+N)$0 levels are only doubly degenerate (Do et al., 2021).

The same formalism clarifies selection rules. Monolayer graphene, Bernal bilayer graphene, and large-angle TBG are dominated by $(M+N)$1 inter-Landau-level transitions. AB–AB TDBG additionally permits $(M+N)$2 channels because the $(M+N)$3 Landau-level wavefunction contains mixed $(M+N)$4 and $(M+N)$5 oscillation modes on its dominant sublattices. In the quantum Hall calculation, these node structures are what produce the step pattern of $(M+N)$6 through the Kubo formula (Do et al., 2021).

Lattice dynamics obey a parallel symmetry logic. The mismatch symmetry identified in twisted bilayers and generalized to TMG unifies two terminologies that had often been presented as distinct: phasons and moiré phonons. At generic incommensurate twist, the mismatch symmetry is continuous but nonlocal, yielding phasons that may be diffusive. In the small-angle continuum limit, the same symmetry becomes continuous and local, and its spontaneous breaking produces propagating moiré phonons as Goldstone modes. For a generic twisted $(M+N)$7-layer system with all relative twists small, the mode count is $(M+N)$8 graphene acoustic phonons, $(M+N)$9 moiré phonons, and $M$0 higher-order phasons (Gao et al., 2022).

This symmetry analysis also exposes nontrivial multilayer mappings. Mirror-symmetric alternating-twist trilayer graphene can be decomposed into even and odd sectors; the even sector maps to TBG at an effective twist

$M$1

and carries the same two gapless moiré phonons, while the odd sector is gapped because the two moiré patterns are perfectly commensurate (Gao et al., 2022). A common misconception is that multilayer lattice soft modes are merely replicated graphene phonons; the symmetry analysis shows that moiré phonons and higher-order phasons are intrinsic collective coordinates of the twisted superstructure.

5. Correlated insulators and superconducting reductions

The many-body literature on TMG splits into several complementary reductions. One route starts from a single correlated moiré miniband on a triangular superlattice, with two valley orbitals and spin-1/2 per site, giving an approximate SU(4) Hubbard model

$M$2

At strong coupling, superexchange $M$3 favors even-parity bond pairing, and a local Hund term selects the spin-triplet, valley-singlet sector. On the triangular lattice, the favored fully gapped states are $M$4-type: a time-reversal-breaking Type A phase and a time-reversal-invariant Type B phase. These support half-quantum vortices carrying $M$5, and the Type A state has total chiral central charge $M$6 with eight co-propagating Majorana edge modes (Xu et al., 2018).

A second route focuses on correlated insulating states at half filling of the flat bands in specific TMG platforms. In twisted bilayer–monolayer graphene and AB–AB TDBG under finite displacement field, Hartree–Fock theory with projected Coulomb interactions finds spin-polarized, $M$7-broken nematic insulators with total Chern number zero. When the occupied states in each valley carry nonzero net Chern, these are quantum valley Hall states. The dominant intravalley Coulomb interaction leaves spin-polarized and valley-polarized states nearly degenerate, but atomic on-site Hubbard terms lower the energy of the spin-polarized nematic state by roughly $M$8–$M$9 meV per electron in the cases studied. Vertical magnetic fields then compete with this selection through orbital Zeeman splitting and the Středa density response, driving transitions to valley-polarized quantum anomalous Hall states at onset fields of order $N$0–$N$1 T (Zhang et al., 2021).

A third route keeps the continuum minibands but studies superconductivity from long-range charge fluctuations. In TBG, TDBG, and helical trilayer graphene, a diagrammatic KL-RPA treatment with moiré Umklapp processes and phonon-dressed screened Coulomb interactions gives markedly different critical temperatures. For TBG at $N$2, including both Umklapp and phonon dressing yields $N$3 K, while removing either drops $N$4 below $N$5 mK. In TDBG at $N$6, superconductivity appears only under displacement field and remains small, $N$7 mK in the reported window. In helical trilayer graphene at $N$8, the same framework gives $N$9 mK near the van Hove singularity, with weaker sensitivity to Umklapp than in TBG but strong dependence on phonon dressing (Long et al., 2024). A plausible implication is that the role of moiré Umklapp is stack-dependent because it tracks the degree of wavefunction inhomogeneity inside the moiré cell rather than simply the presence of narrow bands.

6. Realistic corrections, incommensurability, and computational implementations

Realistic TMG modeling requires several layers of correction beyond the minimal continuum picture. For Bernal-aligned interfaces, the full Slonczewski–Weiss–McClure parametrization adds momentum-dependent tunneling, trigonal warping, sublattice asymmetry, and next-nearest-layer couplings. In twisted bilayer graphene near $a \approx 2.46\,\text{\AA}$00, these corrections generate electron–hole asymmetry and a conduction-side gap of about $a \approx 2.46\,\text{\AA}$01 meV above the flat bands, while in twisted trilayers and related multilayers they change the miniband spectra only subtly. This has been taken as evidence that the minimal moiré tunneling model remains broadly applicable for low-energy multilayer flat-band calculations, even when a refined description is needed for quantitative asymmetries (Garcia-Ruiz et al., 2021).

Ab initio Wannier studies add another correction that minimal continuum models often omit: intrinsic crystal-field polarization and nonlinear screening. Ordered ABA and ABC multilayers exhibit out-of-plane intrinsic symmetric polarization toward the outer layers, with interlayer on-site differences up to $a \approx 2.46\,\text{\AA}$02 meV in ABA trilayers. In twisted multilayers the out-of-plane polarization is weaker, typically $a \approx 2.46\,\text{\AA}$03–$a \approx 2.46\,\text{\AA}$04 meV and generally directed inward, while the in-plane modulation at the twist boundary can reach about $a \approx 2.46\,\text{\AA}$05 meV and is maximal in AA regions. The same study finds out-of-plane dielectric constants $a \approx 2.46\,\text{\AA}$06 for twisted $a \approx 2.46\,\text{\AA}$07 stacks and $a \approx 2.46\,\text{\AA}$08 for twisted $a \approx 2.46\,\text{\AA}$09 at $a \approx 2.46\,\text{\AA}$10, with weak field dependence across $a \approx 2.46\,\text{\AA}$11–$a \approx 2.46\,\text{\AA}$12 mV/\AA. These scales are comparable to correlation and disorder energies near the magic angle, so they materially affect faithful model building (Tepliakov et al., 2021).

Incommensurate and arbitrarily stacked systems require different machinery. A momentum-space generalized Umklapp formalism expands eigenstates in Bloch states of each isolated layer and couples them whenever momenta match modulo reciprocal vectors of the participating layers. Applied to twisted trilayer graphene with $a \approx 2.46\,\text{\AA}$13, $a \approx 2.46\,\text{\AA}$14, and $a \approx 2.46\,\text{\AA}$15, it shows that the low-energy spectrum cannot be reconstructed from pairwise bilayer hybridizations alone: the full three-layer coupling shifts van Hove singularities and generates a visible effective $a \approx 2.46\,\text{\AA}$16 modulation in the local density of states (Amorim et al., 2018).

At larger scales, two complementary computational programs have emerged. High-throughput atomistic tight binding based on maximally localized Wannier functions constructs commensurate supercells for arbitrarily layered assemblies, using in-plane hoppings up to eight neighbors and angle-dependent interlayer $a \approx 2.46\,\text{\AA}$17-$a \approx 2.46\,\text{\AA}$18 couplings without adjustable parameters. This framework identifies low-dispersion candidates across bilayers, trilayers, quadrilayers, and higher alternating-twist stacks, including magic-angle shifts consistent with $a \approx 2.46\,\text{\AA}$19, golden-ratio, and large-$a \approx 2.46\,\text{\AA}$20 hierarchy arguments (Tritsaris et al., 2020). For structural relaxation and dynamics, a separate line of work builds machine-learned interatomic potentials. An Atomic Cluster Expansion potential trained on DFT-FE data across local twists, disregistries, multilayer thicknesses, and interlayer spacings achieves force RMSEs of $a \approx 2.46\,\text{\AA}$21, $a \approx 2.46\,\text{\AA}$22, $a \approx 2.46\,\text{\AA}$23, $a \approx 2.46\,\text{\AA}$24, and $a \approx 2.46\,\text{\AA}$25 eV/\AA\ for mono-, bi-, tri-, quad-, and pentalayer test sets, respectively, and remains stable in $a \approx 2.46\,\text{\AA}$26 K molecular dynamics of a $a \approx 2.46\,\text{\AA}$27-atom twisted bilayer with an interlayer interstitial. This makes it possible to close the loop between mechanical reconstruction and subsequent electronic modeling (Wang et al., 18 Jun 2025).

Taken together, these developments show that the “Twisted Multilayer Graphene Model” is best understood as a layered modeling ecosystem rather than a single Hamiltonian. Atomistic tight-binding captures registry, magnetic phases, and incommensurability; continuum and chiral models expose exact flat-band and topological structure; projected effective theories organize superconductivity and ordered phases; and ab initio or machine-learned structural models supply the polarization, screening, and relaxation data needed for quantitative fidelity.