Minimal model for low-energy electronic states of twisted bilayer graphene

We introduce a physically motivated minimal model for the electronic structure of twisted bilayer graphene (tBLG), which incorporates the crucial role of lattice relaxation. Our model, based on $k \cdot p$ perturbation theory, combines the accuracy of DFT calculations through effective tight-binding Hamiltonians with the computational efficiency and complete control of the twist angle offered by continuum models. The inclusion of relaxation significantly changes the bandstructure at the first magic-angle twist corresponding to flat bands near the Fermi level (the"low-energy"states), and eliminates the appearance of a second magic-angle twist. We argue that the minimal model for the low-energy states of tBLG consists of ten bands, necessary to capture the changes in electronic states as a function of twist angle. We also provide information on the nature of these bands through their wavefunctions, which is closely tied to the features of the atomic relaxation.


I. INTRODUCTION
The discovery of correlated phases in twisted bilayer graphene (TBLG) has generated much interest in this structurally and compositionally rather simple system; it has emerged as a new platform for tunable electronic correlations, and for exploring the nature of unconventional superconductivity [1,2]. The challenge in modeling these phenomena from an atomistic perspective is that the actual structure of TBLG near the magic-angle twist (∼1.1 • ) where correlated behavior is observed, consists of a large number of atoms, exceeding 10 4 . To make progress from the theoretical point of view, a minimal model is needed that can capture the essence of single-particle states near the Fermi level ("low-energy" states). Such a model should reproduce the energy spectrum as a function of their relative twist angle with reasonable accuracy and with the required fidelity in capturing the nature of low-energy states. The appearance of correlated behavior is related to bands with very low dispersion ("flat" bands) caused by interlayer hybridization between the two Dirac cones from the different layers [3][4][5][6].
Existing models based on density functional theory (DFT) calculations [7,8] or large supercell tight-binding Hamiltonians [9][10][11] are too complex to form the basis of a realistic many-body theory. At the other extreme, simplified continuum models allow for efficient calculations, but are based on heuristic arguments about the nature of the relevant electronic states [12][13][14][15]. An important feature of the physical system is the presence of atomic relaxation near the magic-angle twist, which has significant effects on the low-energy band structure [10,[16][17][18][19]. Many simplified models for the flat bands of magic-angle TBLG have been proposed based on symmetry analysis, but they rely on empirical parametrization and are designed for only the magic-angle twist configuration [20][21][22], typically ignoring atomic relaxation.
Here, we present an ab initio k·p perturbation continuum model for TBLG which accurately accounts for the effects of atomic relaxation. Our model reproduces the results of DFT-quality tight-binding Hamiltonians but at a smaller computational cost and, more importantly, it applies to all twist angles near the magic-angle value. Such a single-particle model is a prerequisite for a physically meaningful prediction of correlation effects, as the presence of unphysical features in the single-particle band structure causes uncontrolled errors in many-body calculations. We draw conclusions on the lowenergy electronic states at small twist angles, including the interesting result that there are no additional vanishings of the Fermi velocity in the range of the previously expected second and third magic angles. As a complementary perspective on how the electronic structure is affected by atomic relaxations, we point the reader to Ref. [23], which studies a single commensurate twisting angle (1.05 • ) under various empirical TBLG models. For reference, we compare our continuum model to the seminal and widely employed k·p model of Bistritzer and MacDonald [13] (BMD in the following), and we adopt their dimensionless parameter α = ω/v F k θ for describing the twist angle θ , where v F is the Fermi velocity, k θ is the wave vector set by the moiré length scale, and ω is their effective interlayer coupling strength (0.11 eV).
Within k·p perturbation theory, the set of Bloch states of the two graphene layers is augmented by the addition of interlayer couplings due to the twist-angle-induced umklapp scattering process. As the low-energy electronic structure of TBLG is dominated by a pair of Dirac cones, the momentum expansion can be carried out about one copy of the cone at a valley K point. Taking also into account spin degeneracy, each band represents four electronic states in a real system [24].
Here, we introduce an expanded ab initio k·p model which gives a more complete physical picture of the TBLG system. Our model has three key ingredients: (1) relaxation of the bilayer system [25], including the out-of-plane relaxation of different regions as well as the in-plane strain corrections to the Hamiltonian of the individual monolayers; (2) terms beyond the first shell of couplings in the k·p continuum model, which are necessary to capture the changes in stacking order at small angles; and (3) inclusion of k-dependent terms, which allow the k·p model to reproduce more accurately the particlehole asymmetry of realistic ab initio band structures.
The k·p terms are directly computed from an ab initio tightbinding Hamiltonian model [26][27][28] for supercells spanning the twist-angle range 0.18 • θ 6 • . These terms have a smooth dependence on θ , allowing for interpolation between the specific twist angles that correspond to finite supercells, to generate a model valid for any desired angle in that range.

A. Continuum expansion
Our low-energy Hamiltonian for TBLG contains 2 × 2 block-diagonal elements representing Bloch waves of individual graphene monolayers at various momenta. These blocks are then coupled to one another by in-plane strain or interlayer coupling that break the translation symmetry of the graphene unit cell. These terms have the translational symmetry of the moiré supercell up to a gauge-dependent phase, and thus can be transformed from a real-space (r) basis to a momentum (k) basis as dictated by the reciprocal lattice of the supercell. There is also an implicit truncation on the number of k included in the expansion, understood to be a momentum cutoff parameter that is taken large enough to ensure sufficient convergence of the electronic property of interest. Explicitly, the Hamiltonian takes the form H i D (k) are the Dirac Hamiltonians for each individual graphene monolayer layer while V i (r) is the external potential for each layer. This external potential can include an electric gating potential, sublattice mass terms, or a modulated electric potential from doping or charge redistribution. A i (r) is the inplane pseudogauge field coupled to the Dirac electron. These terms are generated from the geometric deformation and strain for each layer. When expanded in Fourier components with supercell reciprocal lattice vectors, it can be written as A i (r) = p i A i p i e ip i ·r , where p i = nG 1 + mG 2 for some integers m, n and supercell reciprocal lattice vectors G i .
The first part of the interlayer coupling termT (r) gives the scattering terms as in the original BMD model [13]. However, we generalize the expansion to include higher-order terms as T (r) = q iT q i e iq i ·r , where q i = K 1 − K 2 + G for K i the K point of layer i's Brillouin zone and G a reciprocal lattice vector of the supercell. The remaining parts of the interlayer coupling terms, withk ± =k x ± ik y , are the momentumdependent scattering terms. They are relevant for the particlehole asymmetric features of the tight-binding band structures. The anticommutator notation is used to symmetrize the noncommuting operatorsr andk. The coupling constants here are investigated numerically with twist-angle dependence. They are derived from the projection of tight-binding supercells with relaxed geometry obtained from elastic theory. We relegate the detailed description of the extended Hamiltonian and the procedure for obtaining the relevant terms of the continuum model to a companion paper [29].

B. Validation against tight-binding models
In Fig. 1 we compare band structures of a full tight-binding model, our continuum k·p model, and a fitted form of a ten-band model [see Ref. [20], and Supplemental Material (SM) [30]] across four representative angles. The continuum model perfectly recreates the low-energy tight-binding bands if we sample over both valleys of the monolayers. The tenband model for each angle is obtained by first computing the k·p band structure along the high-symmetry lines, and then optimizing the ten-band model's 18 parameters to minimize where E b are the eigenvalues of the ten-band model for band b and b are the precomputed eigenvalues for the k·p model. As we want the model to be most accurate near the Dirac cone energy (E = 0), we use w b to weigh the central bands higher than the outer bands during the optimization procedure.

C. Atomic relaxations
Our continuum model affords a natural interpretation of the electronic structure of TBLG at small twist angles, which is derived directly from the atomic relaxation so we describe this aspect first. For twist angle θ smaller than a critical value θ c ≈ 1 • , the local atomic structure near the AA and AB stackings of the two layers becomes independent of θ . This creates a pattern of small circular domains of AA stacking and large triangular domains of AB/BA stackings. Domain walls (DWs) of intermediate stacking separate the AB and BA domains and connect the AA regions. This creates local electronic environments which are locked in with respect to changing twist angle for θ < θ c , where the TBLG system consists of a few fixed elements [10,18,25], and only their length scale changes for decreasing twist angle. These elements are the AA regions which have a local twist of θ AA = 1.7 • , which is independent of the overall twist angle θ between the two layers, the AB and BA regions with a negligible local twist. Moreover, the diameter of the AA regions and the width of the DW regions are approximately equal and remain unchanged for θ < θ c [18]. These features are shown in Fig. 2 for θ = 0.9 • . The relaxation in TBLG is described by two simultaneous effects. In-plane relaxation decreases the area of the high stacking energy AA region while it increases that of low stacking energy AB/BA regions. Out-of-plane relaxation causes corrugation, increasing the vertical separation between the AA regions from the equilibrium distance in AB stacking of 3.35-3.59 Å, a substantial change (>7%). These relaxations are included in a tight-binding model via their Fourier components [30]. The reduction in the size of AA stacking can be understood as a minimization of planar stress FIG. 2. Left: Structure of relaxed TBLG at θ = 0.9 • with exaggerated vertical relaxation (top). AA, AB, and BA stackings, and domain walls (DWs) are labeled along with a schematic representation (bottom) of the ten orbitals per unit cell of the moiré pattern required to describe the low-energy electronic states: three at the AA region, one at each of the three DW regions, and two at each of the AB and BA regions. Right: Wave-function magnitudes, |ψ l | 2 , l = AA ± , AA z , DW, AB/BA, of the ten-band model, at θ = 0.9 • , projected in the two layers (L 1 and L 2 ) and the sublattices A and B of each layer; the total (far-right column) is the sum of all layer and sublattice contributions [30]. The underlying moiré supercell lattice is given by the thin white lines. energy and stacking energies, and has been modeled through various methods [10,17,18,25], leading to a relaxed pattern in agreement with experimental results [16,31]. The role of vertical relaxation in experimental devices is less understood, as only free-standing TBLG has been modeled. Experimental TBLG devices are typically encapsulated in hexagonal boron nitride, so the actual corrugation may be reduced compared to the free-standing case. To take this into account, we consider two limits of the vertical relaxation: a "full" relaxation model (free-standing bilayer result) and a "flat" model with a constant interlayer distance equal to the average of AA and AB interlayer distances (3.47 Å). The magic angle predicted by the fully relaxed model, θ c ≈ 1.0 • , is closer to the angles where correlated phenomena are observed [1,2,32].

A. A ten-band model
The low-energy electronic states are directly associated with and derived from the presence of the relaxation-induced structural elements described earlier. Since the discovery of correlated phases in TBLG, many simplified n-band models have been proposed for the flat bands, usually based on localized functions of the BMD model. One such minimal model consists of ten bands [20], and we argue that it can accurately capture both the band structure (Fig. 1) and the electronic effects of the different stacking regions that emerge after relaxation. This model comprises three orbitals on a triangular lattice formed by the AA sites, one of p z -like character (AA z ) and two of (p x ± ip y )-like character (AA ± ), three orbitals on a kagome lattice formed by the domain walls, and four orbitals on a honeycomb lattice, two for each of the AB and BA domains. The full details of the ten-band tight-binding Hamiltonian are provided in the Supplemental Material [30].
To compare our ab initio k·p results to the ten-band model, we project nonorthogonal wave functions that satisfy the symmetry conditions, shown in Fig. 2, from band-structure calculations. The form of these wave functions is not sensitive to the twist angle, and is robust for twist angles within ±0.2 • of the magic angle. We note that the z and ± indexing of the AA orbitals describes their symmetry properties over the moiré supercell, not their composition in terms of atomic-scale C p z orbitals. As described in Sec. II B, we have fit parameters of the ten-band tight-binding model for θ ∈ [0.8 • , 1.8 • ] to reproduce the bands produced by our continuum model. The flat bands near the magic angle have AA and DW character [see Fig. 3(a)], showing that the coupling between these states is a necessary ingredient of the model if it is to capture the electronic structure as a function of twist angle. In particular, the orbital character of the electron and hole bands at flips as one reduces the twisting angle: The hole band has DW character for θ > θ c and switches to AA z character for θ < θ c , while the electron band has the reverse character. As the AA z and DW orbitals have opposite xy-plane mirror symmetry eigenvalues (−1 and +1, respectively), the magic angle represents a symmetry-protected band inversion.

B. Twist-angle dependence
Two important parameters in the k·p model are the effective interlayer coupling between orbitals of the same sublattice label, A → A or B → B, and that between orbitals of different labels, A → B or B → A. These nearest-neighbor interlayer couplings have been labeled w i , i = 0, 1 in previous studies and have a simple geometric interpretation: w 0 is the interlayer electronic coupling at the AA sites and w 1 is the coupling at AB/BA sites, averaged over the entire moiré cell. The values of these w i parameters depend strongly on the twist angle θ . As the lattice relaxes, the relative size of the AA regions is greatly reduced while that of the AB/BA regions is increased, causing a reduction in the value of w 0 and a modest increase in the value of w 1 . This dependence is shown in Fig. 3(c) for the full and the flat relaxation models. The overall θ dependence of the ratio w 0 /w 1 is not sensitive to the relaxed height assumption. The flat model has a larger ratio as the full relaxation assumption moves the AB/BA sites closer together (increasing their coupling and the w 1 value) while moving the AA sites farther apart (reducing their coupling and the w 0 value).
To elucidate the salient features of the single-particle model, we study three related indicators of the flat-band phenomenon as a function of θ : the Fermi velocity (v F ), the bandwidth (E w ), and the band gap (E g ). These are shown in Fig. 4. All three are calculated for both the electron and the hole sides of the flat-band manifold. The model without relaxation shows large discrepancies between the extrema of the Fermi velocity, gap, and bandwidth, and the electron and hole features have little in common. The two models (flat and full) that include relaxation show a more regular dependence on θ and a closer correspondence between the electron and hole bands. The bandwidth for the hole band is always smaller than that of the electron band, and the hole band achieves its minimum twice. In general, v F = 0 does not coincide with bandwidth minima. We thus draw the important conclusion that the magic angle is not a single value, but rather a range of ≈0.1 • which spans the extrema in these key features. In particular, even if an experimental device has a variation in twisting angle over a probed region, if that variation is ≈0.1 • the flat-band models may still be reliable enough to explain correlation effects. This range for the full relaxed model is θ ∈ [0.95 • , 1.05 • ] and θ ∈ [0.80 • , 0.90 • ] for the flat model. The band structures for both models are similar after accounting for this offset in θ .
An interesting behavior of the full relaxed model occurs at the center of the magic-angle regime: Although the Dirac cone still has symmetric dispersion near the K point, the hole band dispersion is such that near the point its energy is higher than the Fermi level (see Fig. 4). Thus the charge neutrality point does not occur at the Dirac point energy. This effect persists in all of our ab initio k·p models (even without relaxation), and is a behavior that can be observed in other tight-binding models in the literature [10,11,33]. Assuming the bands of TBLG are not perfectly particle-hole symmetric, and that the flat-band regime is defined by a protected θ -tuned band inversion, such a feature is unavoidable. For transport measurement, this behavior would result in a range of 0.1 • in twist angle where the charge neutrality point of the flat bands does not align with the Dirac point of the moiré superlattice, as well as a reduction in the resistivity at the Dirac point energy due to these other bands near . Thus if a clean Dirac point transport signature is used to assess experimental device quality, this angle range will be difficult to observe.

C. Suppression of magic angles
Another important result of our calculations including atomic relaxation in TBLG is the suppression of the second magic-angle twist, defined as a smaller twist angle at which v F = 0 [13]. In Fig. 5 we show the Fermi velocity as predicted from the BMD model and from our unrelaxed and fully relaxed ab initio k·p models. Although our unrelaxed model shows similar behavior to the BMD model with a second magic angle occurring near θ = 0.5 • , the inclusion of atomic relaxation removes this feature near 0.5 • . As the lattice relaxation in TBLG becomes increasingly sharp on the moiré length scale as the twist angle decreases [10,[16][17][18]25], these sharper features in the relaxation introduce additional important couplings in the k·p model at larger momenta. Thus to accurately model the electronic structure of TBLG below 1 • our inclusion of the higher-order k·p couplings terms is necessary.

IV. CONCLUSION
We have presented a k·p expansion of the low-energy electronic states of TBLG that can be extended to arbitrary order in perturbation theory. This exact continuum model facilitates a better understanding of the single-particle features of TBLG's flat bands, and provides a solid foundation on which to build correlated models.
We have made this model openly available in MATLAB, C++, and PYTHON [34].