An exact continuum model for low-energy electronic states of twisted bilayer graphene

TL;DR

This paper presents a comprehensive continuum model for twisted bilayer graphene that includes lattice relaxation effects, accurately capturing low-energy electronic states and their evolution with twist angle, improving upon previous models.

Contribution

The authors develop a complete physical continuum model incorporating lattice relaxation, providing a more accurate description of low-energy states in twisted bilayer graphene.

Findings

Lattice relaxation significantly alters the bandstructure at the first magic angle.

The model eliminates the second magic-angle twist.

Minimal models can be adapted to reflect electronic changes with twist angle.

Abstract

We introduce a complete physical model for the single-particle 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 show that minimal models for the low-energy states of tBLG can be easily modified to capture the changes in electronic states as a function of twist angle.

Full text

††preprint: APS/123-QED

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 of the nature of unconventional superconductivityCao et al. (2018a, b). The challenge in modeling these phenomena from an atomistic perspective is that the actual structure of tBLG near the magic-angle twist ($\sim 1.1^{\circ}$) 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 layersLi et al. (2009); Brihuega et al. (2012); Luican et al. (2011); Wong et al. (2015).

Existing models based on DFT calculations Trambly de Laissardière et al. (2010); Uchida et al. (2014) or large supercell tight-binding Hamiltonians Suárez Morell et al. (2010); Nam and Koshino (2017); Angeli et al. (2018) 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 Mele (2010); Bistritzer and MacDonald (2011); San-Jose et al. (2012); Weckbecker et al. (2016). 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 bandstructure Yoo et al. (2018); Dai et al. (2016); Nam and Koshino (2017); Zhang and Tadmor (2018); Lin et al. (2018). Many simplified models for the flat bands of magic-angle tBLG have been proposed based on symmetry analysis, but they rely on empirical parameterization and are designed for only the magic-angle twist configurationPo et al. (2018); Yuan and Fu (2018); Koshino et al. (2018), typically ignoring atomic relaxation.

Here, we present an ab initio $k\cdot 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 physically meaningful prediction of correlation effects, as the presence of unphysical features in the single-particle band structure causes uncontrolled error in many-body calculations. We draw new conclusions on the low-energy 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. For reference, we compare our continuum model to the seminal and widely employed $k\cdot p$ model of Bistritzer and MacDonaldBistritzer and MacDonald (2011) (BMD in the following), and we adopt their dimensionless parameter $\alpha=\omega/v_{F}k_{\theta}$ for describing the twist-angle $\theta$, where $v_{F}$ is the Fermi velocity, $k_{\theta}$ is the wave-vector set by the moiré length scale and $\omega$ is their effective interlayer coupling strength ($0.11$ eV).

Within $k\cdot 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 Val . Here we introduce an expanded ab initio $k\cdot p$ model which gives a more complete physical picture of the tBLG system. Our model has three new key ingredients:

(1) relaxation of the bilayer system Carr et al. (2018a), 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\cdot p$ continuum model, which are necessary to capture the changes in stacking order at small angles;

(3) inclusion of $k$-dependent terms, which allow the $k\cdot p$ model to reproduce more accurately the particle-hole asymmetry of realistic ab initio bandstructures.

The $k\cdot p$ terms are directly computed from an ab initio tight-binding Hamiltonian model Fang and Kaxiras (2016); Fang et al. (2018); Carr et al. (2018b) for supercells spanning the twist-angle range $0.18^{\circ}\leq\theta\leq 6^{\circ}$. These terms have smooth dependence on $\theta$, 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. 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 Fan .

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 $\theta$ smaller than a critical value $\theta_{c}\approx 1^{\circ}$, the local atomic structure near the $AA$ and $AB$ stackings of the two layers becomes independent of $\theta$. This creates a pattern of small circular domains of $AA$ stacking and large triangular domains of $AB/BA$ stackings. Domain walls ($DW$) 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 $\theta<\theta_{c}$, where the tBLG system consists of a few fixed elementsNam and Koshino (2017); Zhang and Tadmor (2018); Carr et al. (2018a), and only their length scale changes for decreasing twist angle. These elements are: the $AA$ regions which have a local twist of $\theta_{AA}=1.7^{\circ}$, which is independent of the overall twist angle $\theta$ between the two layers, the $AB$ and $BA$ regions with negligible local twist. Moreover, the diameter of the $AA$ regions and the width of the $DW$ regions are approximately equal and remain unchanged for $\theta<\theta_{c}$Zhang and Tadmor (2018). These features are shown in Fig. 1 for $\theta=0.9^{\circ}$.

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 Å to 3.59 Å, a substantial change ($>7$%). The reduction in the size of $AA$ stacking can be understood as a minimization of planar stress energy and stacking energies, and has been modeled through various methods Dai et al. (2016); Nam and Koshino (2017); Zhang and Tadmor (2018); Carr et al. (2018a) leading to a relaxed pattern in agreement with experimental results Alden et al. (2013); Yoo et al. (2018). 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 constant interlayer distance equal to the average of $AA$ and $AB$ interlayer distances (3.47 Å). The magic angle predicted by the fully relaxed model, $\theta_{c}\approx 1.0^{\circ}$, is closer to the angles where correlated phenomena are observed Cao et al. (2018a, b); Yankowitz et al. (2019).

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 10 bands Po et al. (2018), and we argue that it can accurately capture 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}\pm ip_{y})$-like character ($AA_{\pm}$), 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 $10$-band tight-binding Hamiltonian are provided in the supplementary materials.

To compare our ab initio $k\cdot p$ results to the $10$-band model, we project non-orthogonal wavefunctions that satisfy the symmetry conditions, shown in Fig. 1, from band structure calculations. The form of these wavefunctions is not sensitive to the twist-angle, and is robust for twist angles within $\pm 0.2^{\circ}$ of the magic angle. We note that the $z$ and $\pm$ indexing of the $AA$ orbitals describe their symmetry properties over the moiré supercell, not their composition in terms of atomic-scale C $p_{z}$ orbitals. We also fit the parameters of the $10$-band tight-binding model for $\theta\in[0.8^{\circ},1.8^{\circ}]$, to reproduce the bands produced by our continuum model (see SM). The flat bands near the magic angle have $AA$ and $DW$ character (see Fig. 2a), 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 $\Gamma$ flips as one reduces the twisting angle: the hole band has $DW$ character for $\theta>\theta_{c}$ and switches to $AA_{z}$ character for $\theta<\theta_{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.

Two other important parameters in the $k\cdot p$ model are the effective interlayer coupling between orbitals of the same sublattice label, $A\to A$ or $B\to B$, and that between orbitals of different labels, $A\to B$ or $B\to 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 $\theta$. 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. 2c for the Full and the Flat relaxation models. The overall $\theta$ 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 $\theta$: the Fermi velocity ($v_{F}$), the bandwidth ($E_{w}$), and the band gap ($E_{g}$). These are shown in Fig 3. 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 more regular dependence on $\theta$ and 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 $\approx 0.1^{\circ}$ 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 $\approx 0.1^{\circ}$ the flat-band models may still be reliable enough to explain correlation effects. This range for the Full relaxed model is $\theta\in[0.95^{\circ},1.05^{\circ}]$ and $\theta\in[0.80^{\circ},0.90^{\circ}]$ for the Flat model. The bandstructures for both models are similar after accounting for this offset in $\theta$.

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 $\Gamma$ point its energy energy is higher than the Fermi level (see Fig. 3). Thus the charge neutrality point does not occur at the Dirac point energy. This effect persists in all of our ab initio $k\cdot p$ models (even without relaxation), and is a behavior that can be observed in other tight-binding models in the literatureSboychakov et al. (2015); Nam and Koshino (2017); Angeli et al. (2018). Assuming the bands of tBLG are not perfectly particle-hole symmetric, and that the flat-band regime is defined by a protected $\theta$-tuned band inversion, such a feature is unavoidable. For transport measurement, this behavior would result in a range of $0.1^{\circ}$ 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 $\Gamma$. Thus if a clean Dirac-point transport signature is used to assess experimental device quality, this angle-range will be difficult to observe.

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$ Bistritzer and MacDonald (2011). In Fig. 4 we show the Fermi velocity as predicted from the BMD model and from our unrelaxed and fully relaxed ab initio $k\cdot p$ models. Although our unrelaxed model shows similar behavior to the BMD model with a second magic angle occurring near $\theta=0.5^{\circ}$, the inclusion of atomic relaxation removes this feature in near $0.5^{\circ}$. As the lattice relaxation in tBLG becomes increasingly sharp on the moiré length scale as the twist angle decreases Dai et al. (2016); Nam and Koshino (2017); Zhang and Tadmor (2018); Carr et al. (2018a); Yoo et al. (2018), these sharper features in the relaxation introduce additional important couplings in the $k\cdot p$ model at larger momenta. Thus to accurately model the electronic structure of tBLG below $1^{\circ}$ our inclusion of the higher-order $k\cdot p$ couplings terms is necessary.

In conclusion, we have presented a $k\cdot p$ expansion of the low-energy electronic states of tBLG that can be extended to arbitrary order in pertubation 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 publicly available in MATLAB, C++, and Python at https://github.com/stcarr/kp_tblg.