## Angle-Dependent Ab initio Low-Energy Hamiltonians for a Relaxed Twisted Bilayer Graphene Heterostructure

Shiang Fang, 1, 2 Stephen Carr, 1 Ziyan Zhu, 1 Daniel Massatt, 3 and Efthimios Kaxiras 1, 2

1 Department of Physics, Harvard University, Cambridge, Massachusetts 02138, USA 2 John A. Paulson School of Engineering and Applied Sciences, Harvard University, Cambridge, Massachusetts 02138, USA 3 Department of Statistics, The University of Chicago, Chicago, Illinois 60637, USA

We present efficient angle-dependent low-energy Hamiltonians to describe the properties of the twisted bilayer graphene (tBLG) heterostructure, based on ab initio calculations of mechanical relaxation and electronic structure. The angle-dependent relaxed atomic geometry is determined by continuum elasticity theory, which induces both in-plane and out-of-plane deformations in the stacked graphene layers. The electronic properties corresponding to the deformed geometry are derived from a Wannier transformation to local interactions obtained from Density Functional Theory calculations. With these ab initio tight-binding Hamiltonians of the relaxed heterostructure, the low-energy effective theories are derived from the projections near Dirac cones at K valleys. For twist angles ranging from 0.7 ◦ to 4 ◦ , we extract both the intra-layer pseudo-gauge fields and the inter-layer coupling terms in the low-energy Hamiltonians, which extend the conventional low-energy continuum models. We further include the momentum dependent inter-layer scattering terms which give rise to the particle-hole asymmetric features of the electronic structure. Our model Hamiltonians can serve as a starting point for formulating physically meaningful, accurate interacting electron theories.

## I. INTRODUCTION

The physical system consisting of two layers of graphene with a small relative twist between them, referred to as twisted bilayer graphene (tBLG), has emerged as a new platform for studying correlated phases of matter since the discovery of its Mott insulator 1 and the superconducting 2 behavior. The unconventional nature of both the insulating and superconducting phases has prompted further experimental 3-9 and theoretical efforts 10-22 to better understand the physics of tBLG and related van der Waals heterostructures. Unconventional correlated phases are also observed in the heterostructures that include trilayer graphene with nearly aligned hBN substrate 23,24 and twisted double bilayer graphene 25-27 . The hypothesis is that a better understanding of the correlated many-body phases will emerge by exploring the many parameters available to tune the behavior of the stacked-layer hetersotstructures; these parameters include doping, external electric and magnetic fields, temperature and applied pressure 3 . The response of the system to changes in the parameters would serve to set constraints on theoretical models, much as was the case for the isotope effect 28 or the effect of external magnetic fields in understanding conventional superconductivity. In contrast to the conventional threedimensional crystalline solids, one unique adjustable parameter for tBLG and its relatives, is the twist angle between layers: by manipulating the corresponding moiré length scale through the twist angle, which can be controlled to exquisite precision, the characteristic kinetic and interaction energies can be varied without sacrificing the material quality, an effect referred to as 'twistronics' 29 . As a consequence, the existence of the unconventional, correlated phases depends sensitively on twist angle variations 1-3 . The flat bands that emerge in the electronic spectrum at the 'magic angle' ( ∼ 1 ◦ for tBLG) are a characteristic feature of these systems that signals the emergence of strong electron interaction effects, as first pointed out by A. McDonald and coworkers 30 .

At the single-particle theory level, accurate theoretical modeling can capture the angle-dependent effects on the electronic structure and band topology 31-33 , and thus serves as a starting point for formulating minimal models of interacting theories 34-37 . Such theoretical modeling should take into account both the mechanical relaxation and the electronic properties. In single-layer graphene, the structural deformation can affect the electronic structure by inducing pseudo-gauge fields coupled to the Dirac electrons 38-41 . In the case of tBLG, the interaction between the two misaligned layers gives rise to a modulated structural pattern 42,43 . The origin of the modulation is due to the varying local geometric configurations. Among these configurations, the Bernal stacking order as in the bulk graphite structure is most favorable energetically 43-45 . The relaxed tBLG crystal is hence determined by optimizing the stacking energy cost and elastic energy from layer deformation. Near the magic angle, the crystal relaxation can affect the electronic structure around the flat bands 42 . The single-particle gaps on both electron and hole sides of the flat bands can be accounted for by including crystal relaxation when compared with the experiments. When the twist angle is varied, the strength of atomic relaxation can be modified as well 42,43 . It is therefore important to construct electronic models that capture these angle-dependent effects.

Given the relaxed crystal structure of a tBLG, one can model the electronic properties. A numerical scheme to model such van der Waals heterostructures has to be accurate enough to make numerical predictions that

can be compared to experiment. Moreover, it should also give an intuitive physical picture that can facilitate the design of structures to enable electron correlations. The large number of atoms involved in the supercells of twisted bilayers has prevented straightforward numerical approaches. Numerical methods at different levels of sophistications from large scale density functional theory (DFT) calculations 46,47 , empirical tightbinding Hamiltonians 42,48-51 to low-energy k · p continuum theories 30,52-55 have been employed. Among these approaches, the continuum k · p is computationally efficient and allows continuous twist angle control of the electronic band structure unconstrained by commensurate conditions 55,56 ; the latter type of constraints are required for DFT and tight-binding calculations. Each approach has its strengths and weaknesses in terms of numerical accuracy and efficiency.

In our work, we adopt an ab initio multi-scale numerical approach to model a tBLG which takes fully into account crystal relaxation and the twist angle dependence. Briefly, this multi-scale numerical method is designed to combine the strengths of conventional approaches mentioned above. This framework allows us to extract the relevant mechanical and electronic properties at the microscopic length scale based on ab initio total-energy and electronic structure calculations. We thus obtain a computation scheme that allows for both a clear physical picture and efficient numerical implementation, discussed in detail in the main body of the paper. Through this approach, we are able to generalize the conventional continuum k p · theory 30 to capture all relevant band features of ab initio methods 57 . For example, one such feature is the pronounced particle-hole asymmetry in the tight-binding bands in tBLG, a theoretical prediction that remains an open question on the experiment side.

Another interesting aspect of the flat bands near the Fermi level is the topological properties of the manifold, related to the wavefunction texture in the full Brillouin Zone. In contrast with trivial atomic insulator bands which can be viewed as the Hilbert space from spatially localized orbitals, a non-trivial topological structure of the flat bands in a tBLG near the magic angle has implications on the construction of a finite n -band model near charge-neutrality 31,32 . In other heterostructures, flat bands with non-trivial Chern numbers are predicted 58,59 . The non-trivial topology in the manifold impacts the formation of exotic interacting phases. The accurate multi-scale modeling can provide a reliable numerical method to estimate the model parameters and investigate various perturbations and their effects on the topology of the flat bands.

The paper is structured as follows: In Sec. II, we describe the procedures in each step of the multiscale approach, from mechanical and electronic calculations with DFT and Wannier constructions, to abinitio tight-binding models and the derivations of effective low-energy k · p Hamiltonians based on the projection method. In Sec III, we further explore the numeri- cal implications of our effective models which include the mechanical and electronic properties, and the scaling of coupling constants with respect to the twist angle. We conclude and summarize the discussions in Sec IV, and remark on potential applications and extensions.

## II. RELAXED TBLG MECHANICAL AND ELECTRONIC STRUCTURE: A MULTI-SCALE APPROACH

A van der Waals heterostructure as exemplified in a tBLG is a complicated material system to study, with massive number of atoms involved at small twist angles. For a tBLG at around the magic angle, there are about 12,000 atoms in a twisted moiré supercell. Direct simulation of a complete heterostructure requires substantial computational resources and does not yield a clear physical picture in a straightforward manner. Here we take a different approach to van der Waals heterostructure system, by a multi-scale numerical scheme as outlined in Fig. 1. This multi-scale numerical scheme is designed to combine the accuracy of DFT calculations and the transparent physical picture of efficient low-energy k p · continuum theories, connected by a simplified tight-binding Hamiltonian for a moiré supercell. Our studies begin with the mechanical and electronic modeling for the much smaller local systems sampled from the heterostructure supercell. The full physical picture of the entire heterostructure is obtained by 'stitching' together this local information.

FIG. 1. Schematic representation of the multi-scale numerical approach for the tBLG electronic structure, which begins at the microscopic scale where the mechanical and electronic properties of the local configuration are modeled with DFT calculations. These results are employed for the construction of ab initio tight-binding Hamiltonian and low-energy effective k · p theories. The reduced n -band models are further simplifications intended to form the basis for interacting many-body theories. The blocks in blue color are the focus of the present work (included in the light-blue shaded region).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00058v1-figure-ab6c66c0-54aa-4efd-9f5a-6e07c5e678ec.jpg)

The backbone of this approach, the microscopic atomic coupling and energetics, are provided by accurate DFT calculations through both the generalized stacking fault energy (GSFE), which describes the mechanical energy

of the bilayer stacking 43-45 , and Wannier orbital coupling for the electronic properties 39,60 . With the GSFE, one can determine the relaxed geometry of a tBLG within continuum elasticity theory 43 . The ab-initio tightbinding Hamiltonian can be constructed from the extracted Wannier couplings and applied to the relaxed geometry 60 . These ab-initio tight-binding Hamiltonians can be further simplified by projections onto the lowenergy degrees of freedom to derive the effective k · p Hamiltonians, which are efficient ways of describing the electronic structure. Further simplifications of the k · p Hamiltonians can be carried out by projecting the lowest n -bands near the charge neutrality 31,57,61 . Our method allows us to extract the relevant parameters based on well-defined assumptions and enables study of the twist angle dependence. In the following sections, we give more details for each step in this multi-scale numerical approach.

## A. Mechanical Properties for Atomic Relaxation

We start by setting the conventions for a twisted bilayer graphene crystal, and describing the atomic relaxation effects. Two sheets of monolayer graphene (denoted as L1 and L2) are stacked together, with the L2 (L1) layer rotated counterclockwise (clockwise) by θ/ 2 . This small mismatch of the crystal orientations from two sheets due to this twist angle induces a long-wavelength interference pattern in the local atomic registry, alternating between locally AA, AB and BA stacking orders. We define the origin at the center of a locally AA stacking spot, which coincides with a six-fold rotational symmetry axis at the center of a hexagon in the honeycomb structure. To simplify our discussion, here we focus on the commensurate moiré supercell, which is spanned by the supercell primitive vectors R i ( i = 1 2) , and the associated supercell Brillouin zone with reciprocal lattice vectors G i ( i = 1 2) , 46 . From the perspective of each individual layer, the relative twist of L1 and L2 also rotates the reciprocal lattice vectors from each layer. The differences of the reciprocal lattice vectors from each layer give rise to the supercell reciprocal lattice vectors ( G i ).

Relaxed domains of locally AB/BA stacking are energetically favored in the twisted bilayer graphene crystal 43,62 . As a result of minimizing the additional energy due to the twist, the AB/BA spots are enlarged and become uniform stacking regions; these are separated by domain boundaries. Regions of AA stacking spots are reduced in size with the domain boundaries between neighboring AB/BA regions connecting the neighboring AA stacking spots. As the twist angle θ decreases, the area of AB stacking regions increases, but the domain wall width remains unchanged. This domain structure is stabilized beyond a critical twist angle 42,43,62 . Another feature is the puckered out-of-plane crystal relaxations 46 . The vertical layer separation is shortest for the AB/BA stacking regions and largest for the AA regions.

To obtain the mechanical relaxation pattern of a twisted bilayer graphene, we adopted a continuum model in combination with the GSFE 43-45 . To describe the structural deformation, at an unrelaxed position r in a supercell, we define the corresponding in-plane component of the displacement vector U ( ) i ( r ) ( i = 1 2) , and the out-of-plane component h ( ) i ( r ) = h ( ) i ( r )ˆ z after relaxation with layer index i . The undeformed position r before relaxation is then mapped as r → r + U ( ) i ( r ) + h ( ) i ( r )ˆ z . The total mechanical energy of the twisted system has two components, the intra- and inter-layer components. The intra-layer strain energy of a single layer sheet is described by a linear isotropic continuum approximation 42 :

$$E ^ { \text{intra} } ( U ^ { ( i ) } ) = \sum _ { i } \int \frac { G } { 2 } ( \partial _ { x } U _ { x } ^ { ( i ) } + \partial _ { y } U _ { y } ^ { ( i ) } ) ^ { 2 }$$

$$+ \frac { K } { 2 } \left \{ ( \partial _ { x } U _ { x } ^ { ( i ) } - \partial _ { y } U _ { y } ^ { ( i ) } ) ^ { 2 } + ( \partial _ { x } U _ { y } ^ { ( i ) } + \partial _ { y } U _ { x } ^ { ( i ) } ) ^ { 2 } \right \} d ^ { 2 } r ( 1 )$$

where G and K are shear and bulk modulus of a monolayer graphene, which we take to be G = 9 0eV . / Å , 2 K = 13 2 eV . / Å . 2 These values are obtained with DFT by isotropically straining and compressing the monolayer and performing a quadratic fitting of the ground-state energy as a function of the applied strain or shear.

The inter-layer energy is described by the GSFE 43-45 , denoted as V GSFE , which has been employed to explain relaxation in van der Waals heterostructures 43,45 , and depends only on the relative stacking between two adjacent layers. We obtain the V GSFE by applying a 9 × 9 grid sampling of rigid shifts to L1 in the unit cell with respect to L2 and extract the relaxed ground state energy at each shift. The optimal inter-layer separation h ( r ) = | h (1) ( r ) -h (2) ( r ) | is also extracted for each shifted configuration. The GSFE at any position r can then be expressed as a Fourier sum:

$$V ^ { \text{GSFE} } ( r ) = \sum _ { p _ { i } } V _ { p _ { i } } e ^ { i p _ { i } \cdot r },$$

where p i is defined as p i = n 1 G 1 + n 2 G 2 for n , n 1 2 integers, and V p i is the corresponding Fourier coefficient found by fitting the ground state energy at each shift. In terms of the V GSFE , the inter-layer energy E inter can be then written as follows for a relaxed twisted bilayer:

$$\underset { \text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text	$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text
$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text#\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{	\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{#\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text$$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\test{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\test$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text}$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text`$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text'$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text<$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text1$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text.$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$
$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$$$

where b r ( ) is the local stacking order at an atomic position r , which we can take to be the distance from the given atomic position r to the position of the nearest neighbor of the same sublattice type 63 . Note that the V GSFE is a function of the sum of the local stacking order and the displacement vectors to describe the inter-layer stacking energy after relaxation.

The total energy is the sum of the inter-layer and the

intra-layer energies:

$$E ^ { t o t } ( U ^ { ( i ) } ) = E ^ { i n t r a } ( U ^ { ( i ) } ) + E ^ { i n t e r } ( U ^ { ( i ) } ), \quad ( 4 ) \quad \underset { \text{num} } { \circ } \, \stackrel { \circ } { \text{%} \, \widehat { \dots } } \, \right.$$

where U (1) ( r ) = -U (2) ( r ) due to the mirror symmetry between L1 and L2 42 (i.e. a C 2 rotation in threedimensional space with an in-plane rotation axis). We then minimize the total energy as a function of the inplane displacement field U ( ) i ( r ) to obtain the optimal relaxation pattern. The relaxation pattern respects the three-fold rotation symmetry and mirror symmetry of the twisted bilayer, and the functional form can be Fourier expanded:

$$U ^ { ( i ) } ( r ) = - i \sum _ { p _ { i } } U _ { p _ { i } } ^ { ( i ) } e ^ { i p _ { i } \cdot r }, h ^ { ( i ) } ( r ) = h _ { 0 } ^ { ( i ) } + \sum _ { p _ { i } } h _ { p _ { i } } ^ { ( i ) } e ^ { i p _ { i } \cdot r }, \\. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad.. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad.$$

where symmetry requires U ( ) i R 6 p = R 6 U ( ) i p and h ( ) i R 6 p = h ( ) i p with R 6 the 2 π/ 6 rotation matrix.

## B. Ab-initio Tight-Binding Hamiltonian

Given a relaxed tBLG crystal with a commensurate supercell structure, the conventional Bloch theorem applies due to the presence of supercell translation symmetries. Modeling such a crystal with full DFT simulations at small twist angles is computationally demanding so we instead employ the ab initio tight-binding Hamiltonian method 39,60 . In this model, atomic orbital couplings are short-ranged and determined by local geometry such as the atomic registry 60 , layer separation 64 , strain 39 and orientation in the heterostructure. Calculations of the much smaller aligned bilayer structure in DFT allow us to extract the relevant tight-binding parameters and their dependence on the local geometry 60,64 . This is based on the Wannier transformation of DFT calculations. The large twisted supercell structure can then be parametrized by these ab initio tight-binding Hamiltonians. We have validated the method with the full DFT simulation of a tBLG at larger twist angles 60 . This approach has been applied to the rigid twisted bilayers and the tBLG under external pressure 64 . The scaling with pressure is consistent with the recent experiment of the tBLG in a pressure cell . 3 The resulting accurate ab initio tight-binding Hamiltonians are still rather complicated but can be expanded to derive the simpler low-energy effective theories in the next section.

## C. Low-Energy Effective k.p Hamiltonians

To obtain insights on the electronic structure and simpler efficient computational models, we derive the continuum low-energy effective model based on the expansion of the ab-initio tight-binding Hamiltonians of a relaxed tBLG. We begin with a brief review of the symmetries of the crystal, the low-energy effective Hamiltonian for the unrelaxed twisted bilayer graphene 30 , and then the generalization to a relaxed crystal. This is followed by a numerical projection method to derive the effective models from ab initio tight-binding Hamiltonians, which is described in the following sections.

FIG. 2. Conventions for a tBLG crystal and the symmetries. (a) In the real space, the twisted bilayer is generated by the rotation around the common hexagon center from an AA-stacking bilayer unit. The symmetries are defined with respect to this hexagon centers. M is a valley-preserving mirror symmetry which is a 180 ◦ rotation with respect to an in-plane axis shown as dotted line. (b) The monolayer Brillouin zones from both layers, which are rotated from each other. The effective long-wavelength theories are expanded at either the K or K' valleys. (c) Illustration of the low-energy k · p effective model in momentum space: the filled (empty) circles are momentum states from L1 (L2). The inter-layer (intra-layer) couplings terms with momentum q i ( p i ) are represented by the red (blue) lines which connect the coupled momentum states.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00058v1-figure-812a2c3f-f195-421a-8bbd-3c2ebd7ac079.jpg)

Monolayer graphene features relativistic Dirac electrons at K, K' valleys in the Brillouin zone corners. The twist angle introduces a relative displacement of Dirac cones, which originates from the same valley of each individual layer. These two copies of the Dirac electrons are then coupled through the inter-layer interaction, which is spatially varying with the moiré pattern. The scat-

tering between Dirac electrons from the opposite valleys (K and K') is negligible due to the much larger scattering momentum required compared to the characteristic momentum scale of the moiré pattern at small twist angle. Therefore, Dirac electrons from K and K' valleys are essentially decoupled in the single particle description. Such valley protection and degeneracy, or the emergent valley U v (1) symmetry 56 , allows for the construction of effective models involving only one single K valley type with the opposite valley related by time reversal. Each copy of the effective model involves low-energy states enclosed in one circle in Fig. 2 (b) with the rotated monolayer Brillouin zones.

The above view of the effective k · p low-energy theory stresses its role as the expansion of the full model. On the other hand, symmetry considerations have been used as the guiding principle in constructing these low-energy models. The relevant symmetry operations for a tBLG effective k p · model include the lattice translations, C 3 rotations, valley-preserving M mirror symmetry and antiunitary C 2 T symmetry 56 . In our convention, the origin (rotation center) is chosen to be the hexagon center of a locally AA-stacking spot with the symmetries shown in Fig. 2 (a). Under these symmetries, the wavefunction Ψ( ) = [Ψ r 1A ( r ) , Ψ 1B ( r ) , Ψ 2A ( r ) , Ψ 2B ( r )] transforms as

$$C _ { 3 } \colon & \quad \mathcal { M } \colon & \quad \cdots \quad \cdots \quad \cdots \quad \begin{array} { c c c c } \text{setc} \\ C _ { 3 } \colon & \mathcal { M } \colon & C _ { 2 } \mathcal { T } \colon & \text{this } \iota \\ \widetilde { \phi } \Psi _ { 1 A } ( \mathcal { R } _ { 3 } r ) \\ \widetilde { \phi } \Psi _ { 1 B } ( \mathcal { R } _ { 3 } r ) \\ \widetilde { \phi } \Psi _ { 2 A } ( \mathcal { R } _ { 3 } r ) \\ \widetilde { \phi } \Psi _ { 2 B } ( \mathcal { R } _ { 3 } r ) \end{array} \right ], & \begin{bmatrix} \Psi _ { 2 B } ( \mathcal { M } r ) \\ \Psi _ { 2 A } ( \mathcal { M } r ) \\ \Psi _ { 1 B } ( \mathcal { M } r ) \\ \Psi _ { 1 A } ( \mathcal { M } r ) \end{bmatrix}, & \begin{bmatrix} \Psi _ { 1 B } ^ { * } ( - r ) \\ \Psi _ { 1 A } ^ { * } ( - r ) \\ \Psi _ { 2 B } ^ { * } ( - r ) \\ \Psi _ { 2 A } ^ { * } ( - r ) \end{bmatrix} & \text{$(6)$} \quad \text{fit} \, \iota \\ \text{foldec} \\ \text{where $\theta = e^{1}2\pi/3(\mathcal{A}=e^{-i2\pi/3\mathcal{M}}\mathcal{R}$} \text{rotate the vector $r$} \quad \text{this $r$} \end{bmatrix}$$

where φ = e i 2 π/ 3 ( ¯ = φ e -i 2 π/ 3 ), R 3 rotates the vector r clockwise by 2 π/ 3 and M flips the x coordinate of the vector r as in Fig. 2 (a). The crystal relaxation retains these relevant symmetries.

Such an effective low-energy theory expanded around the K valley has already been derived for the unrelaxed case of tBLG 30,52-54 . Since here we want to eventually include the electronic effects of relaxation, we augment this model to include additional terms that can capture these effects. With the gauge convention chosen such that the origin coincides with a locally AA-stacking spot, the augmented low-energy Hamiltonian takes the form:

$$\tilde { H } _ { \text{K} } = \begin{bmatrix} H _ { D } ^ { ( 1 ) } ( \vec { k } ) + A ^ { ( 1 ) } ( \vec { r } ) & \tilde { T } ^ { \dagger } ( \vec { r } ) \\ \tilde { T } ( \vec { r } ) & H _ { D } ^ { ( 2 ) } ( \vec { k } ) + A ^ { ( 2 ) } ( \vec { r } ), \end{bmatrix} \quad ( 7 ) & \begin{pmatrix} \text{tors.} \\ \text{cross} \alpha \\ \text{Hami} \end{pmatrix}$$

where H ( ) i D ( k ) is the Dirac Hamiltonian for each individual layer ( i = 1 2 , ), and ˜ ( T r ) the inter-layer coupling matrix, which varies with the spatial moiré pattern. The Dirac Hamiltonian is given by:

$$H _ { D } ( { \boldsymbol k } ) = v _ { F } \begin{bmatrix} 0 & - i k _ { + } \\ i k _ { - } & 0 \end{bmatrix}, \quad \quad ( 8 ) \quad \text{A.\ } \text{I\ }$$

with additional Pauli matrix rotation terms to account for the small twist angle.

In the unrelaxed case, the terms A ( ) i ( r ) , i = 1 2 , are set to zero, and the inter-layer coupling ˜ ( T r ) is simplified to:

$$T _ { \alpha \beta } ( r ) = \omega _ { 0 } \begin{bmatrix} T _ { 0 } ( r ) & T _ { + } ( r ) \\ T _ { - } ( r ) & T _ { 0 } ( r ) \end{bmatrix}, \quad \quad ( 9 )$$

where T αβ ( r ) specifies the coupling between various momentum states from the sublattice β of L1 to sublattice α of L2, α, β = A B , . In the above equations, the various symbols that appear have the following meaning: k ± = k x ± ik y , T 0 ( r ) = e i q 1 · r + e i q 2 · r + e i q 3 · r , T -( r ) = e i q 1 · r + φe i q 2 · r + ¯ φe i q 3 · r , T + ( r ) = e i q 1 · r + ¯ φe i q 2 · r + φe i q 3 · r . q 1 = k D ˆ x , q 2 = k D ( -ˆ x -√ 3ˆ) y / 2 and q 3 = k D ( -ˆ + x √ 3ˆ) y / 2 with k D = (8 π/ a 3 G ) sin( θ/ 2) and a G being the graphene lattice constant, as shown in Fig. 2 (c). The inter-layer coupling strength is ω 0 ≈ 110 meV. The moiré supercell has reciprocal lattice vectors G 1 = q 1 -q 2 and G 2 = q 3 -q 2 . The Hamiltonian can be shown to preserve the symmetries above. One interesting observation is that the coupling T ( r ) is not periodic under moiré supercell translations 65 due to the field expansion gauge choice around the Dirac point of each individual layer. The two shifted Dirac cones are connected by q 1 as in Fig. 2 (b) and (c). We will re-examine this labeling of momentum states which would facilitate the projection from the full tight-binding Hamiltonian.

To determine the set of coupled momentum states, we first look for the momentum states in each layer that are folded onto the same supercell momentum label. With k i being the momentum measured from the K ( ) i valley, this requires K (1) + k 1 = K (2) + k 2 + G , where G is a reciprocal lattice vector of the moiré supercell spanned by G i . The set of k i momentum states of each layer (filled and empty circles) that are folded to the supercell Γ point are shown in Fig. 2 (c), which is a bipartite lattice in momentum space. In the absence of the deformation from lattice relaxation, the exact translation symmetry within each single layer eliminates the direct intra-layer coupling terms between momentum states. This leaves only the inter-layer coupling T ( r ) to be determined. The inter-layer coupling T ( r ) describes the coupling between states at momentum k 2 -k 1 = (K (1) -K (2) ) -G , which does not belong to moiré supercell reciprocal lattice vectors. In a rigid tBLG, the dominant contributions in T ( r ) can be derived from the orbital coupling in the microscopic tight-binding Hamiltonian 30 . The low-energy Hamiltonian can be viewed as a momentum lattice with only nearest neighbor couplings. The 2 × 2 coupling matrices to the three nearest neighbors direction q i are:

$$\check { \ h e } ^ { 1 } \, = \, \omega _ { 0 } \begin{bmatrix} 1 & 1 \\ 1 & 1 \end{bmatrix}, \ T _ { 2 } = \omega _ { 0 } \begin{bmatrix} 1 & \bar { \phi } \\ \phi & 1 \end{bmatrix}, \ T _ { 3 } = \omega _ { 0 } \begin{bmatrix} 1 & \phi \\ \bar { \phi } & 1 \end{bmatrix}. \\ \hat { \nu } \quad \ \lambda \quad \ \text{and} \quad \, \quad \, \sin \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colots \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \col. \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \col.$$

A momentum cutoff is also imposed for the momentum lattice within the linear Dirac cone regime.

To generalize the above effective Hamiltonian for a

twisted bilayer crystal with atomic relaxation within the plane (in-plane strain) as well as height variations (out-of-plane strain), we start with a deformed monolayer graphene with a strain field U ( ) i ( r ) , which describes the displacement vector of the underlying constituent atoms, moving the atom to a new position r → r + U ( ) i ( r ) + h ( ) i ( r )ˆ z on the i -th layer. The in-plane strain is known to appear in the low-energy theory as the pseudo-gauge field with A ( ) i ( r ) a 2 × 2 matrix which couples with the Dirac Hamiltonian H ( ) i D ( k ) + A ( ) i ( r ) 38 . The isotropic part of the strain ( ∂ U x ( ) i x + ∂ U y ( ) i y ) or the deformation potential contributes to the diagonal part of A ( ) i ( r ) , while the off-diagonal elements are given by the coupling to ( ∂ U x ( ) i x -∂ U y ( ) i y ) and ( ∂ U x ( ) i y + ∂ U y ( ) i x ) strain components dictated by symmetry. The spatially varying A ( ) i ( r ) matrix introduced by the atomic deformation induces coupling between different momentum states within the same single layer. The A ( ) i ( r ) matrix can be decomposed into Fourier components at moiré supercell reciprocal lattice vectors. Strain perturbations are also known to renormalize the local Fermi velocity of the Dirac electron and give anisotropic velocity corrections 38 . Here, we only retain the contributions in the form of the A ( ) i ( r ) matrix, which causes scatterings between the momentum states. For a twisted bilayer, the mirror symmetry will relate the A ( ) i ( r ) between the two layers.

Regarding the inter-layer coupling, the deformation field U ( ) i ( r ) + h ( ) i ( r )ˆ z from relaxation causes additional relative shifts in the local atomic registry, which modifies the T ( r ) coupling matrix. The height variations in the moiré supercell weaken the inter-layer coupling in the locally AA stacking region (due to its larger vertical separation) and enhance the coupling in the enlarged AB/BA domains. This causes the off-diagonal coupling constants in ˜ ( T r ) to be larger than the diagonal components. The dominant T 1 above with atomic relaxation becomes ω σ 0 0 + ω σ 1 x with ω 0 &lt; ω 1 . These corrections are relevant for the single-particle gaps above and below the flat bands near the magic angle. Furthermore, the domain line formation and the finer structure from the atomic relaxation enhances the scattering with higher momentum components. The dependence of the coupling parameters on the twist angle will be discussed later.

In short summary, the intra-layer deformation of each of the two layers introduces the non-zero scattering terms, A ( ) i ( r ) , in Eq. (7). The Dirac electrons from the two layers are then coupled through the inter-layer ˜ ( T r ) term, which is also modified by the atomic relaxation. Later, we will include another correction term to the inter-layer coupling that involves momentum dependence of the scattered momentum states.

The effects of various coupling terms here can also be visualized from the momentum space representation of the Hamiltonian. For simplicity, we first focus on the supercell Γ point in momentum space. The relevant momentum states from both layers are represented by the circles in Fig. 2 (c). The inter-layer coupling from ˜ ( T r ) links the states from both layers, while the intra-layer A ( ) i ( r ) connects states reside within the same layer. We can expand these coupling terms into the Fourier momentum components

$$\intertext { a } _ { \substack { a } _ { i } } ^ { ^ { \prime } } & A ^ { ( i ) } ( r ) = \sum _ { p _ { i } } A _ { p _ { i } } ^ { ( i ) } e ^ { i p _ { i } \cdot r }, \ \tilde { T } ( r ) = \sum _ { q _ { i } } \tilde { T } _ { q _ { i } } e ^ { i q _ { i } \cdot r }, \quad ( 1 1 )$$

with the twelve q i ( p i ) momentum vectors illustrated in Fig. 2 (c). The p i vectors are the reciprocal lattice vectors of the moiré supercell, p i = n 1 G 1 + n 2 G 2 and the general q i vectors are (K (1) -K (2) ) -G ′ . These correction terms introduce further neighbor couplings in the momentum lattice to the model as in Fig. 2 (c).

We now discuss how the symmetries in the low-energy theory constrain the forms of the coupling terms. For the C 3 and C 2 T symmetries, which do not flip the layer indices, we let F ( r ) represent either of A ( ) i ( r ) or ˜ ( T r ) 2 × 2 matrices. These symmetries dictate

$$C _ { 3 } \colon \ e ^ { i ( 2 \pi / 3 ) \sigma _ { z } } F ( \mathcal { R } _ { 3 } ^ { - 1 } r ) e ^ { - i ( 2 \pi / 3 ) \sigma _ { z } } & = F ( r ), \quad ( 1 2 ) \\ C _ { \lambda } \mathcal { T } \cdot \quad \sigma _ { \ } F ^ { * } ( - \mathbf r ) \sigma _ { \ } = F ( \mathbf r ) \quad ( 1 3 )$$

$$\mathfrak { r } _ { \text{he} } ^ { + } \quad C _ { 2 } \mathcal { T } \colon \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

( R -1 3 rotates r counterclockwise) with the Pauli matrices σ = ( σ , σ x y , σ z ) acting on the sublattice degree of freedom. As for the M mirror symmetry that flips the layer, we can derive the following relations

$$\underset { \text{local} } { \text{aou} } \quad \mathcal { M } \colon \ \sigma _ { x } \tilde { T } ( \mathcal { M } r ) \sigma _ { x } = \tilde { T } ^ { \dagger } ( r ), \ \sigma _ { x } A ^ { ( 1 ) } ( \mathcal { M } r ) \sigma _ { x } = A ^ { ( 2 ) } ( r ), \\ \text{mod} _ { \cdot }$$

and similarly for A (1) ↔ A (2) in the last equation. When these 2 × 2 matrices are decomposed into Fourier momentum components, these symmetry conditions yield the following equivalent constraints on the expanded matrix components:

$$\begin{array} {$$

$$\mathcal { M } \colon \ \sigma _ { x } A _ { \math$$

$$\underset { \mathfrak { r }, \quad } { \mathfrak { r$$

In the original effective model 30 , only the first three q i 's (connections to the nearest neighbors) from ˜ T q i are included. Even though these three q i 's are shown to be the dominant contribution in an unrelaxed twisted bilayer, relaxation would not only modify the matrix elements but also enhance the higher momentum components.

## D. Low-Energy Expansion based on Ab-initio Tight-Binding Hamiltonian

Here, we summarize the procedure to extract these relevant matrix elements for both the pseudo-gauge field A ( ) i ( r ) and inter-layer coupling ˜ ( T r ) in the effective lowenergy theory, derived from an ab initio tight-binding

Hamiltonian of relaxed twisted bilayer graphene.

(i) For a relaxed commensurate twisted bilayer with atomic relaxation, the supercell is spanned by the translation vectors R i . Two integers M and N are used to specify the twist angle 46 . The origin is chosen to be at the center of the hexagon as in Fig. 2 (a). We adopt the series M = N + 1 commensurate supercells which have vanishing twist angles with increasing M and exactly one AA-stacking region in a supercell unit. For an unrelaxed twisted bilayer, these specify all the atomic positions r l in a supercell. The relaxation moves an atom to a new position by r l → r l + U ( ) i ( r l ) + h ( ) i ( r l )ˆ z and modifies the coupling strength between atomic pairs. Given a supercell momentum k SC , we define the Bloch wave as

$$| \Psi _ { l } ( k _ { \text{SC} } ) \rangle$$

where N SC is the total number of supercells. Note that the original unrelaxed position r l + R is used in the definition even though the actual relaxed position is r l + U ( ) i ( r l )+ h ( ) i ( r l )ˆ+ z R . The Hamiltonian H TBH ( k SC ) in the supercell reciprocal space can be derived as:

$$H _ { l m } ^ { \text{TBH} } ( k _ { \text{$$

(ii) The relevant states at low-energy are the Bloch waves near the K valley of both layers.

$$\begin{array} { c c c c } \text{wavco in} \ a v e n \ u n \ u n \ u n \ u n \ u n \colon & \ a v e n \\ | \Phi _ { \alpha } ( K ^ { ( i ) } + k ^ { ( i ) } ) \rangle & & & q _ { j } ) | H \\ = \frac { 1 } { \sqrt { N _ { \text{SC} } N _ { \alpha } } } \sum _ { R } \sum _ { r _ { l } \in \alpha } e ^ { i ( K ^ { ( i ) } + k ^ { ( i ) } ) \cdot ( r _ { l } + R ) } | r _ { l } + R \rangle & & \text{the c} \\ & \text{pling} \\ = \frac { 1 } { \sqrt { N _ { \alpha } } } \sum _ { r _ { l } \in \alpha } e ^ { i G ^ { \prime } \cdot r _ { l } } | \Psi _ { l } ( k _ { \text{SC} } ) \rangle, & ( 2 0 ) & & k \cdot p \\ & & \text{tight} \\ & & 3 ( c ). \\ \text{with } i = 1, 2 \text{ being the layer index, } \alpha = A, B \text{ the sub} \\ & & \text{the b} \end{array}$$

with i = 1 2 , being the layer index, α = A,B the sublattice index, and k ( ) i is the momentum measured rel-

## E. Momentum-Dependent Inter-layer Coupling

Experimentally, it is interesting to investigate whether the tBLG is particle-hole symmetric or not under electrical gating and the implications for the correlated insulating and superconducting phases. For the electronic structures, higher order terms in the Dirac Hamiltonian, Pauli matrix rotations, due to the twist angle 32 and electric potential from doping 66 are known to give rise to particle-hole asymmetric bands. In our model, we have included the Pauli matrix rotations but they are not adequate to fully capture the asymmetric bands in the tight-binding calculations as shows in Fig. 3 (c). The projection method introduced above from the full tight- ative to K ( ) i . These states are defined with the unrelaxed original positions and momentum labeling from graphene primitive unit cell translations. To be compatible with the translational symmetry of the above Hamiltonian H TBH ( k SC ) , the momentum states must be such that they can be folded to the supercell k SC point in momentum space. This yields the condition K ( ) i + k ( ) i = G ′ + k SC with G ′ being a supercell reciprocal lattice vector. Given a supercell momentum k SC , we can then construct a set of Bloch plane wave states for both layers to sample the low-energy sector of the Hamiltonian. The unrelaxed r l + R positions of atoms used in the Bloch phase factors ensure the orthogonality of the projection basis states, for a crystal with or without atomic relaxation. The state | Φ (K α ( ) i + k ( ) i ) 〉 in the tight-binding orbital picture is mapped to the plane wave state on i -th layer with momentum k ( ) i in the low-energy k · p expansion. We will use this picture to bridge the tight-binding and low-energy k · p pictures.

(iii) Having established the tight-binding Hamiltonian and the relevant low-energy states | Φ (K α ( ) i + k ( ) i ) 〉 , we can project out the coupling matrix elements in A ( ) i p j and ˜ T q j for the low-energy Hamiltonian. For example, the intra-layer scattering terms from A ( ) i αβ with momentum p j can be inferred from 〈 Φ (K α ( ) i + k ( ) i + p j ) | H TBH ( k SC ) Φ | β (K ( ) i + k ( ) i ) 〉 with K ( ) i + k and K ( ) i + k + p j both folded to k SC . The inter-layer coupling ˜ T αβ at momentum q j can be obtained from 〈 Φ (K α (2) + k + q j ) | H TBH ( k SC ) Φ | β (K (1) + ) k 〉 with the momentum transfer q j and compatible momentum states. Numerically, the coupling matrix is obtained from the average of sampling k SC and k ( ) i momenta near the K valley that are scattered into new states. The reconstructed low-energy k · p Hamiltonian captures most essential features of the tight-binding electronic band structure as shown in Fig. 3 (c). To further improve the particle-hole asymmetry of the bands, we identify the relevant terms to be added to the low-energy theory in the next section.

binding Hamiltonian enables us to systematically extract and identify various terms and approximations. In the above k · p low-energy Hamiltonian, the inter-layer couplings between | Φ (2) α (K (2) + k + q j ) 〉 and | Φ (1) β (K (1) + k ) 〉 states are assumed to only depend on the momentum transfer q j . However, there is no protecting symmetry for this property and the coupling constants could generally depend on the momentum k as well. This dependence is explicitly seen in our numerical projection method. To elucidate these momentum-dependent correction terms, we focus here on the first three dominant q j contribu-

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00058v1-figure-58c954da-66ed-45af-bc1b-9d509605307c.jpg)

FIG. 3. (a) A comparison of the band structures from decoupled bilayers, one without relaxation (thick green lines) and one with relaxation (purple lines). The in-plane and out-of-plane relaxation can open band gaps at high-symmetry points. (b-d) Band structures benchmarked for effective k p · bilayer Hamiltonians at θ = 1 05 . ◦ with a momentum lattice cutoff radius of 6 | G | . The full tight-binding reference bands are in thick grey lines and the k · p bands are in red. We show examples for both the unrelaxed bandstructure (b) and the relaxed bandstructure, (c) without momentum-dependent and (d) with momentum-dependent ˜ ( T ′ r , k ˆ ± ) terms in Eq. (21).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00058v1-figure-5db011bb-e099-4df1-a811-d50c8f80f336.jpg)

tions, which gives the following form

$$\tilde { T } ^ { \prime } ( r, \hat { k } _ { \pm } ) = & \frac { i } { 2 } \begin{bmatrix} \lambda _ { 1 } \{ T _ { - } ( r ), \hat { k } _ { + } \} & \lambda _ { 3 } \{ T _ { 0 } ( r ), \hat { k } _ { + } \} \\ \lambda _ { 2 } \{ T _ { + } ( r ), \hat { k } _ { + } \} & \lambda _ { 1 } \{ T _ { - } ( r ), \hat { k } _ { + } \} \end{bmatrix} \quad \text{inter} \\ & - \frac { i } { 2 } \begin{bmatrix} \lambda _ { 1 } \{ T _ { + } ( r ), \hat { k } _ { - } \} & \lambda _ { 2 } \{ T _ { - } ( r ), \hat { k } _ { - } \} \\ \lambda _ { 3 } \{ T _ { 0 } ( r ), \hat { k } _ { - } \} & \lambda _ { 1 } \{ T _ { + } ( r ), \hat { k } _ { - } \} \end{bmatrix}, \quad \begin{array} { c } \text{inter} \\ \text{super} \\ \text{and} \ \text{cases} \\ \text{ences} \end{array}$$

where λ 2 ≈ 2 λ 1 ≈ 0 18 . eV Å, · λ 3 ≈ 0 around θ = 1 ◦ and the anti-commutator form for non-commuting operators. We have added this leading momentum-dependent scattering terms in the inter-layer coupling, ˜ ( T r )+ ˜ ( T ′ r , k ˆ ± ) , and they indeed can capture the particle-hole asymmetry as shown in Fig. 3 (d). The projected low-energy bands are in good agreement with the full ab initio tight-binding electronic structure.

## III. NUMERICAL RESULTS FOR ANGULAR DEPENDENCE

The effective models presented above are derived for twist angles in the range of 0 7 . ◦ to 4 ◦ . A series of commensurate ( M,N ) supercells with M = N +1 are used to extract the model parameters in this range. For the general twist angles, the parameters can be obtained by interpolating the results of the commensurate cases. The numerical results are presented in their Fourier component version. There are two sets of crystal momenta involved in the expansions. The first is denoted as p i = n 1 G 1 + n 2 G 2 which can be expressed by the moiré supercell reciprocal lattice vectors G i . The other second can be written as q i as in Fig. 2. The differences between q i vectors are reciprocal lattice vectors of the moiré supercell. Hence, any q i can be written as q i = q 1 + m i 1 G 1 + m i 2 G 2 . The q i are the basis for the inter-layer coupling expansion, while the p i are for the intra-layer expansion (strain). Terms evaluated from momenta of identical magnitude tend to only differ by a complex phase dictated by symmetries. For this reason, we organize the p and q dependent terms by 'shells' of equidistant momenta. For example, in Fig. 2, there are three shells of q (in red) and two shells of p (in blue). In summarizing the strength of terms in the expansion, it is more convenient to compare the average magnitudes of different shells instead of the complex values of individual momenta, so we label these terms as T s j and A s j instead of T q i and A ( ) i p i

We begin with the mechanical properties of the relaxed twisted bilayer graphene. The competition between the stacking energy and the strain energy will depend on the

twist angle which determines the moiré length scale. The atomic structure will tend to reduce to the area of the unfavorable AA stacking, while increasing the area of the favorable Bernal stacking ( AB and BA ). The displacement field in Eq. (5) can be expanded into the Fourier components at momenta p i . A real-space image of the atomic relaxation is provided in Fig. 4, along with the angledependent magnitudes of the representative p i from the first three shells.

FIG. 4. (Left) Relaxation patterns of L2 for bilayers with twist angles θ = 1 2 . ◦ and 0 7 . ◦ from the continuum model. The local relaxation is given by the white arrows, and the curl of the relaxation (local change in twisting angle) is given by color. (Right) Twist-angle dependence of the relaxation coefficients u p i , where the label [ n , n 1 2 ] corresponds to p i = n 1 G 1 + n 2 G 2 .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00058v1-figure-0de33c59-a2e7-4b2c-baeb-dfc7dcd4214d.jpg)

The atomic deformation induces a pseudo-gauge field coupled to the Dirac electrons in each layer. The pseudogauge can also be Fourier expanded at p i momenta, and then easily included in the k · p Hamiltonian. The crystal relaxation also modifies the inter-layer coupling. Both the ˜ ( T r ) and the momentum-dependent corrections ˜ ( T ′ r ) are affected. The diagonal inter-layer coupling term (labeled ˜ T ( AA ) ) is generally smaller than the offdiagonal coupling (labeled ˜ T ( AB ) ), as the AA stacking has a larger inter-layer distance (and thus smaller electronic coupling) than the AB and BA regions. The AA coupling is further reduced at small angles by the reduction of the overall AA stacking area, while the AB coupling is increased. The dependence of all electronic terms on twist angle for the relaxed tBLG system is given in Fig. 5.

Using the fitted k · p parameters, band structures can be calculated for tBLG systems at an arbitrary twist angle (Fig. 6). These band structures are essentially identical to those of the full tight-binding model (see Fig. 3d). Our model predicts an inversion of the two low-energy bands at 0 98 . ◦ , but has nearly flat bands in a range of angles from 0 95 . ◦ to 1 05 . ◦ . For further analysis and discussion of the low-energy electronic structure and its angle-dependence we defer to a companion work, Ref. 57.

FIG. 5. Twist-angle dependence of various parameters in the k · p model. (a) Magnitudes of inter-layer coupling terms ˜ T i for the three nearest shells ( s j ) of couplings (red lines in Fig. 2 (c).). The off-diagonal elements of the ˜ T matrices (AB) are shown by solid lines, while the diagonal elements (AA) by dashed lines. (b) The three λ i momentum-dependent parameters. (c) Magnitudes of the intra-layer coupling terms, A ( ) i for the two nearest shells ( s j ) of couplings (blue lines in Fig. 2(c).).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00058v1-figure-b6c768fd-4776-480f-8331-2f2955ad8944.jpg)

## IV. CONCLUSION

In this work, we adopted a systematic multi-scale approach to the modeling of tBLG with an angle dependence and relaxation effects included. Starting with DFT calculations, the mechanical and electronic properties are extracted to model the relaxed geometry through ab-initio tight-binding Hamiltonians. We then simplify the models by low-energy k · p expansions to capture the essential electronic features of the derived ab initio tight-binding Hamiltonians 57 . Our approach provides an unified and coherent method to derive and simplify the modeling process while preserving the important physics, which paves the way for formulating interacting-electron theories (The model scripts are available at https://github.com/stcarr/kp\_tblg ).

There are also other factors that may affect the electronic structure of tBLG that are not included in our current modeling. Such factors are the hBN substrate effects on tBLG crystal relaxation and electronic cou-

FIG. 6. Twist-angle dependence of the low energy band structure of relaxed tBLG, calculated using the relaxed k p · model. The black lines are for an expansion around the monolayer K valley, while the red lines are for the monolayer K' valley. At θ = 0 98 . ◦ , the upper and lower flat bands at the Γ point pass smoothly through one another due to opposite mirror symmetry eigenvalues 57 .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00058v1-figure-d928c0c2-a1d1-4e5f-9c63-9fa4b3301e47.jpg)

plings, screening effects from the substrate used 67 , electrical gating, and self-consistent electric potential from the electron doping 66 . Disorder and additional strain variations 68 across the samples introduced by the fabrication process might also complicate the interpretation of the observed behavior . 5 Our multi-scale method provides a framework for generalizing the models to include corrections or to estimate the energy scale for such perturbations.

Our multi-scale numerical framework can be easily extended to other van der Waals heterostructures such as a few-layer graphene stacks 58,69 , and transition metal dichalcogenides stacks 59,70 . The effect of mechanical relaxation can also be included in modeling the band structure of twisted trilayer graphene stacks 71 . Following the scheme of Fig. 1, one can derive models at various levels of approximation and scaling of the parameters in conjunction with the relevant symmetry analysis. Further experimental probes to the electronic and mechanical properties can also be used to refine the numerical models. This will close the design loop in predicting the properties of van der Waals heterostructures, towards establishing these systems as a new platform for exploring the physics of correlated electrons.

## ACKNOWLEDGMENTS

We thank Yuan Cao, Alex Kruchkov, Jong Yeon Lee, Hoi Chun Po, Grigory Tarnopolsky, Yao Wang, Pablo Jarillo-Herrero, and Ashvin Vishwanath for fruitful discussions. This work was supported by the STC Center for Integrated Quantum Materials, NSF Grant No. DMR1231319 and by ARO MURI Award W911NF-14-0247. The computations in this paper were run on the Odyssey cluster supported by the FAS Division of Science, Research Computing Group at Harvard University.

## Appendix A DENSITY FUNCTIONAL THEORY CALCULATIONS

The DFT calculations in this work were carried out using the Vienna Ab initio Simulation Package (VASP) 72,73 with Projector Augmented-Wave (PAW) type of pseudopotentials, parametrized by Perdew, Burke and Ernzerhof (PBE) 74 . A slab geometry with a 20 Å vacuum region is used to reduce the interactions between periodic images. The DFT calculations are converged with plane-wave energy cutoff 500 eV and a reciprocal space Monkhorst-Pack grid sampling of size 17 × 17 × 1 .

For the mechanical relaxation, the Generalized Stacking Fault Energy (GSFE) was calculated by performing rigid shift of the unit cell and calculate the ground state energy from DFT, and we used a 9 × 9 grid in a unit cell for the computation. We fix the in-plane positions and allow the out-of-plane positions to relax. We implemented the van der Waals force through the vdW-DFT method using the SCAN+rVV10 functional 75,76 .

The extended Bloch wavefunction basis can be transformed into the maximally-localized Wannier functions (MLWF) basis as implemented in the Wannier90 code 77 . With this transformation, the effective tight-binding Hamiltonian for a designated group of bands of the material can be constructed. This not only gives an efficient numerical method to reproduce DFT results but also provides a physically transparent picture of localized atomic orbitals and their hybridizations. Our work is based on the systematic analysis of such tight-binding Hamiltonians, which inherit the ab initio information without fitting procedures for the numerical parameters 60 . Further corrections for band gaps from advanced GW calculations or other choices of exchange correlation functionals are also compatible with Wannier constructions.

- 1 Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, R. C. Ashoori, and P. JarilloHerrero, Nature 556 , 80 EP (2018).
- 2 Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556 , 43 EP (2018).
- 3 M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean, Science (2019), 10.1126/science.aav1910,
- 19 H. Isobe, N. F. Q. Yuan, and L. Fu, Phys. Rev. X 8 , 041041 (2018).
- 20 C.-C. Liu, L.-D. Zhang, W.-Q. Chen, and F. Yang, Phys. Rev. Lett. 121 , 217001 (2018).
- 21 M. Ochi, M. Koshino, and K. Kuroki, Phys. Rev. B 98 , 081102 (2018).
- 22 H. Guo, X. Zhu, S. Feng, and R. T. Scalettar, Phys. Rev. B 97 , 235453 (2018).
- http://science.sciencemag.org/content/early/2019/01/25/science.aav1910.full.pdf. ence of gate-tunable mott insulator in tri23 G. Chen, L. Jiang, S. Wu, B. Lv, H. Li, B. L. Chittari, K. Watanabe, T. Taniguchi, Z. Shi, Y. Zhang, and F. Wang, 'Evid layer graphene-boron nitride moirÃľ superlattice,' (2018), arXiv:1803.01985.
- 4 Y. Cao, D. Chowdhury, D. Rodan-Legrain, O. RubiesBigordÃă, K. Watanabe, T. Taniguchi, T. Senthil, and P. Jarillo-Herrero, 'Strange metal in magic-angle graphene with near planckian dissipation,' (2019), arXiv:1901.03710.
- 5 A. Kerelsky, L. McGilly, D. M. Kennes, L. Xian, M. Yankowitz, S. Chen, K. Watanabe, T. Taniguchi, J. Hone, C. Dean, A. Rubio, and A. N. Pasupathy, 'Magic angle spectroscopy,' (2018), arXiv:1812.08776.
- 24 G. Chen, A. L. Sharpe, P. Gallagher, I. T. Rosen, E. Fox, L. Jiang, B. Lyu, H. Li, K. Watanabe, T. Taniguchi, J. Jung, Z. Shi, D. Goldhaber-Gordon, Y. Zhang, and F. Wang, 'Signatures of gate-tunable superconductivity in trilayer graphene/boron nitride moirÃľ superlattice,' (2019), arXiv:1901.04621.
- 6 Y. Choi, J. Kemmer, Y. Peng, A. Thomson, H. Arora, R. Polski, Y. Zhang, H. Ren, J. Alicea, G. Refael, F. von Oppen, K. Watanabe, T. Taniguchi, and S. NadjPerge, 'Imaging electronic correlations in twisted bilayer graphene near the magic angle,' (2019), arXiv:1901.02997.
- 7 A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney, K. Watanabe, T. Taniguchi, M. A. Kastner, and D. Goldhaber-Gordon, 'Emergent ferromagnetism near three-quarters filling in twisted bilayer graphene,' (2019), arXiv:1901.03520.
- 8 X. Lu, P. Stepanov, W. Yang, M. Xie, M. A. Aamir, I. Das, C. Urgell, K. Watanabe, T. Taniguchi, G. Zhang, A. Bachtold, A. H. MacDonald, and D. K. Efetov, 'Superconductors, orbital magnets, and correlated states in magic angle bilayer graphene,' (2019), arXiv:1903.06513.
- 9 E. Codecido, Q. Wang, R. Koester, S. Che, H. Tian, R. Lv, S. Tran, K. Watanabe, T. Taniguchi, F. Zhang, M. Bockrath, and C. N. Lau, 'Correlated insulating and superconducting states in twisted bilayer graphene below the magic angle,' (2019), arXiv:1902.05151.
- 10 B. Lian, Z. Wang, and B. A. Bernevig, 'Twisted bilayer graphene: A phonon driven superconductor,' (2018), arXiv:1807.04382.
- 11 F. Wu, A. H. MacDonald, and I. Martin, Phys. Rev. Lett. 121 , 257001 (2018).
- 12 B. Padhi, C. Setty, and P. W. Phillips, Nano Letters 18 , 6175 (2018), pMID: 30185049, https://doi.org/10.1021/acs.nanolett.8b02033.
- 13 M. Xie and A. H. MacDonald, 'On the nature of the correlated insulator states in twisted bilayer graphene,' (2018), arXiv:1812.04213.
- 14 A. Thomson, S. Chatterjee, S. Sachdev, and M. S. Scheurer, Phys. Rev. B 98 , 075109 (2018).
- 15 J. W. F. Venderbos and R. M. Fernandes, Phys. Rev. B 98 , 245103 (2018).
- 16 J. González and T. Stauber, Phys. Rev. Lett. 122 , 026801 (2019).
- 17 D. M. Kennes, J. Lischner, and C. Karrasch, Phys. Rev. B 98 , 241407 (2018).
- 18 T. J. Peltonen, R. Ojajärvi, and T. T. Heikkilä, Phys. Rev. B 98 , 220504 (2018).
- 25 C. Shen, N. Li, S. Wang, Y. Zhao, J. Tang, J. Liu, J. Tian, Y. Chu, K. Watanabe, T. Taniguchi, R. Yang, Z. Y. Meng, D. Shi, and G. Zhang, 'Observation of superconductivity with tc onset at 12k in electrically tunable twisted double bilayer graphene,' (2019), arXiv:1903.06952.
- 26 X. Liu, Z. Hao, E. Khalaf, J. Y. Lee, K. Watanabe, T. Taniguchi, A. Vishwanath, and P. Kim, 'Spin-polarized correlated insulator and superconductor in twisted double bilayer graphene,' (2019), arXiv:1903.08130.
- 27 Y. Cao, D. Rodan-Legrain, O. Rubies-Bigorda, J. M. Park, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, 'Electric field tunable correlated states and magnetic phase transitions in twisted bilayer-bilayer graphene,' (2019), arXiv:1903.08596.
- 28 E. Maxwell, Phys. Rev. 78 , 477 (1950).
- 29 S. Carr, D. Massatt, S. Fang, P. Cazeaux, M. Luskin, and E. Kaxiras, Phys. Rev. B 95 , 075420 (2017).
- 30 R. Bistritzer and A. H. MacDonald, Proceedings of the National Academy of Sciences 108 , 12233 (2011), https://www.pnas.org/content/108/30/12233.full.pdf.
- 31 H. C. Po, L. Zou, T. Senthil, and A. Vishwanath, 'Faithful tight-binding models and fragile topology of magic-angle bilayer graphene,' (2018), arXiv:1808.02482.
- 32 Z. Song, Z. Wang, W. Shi, G. Li, C. Fang, and B. A. Bernevig, 'All "magic angles" are "stable" topological,' (2018), arXiv:1807.10676.
- 33 K. Hejazi, C. Liu, H. Shapourian, X. Chen, and L. Balents, Phys. Rev. B 99 , 035111 (2019).
- 34 H. C. Po, L. Zou, A. Vishwanath, and T. Senthil, Phys. Rev. X 8 , 031089 (2018).
- 35 M. Koshino, N. F. Q. Yuan, T. Koretsune, M. Ochi, K. Kuroki, and L. Fu, Phys. Rev. X 8 , 031087 (2018).
- 36 J. Kang and O. Vafek, Phys. Rev. X 8 , 031088 (2018).
- 37 N. F. Q. Yuan and L. Fu, Phys. Rev. B 98 , 045103 (2018).
- 38 J. L. Mañes, F. de Juan, M. Sturla, and M. A. H. Vozmediano, Phys. Rev. B 88 , 155405 (2013).
- 39 S. Fang, S. Carr, M. A. Cazalilla, and E. Kaxiras, Phys. Rev. B 98 , 075106 (2018).
- 40 M. A. H. Vozmediano, F. de Juan, and A. Cortijo, Journal of Physics: Conference Series 129 , 012001 (2008).
- 41 M. Vozmediano, M. Katsnelson, and F. Guinea, Physics Reports 496 , 109 (2010).

- 42 N. N. T. Nam and M. Koshino, Phys. Rev. B 96 , 075311 (2017).
- 43 S. Carr, D. Massatt, S. B. Torrisi, P. Cazeaux, M. Luskin, and E. Kaxiras, Phys. Rev. B 98 , 224102 (2018).
- 44 S. Zhou, J. Han, S. Dai, J. Sun, and D. J. Srolovitz, Phys. Rev. B 92 , 155438 (2015).
- 45 S. Dai, Y. Xiang, and D. J. Srolovitz, Nano Letters 16 , 5923 (2016), pMID: 27533089, https://doi.org/10.1021/acs.nanolett.6b02870.
- 46 K. Uchida, S. Furuya, J.-I. Iwata, and A. Oshiyama, Phys. Rev. B 90 , 155451 (2014).
- 47 S. Latil, V. Meunier, and L. Henrard, Phys. Rev. B 76 , 201402 (2007).
- 48 G. Trambly de LaissardiÃČÅąre, D. Mayou, and L. Magaud, Nano Letters 10 , 804 (2010), pMID: 20121163, https://doi.org/10.1021/nl902948m.
- 49 E. Suárez Morell, J. D. Correa, P. Vargas, M. Pacheco, and Z. Barticevic, Phys. Rev. B 82 , 121407 (2010).
- 50 G. Trambly de Laissardière, D. Mayou, and L. Magaud, Phys. Rev. B 86 , 125413 (2012).
- 51 A. O. Sboychakov, A. L. Rakhmanov, A. V. Rozhkov, and F. Nori, Phys. Rev. B 92 , 075402 (2015).
- 52 J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Phys. Rev. Lett. 99 , 256802 (2007).
- 53 J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Phys. Rev. B 86 , 155449 (2012).
- 54 E. J. Mele, Phys. Rev. B 84 , 235439 (2011).
- 55 E. J. Mele, Phys. Rev. B 81 , 161405 (2010).
- 56 L. Zou, H. C. Po, A. Vishwanath, and T. Senthil, Phys. Rev. B 98 , 085435 (2018).
- 57 S. Carr, S. Fang, Z. Zhu, and E. Kaxiras, 'An exact continuum model for low-energy electronic states of twisted bilayer graphene,' (2019), arXiv:1901.03420.
- 58 Y.-H. Zhang, D. Mao, Y. Cao, P. Jarillo-Herrero, and T. Senthil, Phys. Rev. B 99 , 075127 (2019).
- 59 F. Wu, T. Lovorn, E. Tutuc, I. Martin, and A. H. MacDonald, Phys. Rev. Lett. 122 , 086402 (2019).
- 60 S. Fang and E. Kaxiras, Phys. Rev. B 93 , 235153 (2016).
- 61 S. Carr, S. Fang, H. C. Po, A. Vishwanath, and E. Kaxiras, 'Derivation of wannier orbitals and minimal-basis tightbinding hamiltonians for twisted bilayer graphene: a first-
- principles approach,' (2019), arXiv:1907.06282.
- 62 K. Zhang and E. B. Tadmor, Journal of the Mechanics and Physics of Solids 112 , 225 (2018).
- 63 P. C. amd Mitchell Luskin and D. Massatt, 'Energy minimization of 2d incommensurate heterostructures,' (2018), arXiv:1806.10395.
- 64 S. Carr, S. Fang, P. Jarillo-Herrero, and E. Kaxiras, Phys. Rev. B 98 , 085144 (2018).
- 65 R. Bistritzer and A. H. MacDonald, Phys. Rev. B 84 , 035440 (2011).
- 66 F. Guinea and N. R. Walet, Proceedings of the National Academy of Sciences 115 , 13174 (2018), https://www.pnas.org/content/115/52/13174.full.pdf.
- 67 D. Y. Qiu, F. H. da Jornada, and S. G. Louie, Nano Letters 17 , 4706 (2017), pMID: 28677398, https://doi.org/10.1021/acs.nanolett.7b01365.
- 68 Z. Bi, N. F. Q. Yuan, and L. Fu, 'Designing flat band by strain,' (2019), arXiv:1902.10146.
- 69 J. Y. Lee, E. Khalaf, S. Liu, X. Liu, Z. Hao, P. Kim, and A. Vishwanath, 'Theory of correlated insulating behaviour and spin-triplet superconductivity in twisted double bilayer graphene,' (2019), arXiv:1903.08685.
- 70 F. Wu, T. Lovorn, E. Tutuc, and A. H. MacDonald, Phys. Rev. Lett. 121 , 026402 (2018).
- 71 Z. Zhu, P. Cazeaux, M. Luskin, and E. Kaxiras, 'Modeling mechanical relaxation in incommensurate two-dimensional trilayers,' (2019), manuscript in preparation.
- 72 G. Kresse and J. Furthmüller, Phys. Rev. B 54 , 11169 (1996).
- 73 G. Kresse and J. FurthmÃČÅŠller, Computational Materials Science 6 , 15 (1996).
- 74 J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77 , 3865 (1996).
- 75 J. Klimeš, D. R. Bowler, and A. Michaelides, Physical Review B 83 , 195131 (2011).
- 76 H. Peng, Z.-H. Yang, J. P. Perdew, and J. Sun, Physical Review X 6 , 041005 (2016).
- 77 A. A. Mostofi, J. R. Yates, Y.-S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, Computer Physics Communications 178 , 685 (2008).
Shiang Fang Stephen Carr Ziyan Zhu Daniel Massatt Efthimios Kaxiras
Angle-Dependent {\it Ab initio} Low-Energy Hamiltonians for a Relaxed Twisted Bilayer Graphene Heterostructure
We present efficient angle-dependent low-energy Hamiltonians to describe the
properties of the twisted bilayer graphene (tBLG) heterostructure, based on
{\it ab initio} calculations of mechanical relxation and electronic structure.
The angle-dependent relaxed atomic geometry is determined by continuum
elasticity theory, which induces both in-plane and out-of-plane deformations in
the stacked graphene layers. The electronic properties corresponding to the
deformed geometry are derived from a Wannier transformation to local
interactions obtained from Density Functional Theory calculations. With these
{\it ab initio} tight-binding Hamiltonians of the relaxed heterostructure, the
low-energy effective theories are derived from the projections near Dirac cones
at K valleys. For twist angles ranging from 0.7$^\circ$ to 4$^\circ$, we
extract both the intra-layer pseudo-gauge fields and the inter-layer coupling
terms in the low-energy Hamiltonians, which extend the conventional low-energy
continuum models. We further include the momentum dependent inter-layer
scattering terms which give rise to the particle-hole asymmetric features of
the electronic structure. Our model Hamiltonians can serve as a starting point
for formulating physically meaningful, accurate interacting electron theories.