Derivation of Wannier orbitals and minimal-basis tight-binding Hamiltonians for twisted bilayer graphene: First-principles approach

Twisted bilayer graphene (TBLG) has emerged as an important platform for studying correlated phenomena, including unconventional superconductivity, in two-dimensional systems. The complexity of the atomic-scale structures in TBLG has made even the study of single-particle physics at low energies around the Fermi level, quite challenging. Our goal here is to provide a convenient and physically motivated picture of single-particle physics in TBLG using reduced models with the smallest possible number of localized orbitals. The reduced models exactly reproduce the low-energy bands of ab initio tight-binding models, including the effects of atomic relaxations. Furthermore, we obtain for the first time the corresponding Wannier orbitals that incorporate all symmetries of TBLG, which are also calculated as a function of angle, a requisite first step towards incorporating electron interaction effects. We construct eight-band and five-band models for the low-energy states for twist angles between 1.3◦ and 0.6◦. The models are created using a multistep Wannier projection technique starting with appropriate ab initio k · p continuum Hamiltonians. Our procedure can also readily capture the perturbative effects of substrates and external displacement fields while offering a significant reduction in complexity for studying electron-electron correlation phenomena in realistic situations.


I. INTRODUCTION
Recent experimental advances in nanotechnology fabrication methods [1] have made it possible to stack two or more layers of two-dimensional (2D) materials with exquisite control of the angle of twist (denoted here as θ ) in the relative orientation of the layers. This level of control in relative orientation, demonstrated in detail in twisted bilayer graphene (TBLG) but also applicable to other 2D layered systems, has opened new possibilities in controlling their electronic, transport, optical and other properties [2][3][4][5]-in essence introducing the new field of "twistronics" [6]. One of the most well known of twist-angle dependent phenomena is the emergence of correlated insulators and unconventional superconductivity in TBLG at θ ≈ 1 • , with a superconducting T c of a few degrees Kelvin [7][8][9]. Theoretical studies of TBLG near the critical (often referred to as "magic") angle of ∼1 • predicted uncommonly flat bands near the Fermi energy, providing some tantalizing hints for why strongly correlated phases can arise near that angle. This single-particle band structure for TBLG has been modeled successful by density functional theory [10,11], tight-binding Hamiltonians [6,[12][13][14], and continuum models [15][16][17][18]. All such models feature a large number of degrees of freedom, which render the study of electrons' correlation unfeasible. Even the simplest, a continuum k · p model, requires a basis of hundreds of Bloch states, which is numerically intractable for proper many-body calculations. A number of more simplified N-band models have been proposed, where N is typically in the range 2-10, not taking into account the spin and k-space "valley" degrees of freedom, the latter referring to the point in reciprocal space where the band extrema near the Fermi level appear [19][20][21][22][23][24][25]. These models are generally empirical, and chosen to address the flat bands near the magicangle while taking advantage of the roughly 30-meV band gaps that separate the flat bands and the rest of the energy bands with the inclusion of atomic relaxation [13,26]. This assumes that the correlations responsible for the observed superconductivity can be accurately described with just the flat bands. If the interaction strength in magic angle TBLG is larger than the gaps that separate the flat band manifold, a model which includes additional nearby bands is required.
Here we combine the accuracy of existing ab initio models with the simplicity of a reduced model by directly constructing Wannier functions for a large range of twist angles around θ = 1 • . Our starting point is a microscopic k · p model for TBLG which includes atomic relaxations [27,28], from which we derive three single-valley tight-binding models. The first, and more comprehensive, is an eight-band model, which includes the two flat bands and three bands above and three below them. Two additional reduced models include five bands, and involve the two flat bands and the nearest three bands either above or below them. The two five-band models are only reliable when the flat band manifold is isolated in the band structure. Near 0.85 • , the gaps above and below the flat bands vanish, and only the eight-band model continues to be trustworthy. We find the eight-band model can be reliably constructed down to θ = 0.6 • , with only technical limitations of the projection technique restricting further reduction in twist angle.
The model construction is performed with a multistep Wannier projection technique, which utilizes the chiralsymmetric k · p model [29] to define a proper subspace of the fully relaxed model. In this context, the chiral-symmetric model is not only an interesting theoretical limiting case, but also an extremely useful numerical tool. The reduced models we obtain can serve as a springboard for studying correlated phenomena but are a useful basis for performing other calculations that depend on accurate electronic structure. We also study the response of these models under two common symmetry-breaking perturbations to the TBLG system: sublattice symmetry breaking due to the presence of a hexagonal boron nitride substrate and layer symmetry breaking due to an external electric displacement field. Compared to previous efforts to produce a set of localized orbitals for the low-energy states of TBLG, the models presented here explicitly preserve all symmetries of the full tight binding models. Furthermore, they are based on realistic ab initio tight binding band structures and avoid any empirical choices of coefficients or artificial band fitting procedures.
The inclusion of additional bands in our models provides clear advantages over existing treatments that focus on only the flat bands. First, all bands within 150 meV of the charge neutrality point are retained near the magic angle, allowing for the treatment of an interaction strength which is comparable to the superlattice gap separating the nearly flat bands from the higher bands. Secondly, all symmetries and band topology are manifest which can allow, for example, one to accurately describe how Chern bands emerge from the staggered chemical potential of a hexagonal boron nitride substrate. Third, the inclusion of additional bands allows for the systematic study of angle dependence beyond the small range where the flat band manifold separates from all other bands.
Because of valley and spin symmetry, dealing with a 20 band (5 × 4) or 32 band (8 × 4) model is challenging within any exact multiparticle numerical method. However, simplifications could arise from the fact that most of the weight in the flat bands are from the triangular lattice sites within our final Wannier basis. Since the additional bands could be well approximated as being completely filled or empty, it might be possible to supplement standard numerical methods with such constraints to reduce the computational effort required. In any case, a useful middle ground must be found between completely ignoring the additional bands (which would give the wrong topology) and including them fully within a numerical calculation (which is prohibitively expensive). This challenging step is left to future work.

A. Wannierization difficulties in TBLG
Electronic bands with weak energy dispersions do not always imply trivial Bloch wave function texture throughout the Brillouin zone (BZ). For example, the flat Landau levels formed in the presence of magnetic field feature nontrivial windings which lead to the topological Chern number [30,31]. Such topological character is known to obstruct the construction of real space Wannier states [32]. On the other hand, internal and crystalline symmetries can further constrain the Hilbert space structure for a group of bands, which can have consequences in obtaining Wannier orbitals in the real space. For example, in the Z 2 topological insulator, the group of nontrivial bands can only be described by a set of exponentially localized Wannier orbitals that break the time-reversal symmetry [33].
In TBLG near the magic angle, it is relevant to construct the minimal tight-binding model as the basis for interacting theories. The natural choice is to construct a model for just the group of flat bands. For a single-valley model without spin polarization, these are the two bands closest to the Fermi level, and throughout we will refer to this two-band manifold as the "flat bands" even when considering twist angles far from the magic angle. It has been shown that when all the (emergent) spatial symmetries are taken into account (with point group D 6 ) [21], one cannot have a real space two-band representation of the flat bands [21,22]. This obstruction can be understood from the net chirality of two Dirac cones or from the symmetry representations of the flat bands [21,22]. One way to circumvent the obstruction is to give up on keeping all symmetries of the system manifest, such as by forgoing the U(1) v valley symmetry [22] or the twofold rotation symmetry [23]. These non-manifest symmetries are then realized in a nontrivial manner. The Wannier functions derived by reducing the symmetry form a honeycomb lattice and have a propellerlike (or "fidget spinner") shape, created by the superposition of three charge pockets at neighboring AAstacking spots [20,22,24,25]. These models tend to have long hopping range and more complicated interaction terms due to the peculiar shape of their Wannier orbitals. Another way to resolve the obstruction is to augment the flat band manifold with additional degrees of freedom which are artificial as in Refs. [21,34] which doubled the bands. This new set of topologically nontrivial bands cancels the topology of the original flat bands. However, these new degrees of freedom are not motivated physically from the realistic electronic structure and can complicate the modeling for interactions.
The other way to avoid the obstruction while retaining the symmetry is to include an additional group of bands from realistic electronic band structure to render the overall Hilbert space trivial in topological character. In Ref. [19], it has been pointed out that the root for the obstruction in TBLG's flat bands is a "fragile topology" [35,36]. This topological obstruction is relatively weak (fragile) when compared to the more robust Chern index, in the sense that the nontriviality can be canceled by including additional atomic insulator bands. The resulting Hilbert space can then be represented by exponentially localized Wannier orbitals. This fragile topology resolution is not unique, and different solutions have been discussed which depend on the choice of additional bands. Requiring that the expanded space captures only nearby bands imposes a physical constraint that can be used to narrow this choice. By paying the cost of additional tight-binding orbitals, these resolutions generally yield Wannier functions which 033072-2 are more conventional and short range as compared to the propeller solutions.

B. Maximally localized Wannier functions
A common approach to generate a localized model is to construct maximally localized Wannier functions (MLWF) [37,38]. They define tight-binding wave functions, φ n (r), which minimize the total wave-function spread The simplest construction of MLWF for N bands relies on the calculation of M N Bloch states ψ mk over a uniform sampling of the BZ, and then evaluation of the overlap matrix elements of the periodic parts of the Bloch states [ where b is the momentum connecting each k to its nearest neighbor on the grid sampling of the BZ. The overlap elements are then used to calculate a set of k-dependent unitary rotations U (k) mn that gives the proper mixing of the original Bloch states over the BZ to minimize the real space spread [37]. Then the nth Wannier function at unitcell R is given by Applying this methodology to a k · p model is not entirely trivial however. If one attempts to construct a sampling of the k points of the BZ, one will find that the basis elements of the low-energy expansion are not smoothly periodic across the BZ boundary. That is to say: the nth basis element of the continuum Hamiltonian at k will not correspond to the n th basis element at k + δ if δ crosses a BZ boundary. This is an unavoidable consequence of the non-periodic truncation of the low-energy continuum model. The problem can be alleviated by selecting only a subset of bands, and manually matching basis elements that have nearly identical momentum. This allows for an approximate evaluation of of M k,b mn even when k + b crosses a BZ boundary.
Even with these choices, the eight-band and five-band electronic structures of TBLG are still not easily described with this approach. An unconstrained optimization of the MLWF (using, for example, the package WANNIER90 [39]) results in Wannier functions that have no clear symmetry within the moire supercell, and the required band structure symmetries are impossible to impose as one truncates the range of the resulting tight-binding model. Selectively localized Wannier functions (SLWF) [40] can minimize the spread while constraining the Wannier function centers, but still cause similar symmetry issues without reducing the spread compared to our projection method. There is a symmetry-constrained WF construction technique [41] but currently it cannot be used alongside band disentanglement, a necessary step to properly construct the bands farthest from the Fermi energy.
The standard MLWF techniques perform poorly because while the wave functions of the flat band manifold are well localized in real space, the states just above and below the manifold are not. When constructing a set of Wannier functions for both of these types of bands, the total spread functional can be minimized by making all the Wannier functions of similar spread. Doing so hybridizes the states and makes the interpretation of the WFs as a tight-binding model difficult. In any case, the Wannier functions we obtain by projection have a possible improvement in total spread of less than 10% (obtained by calculating both the gaugedependent and gauge-independent parts of [37]), meaning the projected orbitals are almost maximally localized anyway.

C. Wannier projection
An alternate approach to constructing Wannier functions is by the projection of M Bloch states ψ mk (r) onto N trial guesses g n (r) [37,38,42]. This is done by first defining the inner products where w mk ∈ [0, 1] is an optional weighting of the Bloch states, which can be tuned to allow for smooth band disentanglement. Next, one computesã, the closest unitary approximation of a (for now suppressing k dependence). In practice, this is done by first factorizing a through the singular value decomposition (SVD), a = W V † , where W, V are unitary matrices and is diagonal, and then setting Note that this is equivalent to definingã = a(a † a) −1/2 . Then, to obtain the WFs one simply replaces the matrix U (k) mn in Eq. (3) withã (k) mn . For this work, the choice of trial functions g(r) and the weighting of the Bloch states w mk are very important. If the trial guess and the Hilbert subspace are not compatible enough, the SVD will be unable to create a good unitary approximation of the a (k) mn matrix, resulting in poorly localized Wannier functions and failures in the tight-binding band structure. The functions g(r) can be initial guesses, or created iteratively by doing multiple Wannier projections of the same (or similar) Hamiltonians. The weightings w mk can depend on band index or eigenvalue energy, or ignored if the number of selected bands is identical to the number of targeted Wannier functions. We will later highlight how each of these choices of g(r) and w mk are used in various stages of our complete projection methodology.

D. Choice of initial trial functions
Following Ref. [19], the five-band and eight-band models must satisfy a representation matching equation to resolve the fragile topology of the flat bands. To write down this equation, one must pick a lattice site and orbital symmetry for each Wannier function. The candidate effective lattice sites are: (1) triangular, labeled τ and corresponding to AA stacking regions, (2) honeycomb, labeled η and corresponding to AB/BA stacking regions, and (3) kagome, labeled κ and corresponding to intermediate stacking regions between AB and BA sites. After considering relaxation of the TBLG system [13,28,[43][44][45], for θ < 1.5 • the honeycomb sites become large triangular domains of AB or BA stacking, and the 033072-3 kagome lattice locations can be interpreted as domain walls (DW ) between them.
Along with these three sets of sites, we also consider four possible orbital symmetries, s, p z , p + , and p − (p ± = p x ± ip y ). These correspond to the same symmetry eigenvalues expected of a hydrogen orbital of the same label, for example the s orbital is symmetric under z → −z, but the p z orbital is antisymmetric. In the case of TBLG, the corresponding operation is physically implemented by a layer-exchanging twofold rotation in three-dimensional space about an axis running parallel to the layers.
With this groundwork covered, the representation matching equations can be summarized. We use (f.b.) to represent the desired flat bands. For the eight-band model, For the five-band model including three bands above the flat bands, which we call five-band (top), the proper equation is and for the five-band model including three bands below the flat bands, which we call five-band (bot), the proper equation In every equation, the left-hand side corresponds to a realizable tight-binding model, while the right-hand side is a combination of topologically trivial states and the flat bands with fragile topological obstruction. The two five-band models are related by swapping p z ↔ s, and are defined by the -point mirror-symmetry eigenvalues of the realistic bands just above or below the flat bands.
These representation-matching equations prescribe the correct g(r) trial guesses for a successful Wannier projection. We build the guesses with a tight Gaussian wavepacket at the desired lattice site, and then assign the proper symmetry by hand by manipulating the phase of the layer and sublattice channels. For example, to create a p z type orbital, one creates a Gaussian wave packet with positive sign on the top layer and negative sign on the bottom layer. These simple Gaussian trial functions change as they are projected into the Hilbert subspace of the low-energy band structure. After successful projection, the resulting WFs have lattice centers and symmetry identical to their trail guesses, but they yield new interpretation in terms of stacking order and layer/sublattice character. In doing so, we find that the symmetry-motivated prescription for each model can be connected to a physical "moire orbital" motivated by the local stacking order. The symmetries and interpretations of the WFs of each model are summarized in Table I.

E. Multistep projection
The desired model will have perfectly accurate band structure for the flat-band manifold, with the accuracy of the band energies decreasing as one moves further away from the Fermi energy. Our target Wannier functions therefore require a very careful projection of the low-energy band structure in order to controllably put most of the projection error in the higher-energy bands. We remark that, as the higher-energy TABLE I. Orbitals of the reduced Hamiltonians. The Wannier function index (No.), lattice description, and symmetry properties for the eight-band model are given in the first three columns, and describe both the initial trial Gaussian functions as well as the resulting Wannier functions. The symbols τ , η, and κ correspond to a triangular, hexagonal, and kagome lattices, respectively. The symbols p ± , p z , and s refer to the symmetries of a hydrogenic wave function (see text for more details). The fourth column gives a description of the Wannier function character, including relations to stacking order and any layer or sublattice polarization. The last two columns list the symmetry content of the five-band top and bottom models, with the rest of their Wannier functions' information identical to the eight-band case.
Lat. Sym. Resulting orbital description 5b top 5b bot bands do not separate into well-defined groups defined by band gaps, certain arbitrariness will inevitably exist in how we disentangle these higher-energy bands. Here, we employ the projection technique [38] to select states that are smooth enough to allow for wannierization. Furthermore, in order to generate a full series of effective models over a large span of twist angles we eliminate all manual tuning in the process. We achieve this automated construction with multiple Wannier projection steps, taking advantage of the chiral-symmetric limit [29] and improving the projection in an iterative fashion. The entire process is schematically displayed in Fig. 1.
Step 1. Use N prescribed Gaussian trial guesses on the subspace defined by exactly N bands of the chiral-symmetric (CS) model [29] to generate Wannier functions of the chiral symmetric model, WF cs . This is only possible in the CS model because the flat bands and the three-band manifolds above and below the flat bands do not intersect with any other bands over the entire BZ. Also perform a projection of only the threeband manifolds in the CS model to generate corresponding three-band WFs.
Step 2. Use the three-band WFs as trial guesses for an energy-weighted projection of the realistic model, which disentangles a three-band subspace from the multiple band crossings away from the flat-band manifold. The projection is tuned depending on the angle, but always has full weight near the flat bands, with exponentially decaying edges in the center of the three-band manifolds. Combine the disentangled threeband subspaces with the flat-band Bloch states to generate a robust Hilbert subspace for the realistic model which still ensures perfect agreement along the flat bands. Use WF CS as trial guesses on the subspace just obtained to get the first Wannier functions for the realistic model, WF 0 . For angles near the magic angle (roughly 0.85 • to 1.2 • ), we find that WF 0 is already a very good model for the low-energy structure, but can have some band structure errors near the point. The following step helps remove these problems, but is generally not necessary when there are sizable band gaps between the flat-band manifold and the three-band manifold above and below the flat bands.
Step 3. Use WF 0 as trial WFs on an unweighted subspace of the ab initio low-energy bands, usually defined by a hard energy cutoff. This cleans up spurious band hybridizations at low angles for the eight-band model, and generally makes the band disentanglement more accurate for both the eight-and five-band models. We call these improved Wannier functions WF 1 .
In general, all of these projections can be stored in k-space or real space. The k-dependent n × n models are defined over a Brillouin zone sampling of the moire supercell, so one can easily integrate the result for different Wannier centers to generate an n × n real space model. The k-space models are sometimes difficult to transfer between the chiral-symmetric and realistic model, especially if one uses a different trun-cation radius for the two models (the CS model tends to converge much faster). For this reason, the resulting Wannier projections are usually stored in real space on a uniform grid sampling of the moire superlattice. By sampling in this manner, one can in principle use the WF 1 generated at some angle θ as a trial guess for a one-shot projection at θ + δ. We find this is not always a stable process, and instead use the full multistep projection process over a range of different angles.
As an aside, we note that the flat bands of twisted double bilayer graphene (TDBLG) [46][47][48][49][50] are much simpler to Wannierize properly. In TDBLG, the C 2 symmetry is broken, which prevents its flat bands from entangling with the rest of the band structure, and thus one can accurately discuss them without invoking additional bands. When there is no Chern number for the flat bands, a one orbital model for the active band can be constructed. In the presence of a Chern number, a Haldane like model for the two opposite Chern number bands can instead be built. This is a worthwhile exercise but relatively straightforward and does not require the careful multistep projection technique introduced here.

F. Redistribution of orbital character
While the projection technique described allows us to construct reduced models with the desired properties, the definition of the effective orbitals still has some arbitrariness which allows further refinement. More precisely, within the reduced Hilbert space one can still perform a basis rotation and redefine each of the effective orbitals (subjected to the symmetry constraints). This should be contrasted with conventional solid-state systems, in which the notion of atomic orbitals is fixed by the microscopic lattice.
It is expected that the most important interaction terms for TBLG's flat bands arise due to the electron density that forms on the AA stacking sites near the charge neutrality point. Due to our symmetry prescription, these areas of high electron density can be best captured by the AA ± orbitals in our models. In anticipation of eventually incorporating electron-electron interaction, it is desirable to utilize the remaining freedom in basis transformation to maximize the weights of the AA ± orbitals (WF Nos. 1 and 2) in the nearly flat bands. Doing so can simplify interpretations of the many-body physics and, as a practical consideration, reduce the number of two-body and four-body integrals required to properly define the many-body model.
To achieve the desired weighting, we perform a final "rewannierization" of the reduced models. Note that this is done entirely within the reduced basis obtained from the previous subsection. We take our initial guesses to be the original Wannier Function basis (so g n = e n , a basis element of H), and the subspace is the entire reduced band structure.
For the sake of argument, imagine that the nearly flat bands have the symmetry and topology of the AA ± orbitals. If this was true, we should simply construct Wannier functions for the flat band manifold alone and relabel them as the new AA ± orbitals. In reality, there is a topological obstruction to such a construction, and by considering the symmetry representations we know that the overlap between the nearly flat bands and the AA ± orbitals has to vanish at the point. In other words, the weights of the AA ± orbitals will necessarily 033072-5 have to be transferred to the higher energy bands at (and near) the point. This motivates us to define a weighting of the Bloch states, w mk , which depend on both the band index and the momentum. We let the weighting for the two flat bands be 1 at all k, but for the two bands above the flat bands, as well as the two below, we incorporate them into the projection with weights w mk = e −(k/σ k ) 2 , i.e., they are only incorporated in the vicinity of the point. For every k, we perform a two-orbital projection onto the subspace spanned by the original AA ± orbitals, which allows us to construct a new set of Wannier functions of the same symmetry character. The remaining WFs are then reconstructed by a projection of the orthogonal subspace onto the remaining orbitals.
An example of the rewannierization of the eight-band model is presented in Fig. 2. One can see in both the projected band structures and the projected density of states (PDoS) that the flat bands in the original reduced model ("Before") do not have a clear orbital character, except for certain highsymmetry points at which symmetry constraints are in play. After the rewannierization ("After"), the flat bands are almost entirely AA ± character in the new basis, except for the small "leakage" at the point necessitated by the symmetry constraints. Likewise, the other Wannier functions only have a small weighting in the flat bands, confined near the point. Overall, the flat bands go from 58% AA ± character to 92% AA ± character, a significant improvement that ensures the dominant interaction terms relevant to the correlated phenomena in the flat bands will involve mostly the AA ± Wannier functions. As the rewannierization is simply a change of basis within the reduced Hilbert space, it has no effect on the band structure. Furthermore, we find that it does not increase nor decrease the reduced model hopping range for suitable choices of the weighting parameter σ k . We note in passing that, in principle, if different orbital character was desired in various bands, the above procedure could be modified to create such a model, but we will not pursue this point in our reduced models.

A. Reduced tight-binding Hamiltonians
An overview of the three models obtained from the multistep projection and targeted rewannierization is shown for θ = 1.1 • in Fig. 3. All three models perfectly reproduce the flat-band manifold of the k · p model, and also capture the relevant band gap(s). Even the band dispersion of the first two bands above and below the flat-band manifold are robustly captured. Since these models are based on a single valley in a k · p expansion, the underlying model lacks proper timereversal symmetry, and as such the tight-binding Hamiltonians are generically complex. Combining two copies of the Hamiltonian, one from the monolayer K (already obtained) and one from K (which can be derived from K), would yield a tight-binding Hamiltonian that can be written with purely real coefficients after a linear transformation of the basis elements. Such a model would have band structure identical to a fully atomistic tight-binding treatment, but naturally with twice as many orbitals as the single-valley model. Approximating the intervalley coupling to be zero, there is usually no reason to use the two-valley tight-binding model. If intervalley coupling is desired, it can be easily added in an empirical fashion to this two-valley model. The Wannier functions have the proper symmetry, and most of their respective density remains close to their center (unlike the "fidget spinner" Wannier functions of the previous effective tight-binding models [20,22,24,25]). The layer and sublattice projections of these functions is summarized in the last column of Table I. The images of the wave functions in Fig. 3 show only their magnitudes averaged over both layer and sublattice indices, but we have checked that they have underlying index polarization and complex phase consistent with their symmetry descriptions. These interpretations of the model are similar to our previous and less-sophisticated treat-ment of a ten-band model [27]. We note that the interpretation of the projected density of states presented in the ten-band study are related, but not identical, to the electronic structure here. In this model, the protected s versus p z band-crossing of the flat bands at the magic-angle still exists, but now one of the flat bands at the point is a mixture of AA and DW (s) character, while the other is predominately AB (p z ) character. This is due to changes in the symmetry prescription between the 10 and 8 band models.
To accurately reproduce the flat bands, the five-band Hamiltonian needs to retain coupling terms only up to two 033072-7 moire lengths (2λ θ ) and the eight-band needs only up to three (3λ θ ). Thus for studying flat-band correlations, the five-band models are preferable, but for studying optical properties (or for correlations at twist angles well away from the magic angle) the eight-band model is still usable. Our approach differs from previous empirical treatment of a similar fiveband model (see Appendix of Ref. [19]), as we obtain a generic angle dependence of all tight binding terms and couplings beyond nearest neighbor. Overall, these models are too complicated to give here in detail. Instead, we provide tables of the values of the Hamiltonian as a function of twist angle, and tools to generate the reduced models presented here from the realistic k · p model [51].

B. Twist angle dependence
The twist angle dependence of the eight-band model is summarized in Fig. 4. We find that the automated multistep projection technique is able to capture the low-energy band structure for θ ∈ [0.6 • , 1.3 • ]. It fails outside of this angle range because of band intersections in the chiral-symmetric model. For θ > 1.3 • , the three-band manifolds above and below the flat bands cross with the bands further from the Fermi energy (the thin black bands near in the left panel of Fig. 1). Similarly, for θ < 0.63 • , the flat bands begin to mix with the three-band manifolds above and below the flat bands, meaning the construction of separate two-band and three-band manifolds is no longer well defined. The construction of Wannier functions past these twist angles is not impossible, but a projection technique different than the one we present here is necessary.
The dependence of the spread of the individual Wannier functions is also shown in Fig. 4(b). We see that the Wannier functions 1-3 (on the triangular lattice) do not scale like the moire length λ θ = 1/θ near the magic angle. This is because for θ 1 • , atomic relaxations cause the formation of AA stacking regions that are fixed in size even as the angle decreases [13,28,[43][44][45]. The AA stacking regions can be thought of as the intersection of the 1D strain solitons on the domain walls, and their size is determined through a balancing of in-plane strain energy to interlayer stacking energy [44]. The other orbitals on the kagome lattice (domain walls) and Hexagonal lattice (AB/BA) increase in spread as the twisting-angle decreases, which is not surprising as they are pinned to the geometry of the expanding domains and domain-walls. When the gap between the flat band manifold and the nearby band manifolds closes (θ < 0.85 • ) we see Wannier functions 1 and 2 (AA ± ) begin to grow in size, due to increasing hybridization with electronic structures not associated with the AA stacking spots.

C. Representing symmetry lowering perturbations in the reduced model
Two of the important symmetries of bilayer graphene system can be easily broken in experimental devices: the layer symmetry and sublattice symmetry. Layer symmetry can be lifted by applying an external vertical electric displacement field. This places the top and bottom layers at different electrostatic potentials, and to first order can be included in a k · p model by adding an onsite energy difference ±δ E for the top and bottom layer degrees of freedom (determined by the magnitude and sign of E ).
The sublattice symmetry can be broken experimentally when the bilayer is encapsulated by nearly-aligned hexagonal boron nitride (hBN). hBN has a similar crystal structure to graphene, but is composed of boron and nitrogen atoms instead of only carbon. When graphene is in close proximity of aligned hBN, sublattice symmetry is broken as a carbon atom near boron feels a slightly different electrostatic potential than for the one near nitrogen. This can be similarly included in a k · p model by including a ±δ S for the A and B sublattice degrees of freedom (determined by α or β alignment with hBN). In general, this energy difference does not have to be the same for the two layers, as physically the top and bottom hBN encapsulations do not have to be simultaneously aligned in the same way. To study any combination, we consider both layer symmetric sublattice energy, δ S , and layer 033072-8  Fig. 3 for the eight-band model. The red arrows correspond to large matrix elements that "hop" into a Wannier function of the fundamental moire cell. The circular arrow on H S represents a large onsite term for Wannier functions 1 and 2 (AA ± ), but not 3 (AA s ).
anti-symmetric sublattice energy, δ S . For example, the layer symmetric perturbation means the A (B) sublattice of both layers feels a positive (negative) onsite perturbation, while the layer antisymmetric means the A (B) sublattice of layer 1 is positive (negative) while the A (B) sublattice of layer 2 is negative (positive).
For both types of perturbations, the new k · p Hamiltonian is given by H k = H 0 k + δ i H k . One can easily calculate the first order perturbation H for the reduced Wannier models by taking the n-band subspace projector P n and sandwiching it around the perturbation of H k : H = P † n H k P n . The k · p band structures and a summary of the largest terms of H i are presented in Fig. 5.
If the two graphene layers were decoupled, the only effect of electric displacement field would be to move the Dirac cones of each layer away from one another in energy. This can be seen in the k · p band structures as a valley-dependent energy shift of the states at the moiré K points (the red and blue bands are no longer identical near K), which corresponds to a breaking of the layer-exchanging symmetries. The eightband reduced model captures this effect by having mirror-symmetry breaking terms added between the triangular and kagome lattice orbitals.
For a graphene monolayer, the substrate symmetry breaking term opens up a finite gap at the Fermi energy, splitting the Dirac cone into two parabolas. In the bilayer system, the gap remains and is roughly equal in magnitude to the value of δ S . In the reduced eight-band Hamiltonian, the layer-symmetric sublattice perturbation manifests as ideal onsite energies for Wannier functions 1 and 2 (0.5δ S and −0.5δ S , respectively), as well as some smaller coupling terms between the honeycomb and kagome sites. This is not surprising, as WFs 1 and 2 have the highest amount of sublattice polarization, while WF 3 has the least. We remark that, although the Dirac cones are gapped by an onsite energy difference, the resulting nondegenerate band (in each valley) carries a nonzero Chern number. This is a direct consequence of the topological character of the flat bands in TBLG [22]. The anti-symmetric perturbation δ S is not shown in Fig. 5, but it has a similar band structure to δ S with a smaller gap opening (roughly 30% of that for δ S ). The largest terms are couplings between orbitals 1, 2 and 6, 7, 8 (e.g., AA to DW couplings), but they have a smaller overall value of 0.15δ S .

IV. CONCLUSION
We have presented a robust projection scheme for obtaining accurate Wannier functions in twisted bilayer graphene. The construction of these Wannier functions takes advantage of the chiral-symmetric limit of TBLG's low-energy model, and can create Hamiltonians that reproduce realistic band structures even outside the magic angle regime. Although we did not maximize the localization of these Wannier functions, their spread is close to the smallest possible for the chosen Hilbert subspace.
The Wannier functions exist on a triangular, honeycomb, and kagome lattice. This yields a direct interpretation in terms of local stacking order, allowing one to ascribe features of the electronic structure to interactions between effective orbitals defined over the moire superlattice. The low-energy Hamiltonians can be truncated symmetrically, and they need hoppings with lengths only up to three moire lengths to accurately reproduce the flat bands near the magic angle. Studying the Wannier functions as the twist angle changes shows that the functions defined on the AA stacking (orbitals 1-3) do not scale like the moire length λ θ while those defined on the AB or DW geometry scale like λ θ , consistent with the features of atomistic relaxation in the system. The effect of external electric displacement field and sublattice symmetry breaking cause clear changes to the low-energy band structure, and these changes are captured in the reduced model by changing different subsets of hopping terms.
The models presented here can be used as a starting point for unbiased study of correlated phases in TBLG. Our models differ from previous approaches primarily in the resolution of TBLG's fragile topology: instead of breaking a symmetry we include additional bands outside of the flat band manifold. This allows for the comparison of symmetry and symmetrybroken states and removes the unconventional long-range coupling of the propeller orbitals. Furthermore, with complete 033072-9 freedom in the choice of twist angle, one can illuminate the nature of correlations near the magic angle and also predict the presence of other (similar or dissimilar) correlated phases in the range [0.6 • , 1.3 • ]. We leave such study to future work. The reduced models are made publicly available as part of a larger suite of continuum models for relaxed TBLG [51]. These codes include the eight-band and five-band Hamiltonians for various angles highlighted in this work and also allows for the construction of additional models with arbitrary twist angle and external perturbations.