Derivation of Wannier orbitals and minimal-basis tight-binding hamiltonians for twisted bilayer graphene: a first-principles approach

TL;DR

This paper develops minimal-basis tight-binding models and Wannier orbitals for twisted bilayer graphene, accurately capturing low-energy physics and symmetries, facilitating studies of correlated phenomena.

Contribution

It introduces a first-principles multi-step Wannier projection method to derive minimal models and orbitals for TBLG, including effects of relaxations and symmetries.

Findings

Constructed eight-band and five-band models for twist angles 1.3° to 0.6°.

Derived Wannier orbitals incorporating all TBLG symmetries.

Models reproduce low-energy bands of ab initio tight-binding calculations.

Abstract

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 \textit{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^\circ$ and $0.6^\circ$. The models are created using a multi-step Wannier projection technique starting with appropriate $\textit{ab initio}$ $k \cdot 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.

Tables

Table 1. Table 1: Orbitals of the reduced Hamiltonians. The Wannier function index (#), lattice description, and symmetry properties for the 8-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 τ 𝜏 \tau , η 𝜂 \eta and κ 𝜅 \kappa correspond to a triangular, hexagonal, and Kagome lattice, respectively. The symbols p ± subscript 𝑝 plus-or-minus p_{\pm} , p z subscript 𝑝 𝑧 p_{z} , and s 𝑠 s refer to the symmetries of a hydrogenic wavefunction (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 5-band top and bottom models, with the rest of their Wannier functions’ information identical to the 8-band case.

#	Lat.	Sym.	Resulting orbital description	5b top	5b bot
1	$\tau$	$p_{+}$	$B$ orbital dimers at $AA$	$p_{+}$	$p_{+}$
2	$\tau$	$p_{-}$	$A$ orbital dimers at $AA$	$p_{-}$	$p_{-}$
3	$\tau$	$s$	Equal mixture around $AA$	$s$	$p_z$
4	$\eta$	$p_z$	$A$ on $L_1$, $B$ on $L_2$ at $AB$	$p_z$	$s$
5	$\eta$	$p_z$	$B$ on $L_1$, $A$ on $L_2$ at $BA$	$p_z$	$s$
6	$\kappa$	$s$	$AB$ ($BA$) above (below) $DW$	$-$	$-$
7	$\kappa$	$s$	$AB$ ($BA$) left (right) of $DW$	$-$	$-$
8	$\kappa$	$s$	$BA$ ($AB$) left (right) of $DW$	$-$	$-$

Equations

$$\Omega =\sum _n\left[⟨{\varphi }_n|r^2|{\varphi }_n⟩-(⟨{\varphi }_n|\mathbf{r}|{\varphi }_n⟩{)}^2\right].$$

$$M_{mn}^{\mathbf{k},\mathbf{b}}=⟨u_m^{(\mathbf{k})}⟩u_n^{(\mathbf{k}+\mathbf{b})}$$

$$|\mathbf{R}n⟩=\frac{V}{(2\pi {)}^3}{\int }_{BZ}d\mathbf{k}{e}^{-i\mathbf{k}\cdot \mathbf{R}}\sum _{m=1}^N{U}_{mn}^{(\mathbf{k})}|{\psi }_{m\mathbf{k}}⟩$$

$$a_{mn}^{(\mathbf{k})}=w_{m\mathbf{k}}⟨{\psi }_{m\mathbf{k}}⟩g_n,$$

$$\tilde{a}=a\left[V\frac{1}{\sqrt{{\Sigma }^{†}\Sigma }}{V}^{†}\right]$$

$$(\tau ,p_{\pm })\oplus (\tau ,s)\oplus (\eta ,p_z)\oplus (\kappa ,s)=(\kappa ,s)\oplus (\kappa ,p_z)\oplus (\text{f.b.})$$

$$(\tau ,p_{\pm })\oplus (\tau ,s)\oplus (\eta ,p_z)=(\kappa ,p_z)\oplus (\text{f.b.})$$

$$(\tau ,p_{\pm })\oplus (\tau ,p_z)\oplus (\eta ,s)=(\kappa ,s)\oplus (\text{f.b.})$$

Full text

††preprint: APS/123-QED

I Introduction

Recent experimental advances in nanotechnology fabrication methods Frisenda et al. (2018) 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 $\theta$) 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 propertiesLi et al. (2009); Brihuega et al. (2012); Luican et al. (2011); Wong et al. (2015) – in essence introducing the new field of “Twistronics” Carr et al. (2017). One of the most well known of twist-angle dependent phenomena is the emergence of correlated insulators and unconventional superconductivity in TBLG at $\theta\approx 1^{\circ}$, with a superconducting $T_{c}$ of a few degrees Kelvin Cao et al. (2018a, b); Yankowitz et al. (2019). Theoretical studies of TBLG near the critical (often referred to as “magic”) angle of $\sim 1^{\circ}$ 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 Trambly de Laissardière et al. (2010); Uchida et al. (2014), tight-binding hamiltonians Suárez Morell et al. (2010); Carr et al. (2017); Nam and Koshino (2017); Angeli et al. (2018), and continuum models Mele (2010); Bistritzer and MacDonald (2011); San-Jose et al. (2012); Weckbecker et al. (2016).

All such models feature a large number of degrees of freedom, which render the study of electrons’ correlation unfeasible. Even the simplest, a continuum $\bm{k}\cdot\bm{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 Po et al. (2018a); Koshino et al. (2018); Zou et al. (2018); Po et al. (2018b); Yuan and Fu (2018); Kang and Vafek (2018); Goodwin et al. (2019). These models are generally empirical, and chosen to address the flat bands near the magic-angle 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 Nam and Koshino (2017); Carr et al. (2018a). 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 $\theta=1^{\circ}$. Our starting point is a microscopic $\bm{k}\cdot\bm{p}$ model for TBLG which includes atomic relaxations Carr et al. (2019, 2018b), from which we derive three single-valley tight-binding models. The first, and more comprehensive, is an $8$-band model, which includes the 2 flat bands and 3 bands above and 3 below them. Two additional reduced models include $5$ bands, and involve the 2 flat bands and the nearest three bands either above or below them. The two $5$-band models are only reliable when the flat band manifold is isolated in the band structure. Near $0.85^{\circ}$, the gaps above and below the flat bands vanish, and only the $8$-band model continues to be trustworthy. We find the $8$-band model can be reliably constructed down to $\theta=0.6^{\circ}$, with only technical limitations of the projection technique restricting further reduction in twist angle.

The model construction is performed with a multi-step Wannier projection technique, which utilizes the chiral-symmetric $\bm{k}\cdot\bm{p}$ model Tarnopolsky et al. (2019) 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\times 4$) or 32 band ($8\times 4$) model is challenging within any exact multi-particle 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.

II Methods

II.1 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 non-trivial windings which lead to the topological Chern number Xiao et al. (2010); Thouless et al. (1982). Such topological character is known to obstruct the construction of real space Wannier states Thonhauser and Vanderbilt (2006). 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 Soluyanov and Vanderbilt (2011).

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}$) Zou et al. (2018), one cannot have a real space two-band representation of the flat bands Po et al. (2018b); Zou et al. (2018). This obstruction can be understood from the net chirality of two Dirac cones or from the symmetry representations of the flat bands Po et al. (2018b); Zou et al. (2018). 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 Po et al. (2018b) or the two-fold rotation symmetry Yuan and Fu (2018). 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 propeller-like (or ”fidget spinner”) shape, created by the superposition of three charge pockets at neighboring AA-stacking spots Po et al. (2018b); Koshino et al. (2018); Kang and Vafek (2018); Goodwin et al. (2019). 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. Zou et al. (2018); Song et al. (2019) which doubled the bands. This new set of topologically non-trivial 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. Po et al., 2018a, it has been pointed out that the root for the obstruction in TBLG’s flat bands is a “fragile topology” Po et al. (2018c); Cano et al. (2018). This topological obstruction is relatively weak (fragile) when compared to the more robust Chern index, in the sense that the non-triviality 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 are more conventional and short range as compared to the propeller solutions.

II.2 Maximally Localized Wannier Functions

A common approach to generate a localized model is to construct maximally localized Wannier functions (MLWF) Marzari and Vanderbilt (1997); Marzari et al. (2012). They define tight-binding wavefunctions, $\phi_{n}(\bm{r})$, which minimize the total wavefunction spread

[TABLE]

The simplest construction of MLWF for $N$ bands relies on the calculation of $M\geq N$ Bloch states $\psi_{m\bm{k}}$ over a uniform sampling of the BZ, and then evaluation of the overlap matrix elements of the periodic parts of the Bloch states, $u_{m}^{(\bm{k})}$ ,

[TABLE]

where $\bm{b}$ is the momentum connecting each $\bm{k}$ to its nearest neighbor on the grid sampling of the BZ. The overlap elements are then used to calculate a set of $\bm{k}$-dependent unitary rotations $U_{mn}^{(\bm{k})}$ that gives the proper mixing of the original Bloch states over the BZ to minimize the real space spread Marzari and Vanderbilt (1997). Then the $n^{\textrm{th}}$ Wannier function at unitcell $\bm{R}$ is given by

[TABLE]

Applying this methodology to a $\bm{k}\cdot\bm{p}$ model is not entirely trivial however. If one attempts to construct a sampling of the $\bm{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 $n^{\textrm{th}}$ basis element of the continuum Hamiltonian at $\bm{k}$ will not correspond to the $n^{\textrm{th}}$ basis element at $\bm{k}+\bm{\delta}$ if $\bm{\delta}$ 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^{\bm{k},\bm{b}}_{mn}$ even when $\bm{k+b}$ crosses a BZ boundary.

Even with these choices, the 8-band and 5-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 Mostofi et al. (2008)) results in Wannier functions that have no clear symmetry within the moiré 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) Wang et al. (2014) 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 Sakuma (2013) 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 wavefunctions 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 $\Omega$ 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 gauge-dependent and gauge-independent parts of $\Omega$ Marzari and Vanderbilt (1997)), meaning the projected orbitals are almost maximally localized anyway.

II.3 Wannier Projection

An alternate approach to constructing Wannier functions is by the projection of $M$ Bloch states $\psi_{m\bm{k}}(\bm{r})$ onto $N$ trial guesses $g_{n}(\bm{r})$ Marzari and Vanderbilt (1997); Marzari et al. (2012); Mustafa et al. (2015). This is done by first defining the inner products

[TABLE]

where $w_{m\bm{k}}\in[0,1]$ is an optional weighting of the Bloch states, which can be tuned to allow for smooth band disentanglement. Next, one computes $\tilde{a}$, the closest unitary approximation of $a$ (for now suppressing $\bm{k}$ dependence). In practice, this is done by first factorizing $a$ through the singular value decomposition (SVD), $a=W\Sigma V^{\dagger}$, where $W,V$ are unitary matrices and $\Sigma$ is diagonal, and then setting

[TABLE]

Note that this is equivalent to defining $\tilde{a}=a(a^{\dagger}a)^{-1/2}$. Then, to obtain the WFs one simply replaces the matrix $U_{mn}^{(\bm{k})}$ in Eq. 3 with $\tilde{a}_{mn}^{(\bm{k})}$.

For this work, the choice of trial functions $g(\bm{r})$ and the weighting of the Bloch states $w_{m\bm{k}}$ 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_{mn}^{(\bm{k})}$ matrix, resulting in poorly localized Wannier functions and failures in the tight-binding band structure. The functions $g(\bm{r})$ can be initial guesses, or created iteratively by doing multiple Wannier projections of the same (or similar) Hamiltonians. The weightings $w_{m\bm{k}}$ 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(\bm{r})$ and $w_{m\bm{k}}$ are used in various stages of our complete projection methodology.

II.4 Choice of Initial Trial Functions

Following Ref. Po et al., 2018a, the 5-band and 8-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 $\tau$ and corresponding to $AA$ stacking regions, (2) honeycomb, labeled $\eta$ and corresponding to $AB/BA$ stacking regions, and (3) Kagome, labeled $\kappa$ and corresponding to intermediate stacking regions between $AB$ and $BA$ sites. After considering relaxation of the TBLG system Dai et al. (2016); Nam and Koshino (2017); Zhang and Tadmor (2018); Carr et al. (2018b); Yoo et al. (2019), for $\theta<1.5^{\circ}$ the honeycomb sites become large triangular domains of $AB$ or $BA$ stacking, and the 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_{\pm}$ = $p_{x}\pm 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\to-z$, but the $p_{z}$ orbital is anti-symmetric. In the case of TBLG, the corresponding operation is physically implemented by a layer-exchanging two-fold 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 $(\textrm{f.b.})$ to represent the desired flat-bands. For the 8-band model:

[TABLE]

For the 5-band model including three bands above the flat-bands, which we call 5-band (top), the proper equation is:

[TABLE]

and for the 5-band model including three bands below the flat-bands, which we call 5-band (bot), the proper equation is:

[TABLE]

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 5-band models are related by swapping $p_{z}\leftrightarrow{}s$, and are defined by the $\Gamma$-point mirror-symmetry eigenvalues of the realistic bands just above or below the flat-bands.

These representation-matching equations prescribe the correct $g(\bm{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 wavepacket 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 “moiré orbital” motivated by the local stacking order. The symmetries and interpretations of the WFs of each model are summarized in Table 1.

II.5 Multi-step 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 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 Marzari et al. (2012) 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 Tarnopolsky et al. (2019) 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 Tarnopolsky et al. (2019) to generate Wannier functions of the chiral symmetric model, $\textrm{WF}_{\textrm{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 three-band 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 3-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 3-band 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 $\textrm{WF}_{\textrm{CS}}$ as trial guesses on the subspace just obtained to get the first Wannier functions for the realistic model, $\textrm{WF}_{0}$. For angles near the magic angle (roughly $0.85^{\circ}$ to $1.2^{\circ}$), we find that $\textrm{WF}_{0}$ is already a very good model for the low-energy structure, but can have some band structure errors near the $\Gamma$ 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 $\textrm{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 8-band model, and generally makes the band disentanglement more accurate for both the 8 and 5-band models. We call these improved Wannier functions $\textrm{WF}_{1}$.

In general, all of these projections can be stored in $\bm{k}$-space or real space. The $\bm{k}$-dependent $n\times n$ models are defined over a Brillouin zone sampling of the moiré supercell, so one can easily integrate the result for different Wannier centers to generate an $n\times n$ real space model. The $\bm{k}$-space models are sometimes difficult to transfer between the chiral-symmetric and realistic model, especially if one uses a different truncation 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 moiré superlattice. By sampling in this manner, one can in principle use the $\textrm{WF}_{1}$ generated at some angle $\theta$ as a trial guess for a one-shot projection at $\theta+\delta$. 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) Shen et al. (2019); Cao et al. (2019); Liu et al. (2019); Chebrolu et al. (2019); Koshino (2019) 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 multi-step projection technique introduced here.

II.6 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_{\pm}$ 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_{\pm}$ orbitals (WF #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 ”re-Wannierization” 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_{\pm}$ orbitals. If this was true, we should simply construct Wannier functions for the flat band manifold alone and relabel them as the new $AA_{\pm}$ 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_{\pm}$ orbitals has to vanish at the $\Gamma$ point. In other words, the weights of the $AA_{\pm}$ orbitals will necessarily have to be transferred to the higher energy bands at (and near) the $\Gamma$ point. This motivates us to define a weighting of the Bloch states , $w_{m\bm{k}}$, which depend on both the band index and the momentum. We let the weighting for the two flat bands be $1$ at all $\bm{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_{m\bm{k}}=e^{-(k/\sigma_{k})^{2}}$, i.e., they are only incorporated in the vicinity of the $\Gamma$ point. For every $\bm{k}$, we perform a two-orbital projection onto the subspace spanned by the original $AA_{\pm}$ 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 re-Wannierization of the 8-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 high-symmetry points at which symmetry constraints are in play. After the re-Wannierization (”After”), the flat bands are almost entirely $AA_{\pm}$ character in the new basis, except for the small ”leakage” at the $\Gamma$ point necessitated by the symmetry constraints. Likewise, the other Wannier functions only have a small weighting in the flat bands, confined near the $\Gamma$ point. Overall, the flat bands go from 58% $AA_{\pm}$ character to 92% $AA_{\pm}$ character, a significant improvement that ensures the dominant interaction terms relevant to the correlated phenomena in the flat bands will involve mostly the $AA_{\pm}$ Wannier functions. As the re-Wannierization 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 $\sigma_{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.

III Results

III.1 Reduced Tight Binding Hamiltonians

An overview of the three models obtained from the multi-step projection and targeted re-Wannierization is shown for $\theta=1.1^{\circ}$ in Fig. 3. All three models perfectly reproduce the flat-band manifold of the $\bm{k}\cdot\bm{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 $\bm{k}\cdot\bm{p}$ expansion, the underlying model lacks proper time-reversal 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^{\prime}$ (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 inter-valley coupling to be zero, there is usually no reason to use the two-valley tight-binding model. If inter-valley 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 Po et al. (2018b); Koshino et al. (2018); Kang and Vafek (2018); Goodwin et al. (2019)). The layer and sublattice projections of these functions is summarized in the last column of Table 1. The images of the wavefunctions 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 treatment of a $10$-band model Carr et al. (2019). We note that the interpretation of the projected density of states presented in the $10$-band study are related, but not identical, to the the electronic structure here. In this model, the protected $s$ vs $p_{z}$ band-crossing of the flat-bands at the magic-angle still exists, but now one of the flat-bands at the $\Gamma$ 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 5-band Hamiltonian needs to retain coupling terms only up to two moiré lengths ($2\lambda_{\theta}$) and the 8-band needs only up to three ($3\lambda_{\theta}$). Thus for studying flat-band correlations the 5-band models are preferable, but for studying optical properties (or for correlations at twist angles well away from the magic-angle) the 8-band model is still usable. Our approach differs from previous empirical treatment of a similar 5-band model (see Appendix of Ref. Po et al., 2018a), as we obtain 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 $\bm{k}\cdot\bm{p}$ model not .

III.2 Twist angle dependence

The twist angle dependence of the 8-band model is summarized in Fig. 4. We find that the automated multi-step projection technique is able to capture the low-energy band structure for $\theta\in[0.6^{\circ},1.3^{\circ}]$. It fails outside of this angle range because of band intersections in the chiral-symmetric model. For $\theta>1.3^{\circ}$, the three-band manifolds above and below the flat bands cross with the bands further from the Fermi energy (the thin black bands near $\Gamma$ in the left panel of Fig. 1). Similarly, for $\theta<0.63^{\circ}$, the flat bands begin to mix with the the three-band manifolds above and below the flat bands, meaning the construction of seperate 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. 4b. We see that the Wannier functions 1-3 (on the triangular lattice) do not scale like the moiré length $\lambda_{\theta}=1/\theta$ near the magic angle. This is because for $\theta\leq 1^{\circ}$, atomic relaxations cause the formation of $AA$ stacking regions that are fixed in size even as the angle decreases Dai et al. (2016); Nam and Koshino (2017); Zhang and Tadmor (2018); Carr et al. (2018b); Yoo et al. (2019). 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 Zhang and Tadmor (2018). 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 ($\theta<0.85^{\circ}$) we see Wannier functions #1 and #2 ($AA_{\pm}$) begin to grow in size, due to increasing hybridization with electronic structures not associated with the $AA$ stacking spots.

III.3 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 $\bm{k}\cdot\bm{p}$ model by adding an onsite energy difference $\pm\delta_{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 $\bm{k}\cdot\bm{p}$ model by including a $\pm\delta_{S}$ for the $A$ and $B$ sublattice degrees of freedom (determined by $\alpha$ or $\beta$ 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, $\delta_{S}$, and layer anti-symmetric sublattice energy, $\delta^{\prime}_{S}$. For example, the layer symmetric perturbation means the $A$ ($B$) sublattice of both layers feels a positive (negative) onsite perturbation, while the layer anti-symmetric 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 $\bm{k}\cdot\bm{p}$ Hamiltonian is given by $H^{\prime}_{k}=H_{k}^{0}+\delta_{i}\Delta H_{k}$. One can easily calculate the first order perturbation $\Delta H$ for the reduced Wannier models by taking the $n$-band subspace projector $P_{n}$ and sandwiching it around the perturbation of $H^{\prime}_{k}$: $\Delta H=P_{n}^{\dagger}\Delta H_{k}P_{n}$. The $\bm{k}\cdot\bm{p}$ band structures and a summary of the largest terms of $\Delta 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 $\bm{k}\cdot\bm{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 8-band 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 $\delta_{S}$. In the reduced 8-band Hamiltonian, the layer-symmetric sublattice perturbation manifests as ideal onsite energies for Wannier functions $1$ and $2$ ($0.5\delta_{S}$ and $-0.5\delta_{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 non-degenerate band (in each valley) carries a nonzero Chern number. This is a direct consequence of the topological character of the flat bands in TBLG Po et al. (2018b). The anti-symmetric perturbation $\delta^{\prime}_{S}$ is not shown in Fig. 5, but it has a similar band structure to $\delta_{S}$ with a smaller gap-opening (roughly $30\%$ of that for $\delta_{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\delta^{\prime}_{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 moiré superlattice. The low-energy Hamiltonians can be truncated symmetrically, and they need hoppings with lengths only up to 3 moiré 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 moiré length $\lambda_{\theta}$ while those defined on the $AB$ or $DW$ geometry scale like $\lambda_{\theta}$, 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 symmetry-broken states and removes the unconventional long-range coupling of the propeller orbitals. Furthermore, with complete 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^{\circ},1.3^{\circ}]$. 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 not . These codes include the 8-band and 5-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.

Acknowledgements.

HCP and AV thank L. Zou and T. Senthil for earlier related collaborations on TBLG. This work was supported by the STC Center for Integrated Quantum Materials, NSF Grant No. DMR-1231319 and by ARO MURI Award W911NF-14-0247. HCP was supported by a Pappalardo Fellowship at MIT and a Croucher Foundation Fellowship. AV was supported by a Simons Investigator award and NSF-DMR 1411343.