key: cord-0145508-y0vyvzxx
authors: Cea, Tommaso; Pantaleon, Pierre A.; Guinea, Francisco
title: Band structure of twisted bilayer graphene on hexagonal boron nitride
date: 2020-05-15
journal: nan
DOI: nan
sha: 9d93f49145fd3165497d3ce05983da0ee54185a4
doc_id: 145508
cord_uid: y0vyvzxx

The effect of an hexagonal boron nitride (hBN) layer close aligned with twisted bilayer graphene (TBG) is studied. At sufficiently low angles between twisted bilayer graphene and hBN, $theta_{hBN} lesssim 2^circ$, the graphene electronic structure is strongly disturbed. The width of the low energy peak in the density of states changes from $W sim 5 - 10$ meV for a decoupled system to $sim 20 - 30$ meV. Spikes in the density of states due to van Hove singularities are smoothed out. We find that for a realistic combination of the twist angle in the TBG and the twist angle between the hBN and the graphene layer the system can be described using a single moir'e unit cell.

The effect of an hexagonal boron nitride (hBN) layer close aligned with twisted bilayer graphene (TBG) is studied. At sufficiently low angles between twisted bilayer graphene and hBN, θ hBN 2 • , the graphene electronic structure is strongly disturbed. The width of the low energy peak in the density of states changes from W ∼ 5 − 10 meV for a decoupled system to ∼ 20 − 30 meV. Spikes in the density of states due to van Hove singularities are smoothed out. We find that for a realistic combination of the twist angle in the TBG and the twist angle between the hBN and the graphene layer the system can be described using a single moiré unit cell.

Introduction. The discovery of insulating behavior at integer filling and superconductivity in TBG [1-5] motivated a recent effort in the study of a wide class of van der Waals heterostructures, displaying moiré patterns on length scales much larger than the lattice constant of their constituent layers. The periodicity induced by the moiré can affect sensitively the electronic structure of the material, giving rise to narrow, almost dispersionless, flat bands. The kinetic quenching may favor the interactions between the electrons, paving the way for the appearance of strongly correlated phases.

So far, transport measurements on TBG have been performed with the sample either encapsulated between two hBN clapping layers, see e.g. [1, 2, 6-10], or suspended on a substrate of hBN [11, 12] . Recently, it has been observed [13] that the presence of an additional layer of WS 2 between hBN and TBG stabilizes the SC phase in a range of angles wider than reported previously. The presence of hBN breaks the in-plane two-fold rotational symmetry (C 2 ), gapping out the Dirac crossing of monolayer graphene [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] . Placing TBG on top of hBN accounts for two coexisting moiré patterns: that induced by the mismatch between the lattice constants of hBN and graphene [27] [28] [29] [30] [31] , and that induced by the relative orientation between the two graphene layers of the TBG. If the hBN is aligned with the adjacent graphene layer and the relative twist between the two layers of graphene is close to 1 • , then the two moires have very similar periods, ∼ 13nm, even if they are not commensurate. This sensitively affects the band structure of TBG close to charge neutrality (CN). To a first approximation, it is expected that the breaking of inversion symmetry induced by the hBN layer gives rise to a gap at the Dirac point of the TBG, separating two flat conduction and valence bands, carrying opposite Chern numbers, C = ±1 [32, 33] . This analysis may explain the observed anomalous Hall effect at the integer band filling, ν = 3 [8, 34, 35] . The existence of flat bands can lead to Chern insulators, with features similar to those found in the Quantum Hall Effect [36, 37] . * These authors contributed equally In the following, we assume that the twist angle, θ T BG , in the TBG is fixed to a value near a magic angle, and study the effect of a hBN layer as function of the angle between this layer and the neighboring graphene layer, θ hBN . The periodicities of the two moiré patterns described above are L T BG ≈ d G /θ T BG and L G,hBN ≈ d G / θ 2 hBN + (d hBN /d G ) 2 , where d G and d hBN are the lattice constants of graphene and hBN. The overall structure resembles the arrangement in a twisted graphene trilayer, or in twisted graphite [38] [39] [40] . The two moiré lattices define a generically incommensurate superstructure.

Interestingly, a realistic choice of parameters allows us to define a single moiré unit cell for the whole system. This possibility permits an accurate study of the electronic properties. Interaction effects are included using the unrestricted Hartree Fock approximation [41] . We analyze the similarities and differences with the electronic structure of TBG decoupled from the substrate, and also with other graphitic systems which show narrow bands [42] .

The model. A sketch of the atomic arrangement to be considered, and of its Brillouin zone is shown in Fig. [1] . We assume, as in the continuum model for TBG [43, 44] that one side of the hexagonal Brillouin zone of the superlattice connects the corners of the Brillouin zones of each pair of layers, see 

where n x and n y are unit vectors along the x and y axes, and we have expanded the exact expressions to lowest order, θ T BG , θ hBN , d hBN /d G − 1 1. In order for the two moires to have the same unit cell, we need the vectors K G,top −K G,bottom and K G,top −K hBN to have the same length, and to make an angle equal to (2π)/3. These two conditions imply that:

For a fixed value of d hBN /d G these equations are satisfied when

For d G = 2.46Å, d hBN = 2.50Å, and d hBN /d G − 1 ≈ 0.017 we obtain θ T BG 1.05 • . This number is reasonably close to the twist angles where TBG shows a non trivial phase diagram. The twist of hBN, θ hBN ≈ 0.52 • is close to perfect alignment. The presence of a unique moiré pattern in hBN/TBG heterostructures is consistent with recent scanning-tunneling-microscopy (STM) maps [12] , where, in some samples, only the moiré pattern identified by the TBG appears.

We model the Hamiltonian of the TBG within the low energy continuum model, see [43] [44] [45] [46] . The effect of the hBN is included by means of an effective periodic potential acting on the nearest graphene layer [31, 47] . In what follows we refer to the parametrization of such potential as given in the Ref. [48] . A detailed description of the model is given in [49] .

We first study the arrangement described by the angles θ T BG = 1.05 • and θ hBN = 0.52 • , where a single moiré unit cell describes the system, as shown in Fig. [1] . Note that three different stacking configurations of the twisted system can be identified, see [49] , Figs. [S3] and [S4]. The three stackings differ in the relative arrangement of layers which are second nearest neighbors. As shown below, the three geometries lead to different electronic structures.

The band structure and density of the states (DOS) per unit cell of the hBN/TBG are shown in the Fig. [2] (black lines), and compared to that of the TBG (red dashed lines). A c = √ 3L 2 T BG /2 is the area of the moiré unit cell. The presence of the substrate strongly affects the spectrum of the TBG. The flat bands at CN of the TBG become dispersive in the hBN/TBG stack, acquiring a finite bandwidth of ∼ 6 − 8 meV, which is almost twice that of the TBG. As a consequence, the peak in the DOS of the TBG at CN is strongly smoothed and also split in the hBN/TBG, giving rise to an insulating structure with a small band gap, which is due to the breaking of C 2 induced by the hBN layer. The effect of a self consistent Hartree potential away from CN is shown in Fig. [3] and [4] . The Chern numbers of these bands are shown in [49] . It is worth mentioning that Chern numbers of up to C = 3 are obtained. The effect of the exchange term at the neutrality point is analysed in [49] . As in the absence of a substrate [50, 51] , the bands are significantly distorted by the Hartree potential.

Generic values of the twist angle between hBN and TBG cannot be described by a simple moiré unit cell. This is, for example, the case for perfect alignment: θ hBN = 0 • . Because of the lack of commensuration, we cannot define a crystal momentum. In order to study the energy spectrum, we project the perturbation induced by the hBN on the low energy states of the TBG and solve a dual lattice in the reciprocal space of the TBG, as detailed in [49] . The scheme follows closely continuum models for a graphene monolayer on hBN, where the perturbation due to the hBN layer is projected onto the graphene Dirac cone. An infinite number of minibands emerge, induced by the periodicity of the poten- tial due to the hBN layer. In a similar manner, the TBG bands are replicated and coupled in our calculation. A similar method has been recently used in the Ref. [52] for studying the quasi-crystalline electronic spectrum of the non-commensurate 30 • -TBG.

The quasi-band structures, obtained by varying the momentum k in the BZ of the TBG, and the DOS are shown in the Fig. [5] for different orientations between the hBN and graphene, and θ T BG = 1.05 • . The black and red hexagons show the two different BZs of the TBG and of the hBN/G, respectively. The red lines refer to the band structure and the DOS of the unperturbed TBG, that also includes the long-wavelength staggered potential induced by the hBN, weakly breaking the C 2 symmetry of graphene (see [49] ). As is evident, the hBN strongly affects the spectrum close to CN at small angles, θ hBN < 1 • , where the moiré identified by the TBG and that identified by the substrate of hBN have similar periods. In contrast to the narrow bands of the TBG, the hBN/TBG exhibits a broad structure, with a bandwidth of approximatively 30meV, which can even overlap with the higher energy bands. This is for example the case of θ hBN = 0.52 • , Fig. [5](b), which is close to commensuration. In addition, the band broadening at CN lowers the DOS of the hBN/TBG as compared to the sharp van Hove singularity of the TBG, that has been cut out of the energy scale in the central panels of Fig. [5] . In general, the high energy spectrum is barely affected by the hBN. Upon increasing θ hBN , the bandwidth at CN gradually shrinks, while the DOS gains intensity and the central bands further separate from the rest of the spectrum. At θ hBN = 2 • the effect of the hBN is almost completely negligible and we recover the narrow band feature of the TBG, except for a constant gap. Even though we are using a small value of the staggered potential, ∆ = 3.62meV, as given by the Ref. [48] , we checked that our results remain valid in a wide range of values of ∆. The details are shown in [49] , where we report the case of ∆ = 40meV.

Finally, we show plots of the changes in the charge density distribution induced by the substrate in [49] . The sixfold symmetry of the TBG is reduced to threefold, in the case where the two moirs coincide. We expect a lower symmetry in the general case. These results are in agreement with the reduced symmetry observed in STM experiments [53] .

Conclusions. We have analyzed the effect of a nearly aligned hBN substrate on the low energy bands of a twisted graphene bilayer. We find that at large enough angles between the TBG and the next hBN layer, θ hBN 2 • , the TBG is effectively decoupled from the substrate.

Only those effects which do not average to zero over the hBN/G moiré unit cell, such as a finite gap due to the lack of inversion symmetry, survive.

On the other hand, perturbations of zero average over the unit cell change significantly the electronic structure. The narrow peak in the DOS associated to the narrow bands of the TBG broadens appreciably. The van Hove singularities are smoothed out, and the width of the peak near the Dirac energy increases from W ∼ 5−10 meV for a decoupled TBG to W 20 − 30 meV for a well aligned hBN/TBG stack. This peak overlaps with higher energy bands. In addition, the gap expected from the lack of inversion symmetry becomes filled with states which arise from the non uniform part of the perturbation due to the substrate.

Acknowledgements. This work was supported by funding from the European Commision, under the Graphene Flagship, Core 3, grant number 881603, and by the grants NMAT2D (Comunidad de Madrid, Spain), SprQuMat and SEV-2016-0686, (Ministerio de Ciencia e Innovacin, Spain).

This paper, main ideas, theory and figures have being developed during the COVID-19 lockdown. Supplementary information for Band structure of twisted bilayer graphene on hexagonal boron nitride

We consider the heterostructure obtained by placing the TBG on top of a substrate with an hexagonal atomic structure and lattice constant d s . Let d G,1 = d G (1, 0) , d G,2 = d G 1/2, √ 3/2 be the primitive vectors of the Bravais lattice of the two unrotated layers of graphene, with d G = 2.46Å, and d s,1 = d s (1, 0) , d s,2 = d s 1/2, √ 3/2 those of the substrate. We assume that at the origin, O = (0, 0), the atomic positions of the three Bravais lattices coincide. Then, given i, l integer numbers, we consider the atomic positions:

of the bottom graphene layer,

of the top graphene layer, and . This rotation moves the atoms C b and C t to the same position, lying on the vertical axis. The period of the moiré is:

Keeping the substrate fixed with respect to the aligned configuration, the moiré superlattice of the TBG is preserved by the presence of the substrate if:

This condition indeed implies that the rotation of the two graphene layers moves the atoms C b and C t to the position occupied by S. The rhs of the Eq. (S2) expresses d s as a function of two integer numbers, i and l, setting a sufficient condition for the moiré of the TBG to persist in the presence of the substrate. The relative twist between the substrate and each of the two graphene layers is: θ s = θ/2. In particular, for i = l = 31, the Eq. (S2) gives: d s 2.50Å, which is actually the value of the lattice constant of hBN. The corresponding angle is: θ 1.05 • , with L 13.4nm.

We describe the TBG within the low energy continuum model considered in Refs. [43] [44] [45] [46] , which is meaningful for sufficiently small angles, so that an approximatively commensurate structure can be defined for any twist. The moiré mini-BZ, resulting from the folding of the two BZs of each monolayer (see Fig.S2(a) ), is generated by the two reciprocal lattice vectors:

shown in green in Fig. S2(b) . For small angles of rotation the coupling between different valleys of the two monolayers can be safely neglected, as the interlayer hopping has a long wavelength modulation. In what follows we describe the model for the behavior of the two K-valleys of the twisted monolayers, where K = 4π(1, 0)/(3d G ) is the Dirac point of the unrotated monolayer graphene. The case corresponding to the opposite valleys, at K = −K, directly follows from time reversal symmetry, by inverting: k → −k.

The Hamiltonian of the TBG is a 4 × 4 matrix, with entries: (A b , B b , A t , B t ) , where A, B denote the sub-lattice and b, t refer to the bottom and top layer, respectively. Without loss of generality, we assume that in the aligned S1 . Sketch of the aligned Bravais lattices of graphene (red points) and the substrate (black points). θ is the angle identifying the TBG. Keeping the substrate fixed with respect to the aligned configuration and rotating the bottom (top) graphene layer by −θ/2(+θ/2) around O, moves the atoms C b and Ct to the position occupied by S, preserving the moiré superlattice of the TBG. configuration, at θ = 0, the two layers are AA-stacked. In the continuum limit, the effective Hamiltonian of the TBG can be generally written as [43] [44] [45] [46] :

where:

is the Dirac Hamiltonian for the layer l = b, t, v F = √ 3td G /(2 ) is the Fermi velocity, t is the hopping amplitude between localized p z orbitals at nearest neighbors carbon atoms, θ b,t = ∓θ/2, τ θ l = e iτzθ l /2 (τ x , τ y ) e −iτzθ l /2 , τ i are the Pauli matrices, and K l = 4π (cos θ l , sin θ l ) /(3d G ) are the Dirac points of the two twisted monolayers, which identify the corners of the mini-BZ shown in Fig. S2. U (r) is the interlayer potential, which is a periodic function in the moiré unit cell. In the limit of small angles, its leading harmonic expansion is determined by only three reciprocal lattice vectors [43] :

where the amplitudes U (G) are given by:

In the following we adopt the parametrization of the TBG given in the Ref. [46] : v F /d G = 2.1354eV, g 1 = 0.0797eV and g 2 = 0.0975eV. The difference between g 1 and g 2 , as described in [46] , accounts for the inhomogeneous interlayer distance, which is minimum in the AB/BA regions and maximum in the AA ones, or it can be seen as a 

where m is the band index, k is the momentum in the mini-BZ, G = n 1 G 1 + n 2 G 2 are reciprocal lattice vectors, with n 1 , n 2 integers, a is the sub-lattice/layer index and φ m,k,a (G) are numerical eigenvectors. Upon varying k, the eigenvalue E m (k) corresponding to |m, k defines the m-th Bloch band. Note that the φ's can be chosen within a gauge degree of freedom. Given the gauge U, the mapping between equivalent points of the reciprocal space is defined according to:

where U k (G 0 ) U † k (G 0 ) = 1. The Eq. (S8a) in turn implies:

In particular, U k (G 0 ) = 1 for a periodic gauge.

In the numerical calculations, the number of Fourier components defining the eigenfunctions |m, k is bounded by a cutoff: |G| < G c , where G c is chosen in order to achieve the convergence of the low energy bands.

The commensurate heterostructure hBN/TBG accounts for three possible non-equivalent configurations of the moiré unit cell, as shown in Fig. S3 , where the green, light blue and red hexagons refer to the hBN, bottom and top graphene layer, respectively, and C denotes the center of the hexagon. Each configuration is determined uniquely by the local stacking between the hBN and the TBG in the center of the cell: AAA, BAA and CAA. Consequently, there are nine possible stacking arrangements, three for each configuration. They are shown schematically in Fig. S4 , where the black points represent the carbon atoms, the points at the right (left) edges of the horizontal lines represent atoms of type A (B), and we assume that the sub-lattice A of the hBN is occupied by nitrogen atoms, without loss of generality. The case of the heterostructure hBN/TBG/hBN is more complex, accounting for nine non-equivalent configurations, which are determined by the local stacking between the bottom and top hBN and the TBG in the center of the cell. We describe the effect of hBN on the nearest graphene layer by means of an effective potential, periodic in the moiré unit cell [47] :

where α = A, B, C labels the three configurations of the Fig. S3(a),(b) ,(c), respectively, w α 0 and ∆ α represent a spatially uniform scalar and mass term (note that hBN breaks inversion symmetry, and allows for a mass term [14]),

identify the first star of reciprocal lattice vectors and we assume that the modulation of V α SL at smaller wavelengths is negligible. The amplitudes v α SL (G j ) are given by:

where are position-dependent mass and gauge terms, respectively. We use the parametrization of V α SL given by the Ref. [48] . In particular, the set of parameters for the configuration: α = A, is: The configurations: α = B, C, are related to α = A by a rotation of ±2π/3 in the parameters space, as detailed in Ref. [48] . The continuum Hamiltonians of the hBN/TBG and of the hBN/TBG/hBN are given by:

respectively.

In this section, we address the topological phases of the bands for the three possible non-equivalent configurations, where the breaking of either time reversal or inversion symmetry in the single valley model in Eq. S12 allows for a finite Berry curvature Ω k,l = 2 Im ∂ kx ψ k,l |∂ ky ψ k,l ,

where l is a band index with energy E l (k) and wavefunctions ψ k,l . Notice that due to time-reversal symmetry, the Berry curvature in each graphene valley has opposite sign and hence the total Chern number is zero. However, by assuming an absence of intervalley scattering, the topological invariants can be defined separately [55, 56] . For the different stacking configurations, the bands are isolated and their Berry curvature is well defined. Therefore, we can assign a valley Chern number C l to the band l which is given by the integral of the Berry curvature about the moiré Brillouin zone:

We use the algorithm in Ref.

[57] to compute the Berry curvature and valley Chern number. The presence of an hBN substrate gives rise to an insulating structure with a small band gap, which is due to the breaking of inversion symmetry due to the hBN layer. As is shown in Fig. 2 in the main text, the resulting isolated conduction and valence bands in each valley and in each stack configuration carry Chern numbers C = ±1, in agreement with Ref. [32, 33] . However, similar to other graphene heterostructures [42] and consistent with experimental observations [6, 34], we find different topological phases induced by interaction effects. This is summarized in Table. I where we display the Chern number as a function of the filling of the conduction band, for the different stack configurations. 

In Fig. S5 we show the real space charge density computed at the points Γ and K of the BZ, and obtained for pristine TBG, panels a) and b), and for TBG with a substrate of hBN, panels c) and d). Interestingly, we find that the C 6 symmetry of the charge density of pristine TBG is lowered to C 3 by the presence of the substrate, an effect that is more evident at the K point.

Using the formalism detailed in [41] , we performed a fully Hartree-Fock calculation of the band structure of the commensurate hBN/TBG heterostructure, in the stacking configuration α = A. The results are shown in the Fig. S6 for different fillings of the conduction band, ν. As it is evident, the spectra display pinning of the Fermi level, E F , at the van Hove singularity at finite filling and are quite similar to the ones shown in the main text, which include only the Hartree potential. Consequently, we argue that the Fock contribution is negligible for this kind of system, in contrast to the case of freely standing TBG, where the Fock terms have been shown to change considerably the band structure and to open spectral gaps at integer filings [41, 58, 59] .

Here we focus on the general case in which the moiré pattern identified by the hBN and its nearest graphene layer and that of the TBG are not commensurate, making the system non-periodic. This happens, for example, when the hBN is perfectly aligned with its nearest graphene layer.

We compute the quasi-band structure and the DOS, shown in the Fig. 5 of the main text, by projecting the perturbation induced by the hBN on the low energy eigenstates of the TBG.

Let us write the local potential induced by the hBN on the nearest graphene layer as:

where V 0 = w 0 τ 0 + ∆τ z is the uniform contribution, while δV SL (r) = 5 j=0 v SL G j e iGj ·r is the contribution at finite wavelengths. In the following we choose the parametrization of V SL given by the Eq. (S11), which corresponds to set the origin in a region of AAA stacking. TheG j are reciprocal vectors of the moiré superlattice identified by the hBN/G, as given by: that of the freely standing TBG. To check the robustness of our findings against larger values of ∆, we computed the spectrum of the non-commensurate hBN/TBG for ∆ = 40meV. The results are shown in the Fig. S8 , for the same values of θ hBN considered in the main text. The red lines refer to the unperturbed spectrum, which indeed displays a sizeable gap at CN. As is evident, for small angles the perturbed spectrum remains gapless, in very good agreement with the Figs. 5(a)-(b) of the main text. This means that the finite wavelengths term, δV SL (r), gives the leading contribution, rather than ∆. This contribution becomes however negligible upon increasing θ hBN , and we recover the spectral gap already for θ hBN 1 • . 

cones in graphene/hexagonal boron nitride

Dynamic band-structure tuning of graphene moirésuperlattices with pressure

Evendenominator fractional quantum hall states at an isospin transition in monolayer graphene

Accurate Gap Determination in Monolayer and Bilayer Graphene/ h-BN Moiré Superlattices

Scanning tunnelling microscopy and spectroscopy of ultra-flat graphene on hexagonal boron nitride

Emergence of superlattice Dirac points in graphene on hexagonal boron nitride

Commensurate-incommensurate transition in graphene on hexagonal boron nitride

Electronic properties of graphene/hexagonal-boron-nitride moiré superlattice

Spontaneous strains and gap in graphene on boron nitride

Twisted Bilayer Graphene Aligned with Hexagonal Boron Nitride: Anomalous Hall Effect and a Lattice Model

Mechanism for anomalous hall ferromagnetism in twisted bilayer graphene

Intrinsic quantized anomalous Hall effect in a moir\'e heterostructure

Reversible modifications of linear dispersion: Graphene between boron nitride monolayers

Ground State and Hidden Symmetry of Magic Angle Graphene at Even Integer Filling

Charged Skyrmions and Topological Origin of Superconductivity in Magic Angle Graphene

Electronic spectral properties of incommensurate twisted trilayer graphene

Flatbands and perfect metal in trilayer moiré graphene

Twists and the Electronic Structure of Graphitic Materials

Band structure and insulating states driven by the coulomb interaction in twisted bilayer graphene

Narrow bands and electrostatic interactions in graphene stacks

Graphene bilayer with a twist: Electronic structure

Moiré bands in twisted double-layer graphene

Continuum model of the twisted graphene bilayer

Maximally Localized Wannier Orbitals and the Extended Hubbard Model for Twisted Bilayer Graphene

Generic miniband structure of graphene on a hexagonal substrate

Moiré band model and band gaps of graphene on hexagonal boron nitride

See Supplementary Information

Electrostatic effects, band distortions, and superconductivity in twisted graphene bilayers

Electronic band structure and pinning of Fermi energy to Van Hove singularities in twisted bilayer graphene: A selfconsistent approach

Quasicrystalline electronic states in 30 • rotated twisted bilayer graphene

(S19)where we used the Eq.s (S8). Thus, the Eq. (S17) finally becomes:which defines a dual tight binding multi-orbital model in the reciprocal space, where the k points act as sites, E m (k) are the onsite energies and only the overlaps between nearest neighbor sites in the triangular lattice are allowed. The corresponding hopping integrals are:Because the two moires are not commensurate, the above overlaps span ergodically the BZ of the TBG. Note that this dual tight binding framework is similar to that considered in the Ref. [52] for studying the quasi-crystalline electronic structure of the 30 • -TBG. In our calculations, for each k we consider the reduced dual model obtained by including only the 19 sites lying in the first two stars ofG vectors surrounding k, accounting for 42 overlaps. This is schematically shown in the Fig.  S7(a) , where the blue points represent the sites, the blue lines the overlaps between them, the black hexagons are the equivalent BZs of the TBG and θ hBN = 0 • . The Fig. S7(b) shows the projection of the overlaps in the first BZ of the TBG. Note that, in the commensurate case, all the neighboring sites would collapse into the same point, k, upon projection. We checked that this approximation is sufficient to achieve the convergence of the DOS. Within this framework, the spectrum at the wave vector k is obtained by diagonalizing a matrix of size (19N b where N b is the number of bands of the TBG taken into account in the projected Hamiltonian of the Eq. (S20). We consider the N b = 14 bands closest to the CN point. The multiple sets of bands shown in the left panels of the Fig. 5 of the main text arise from the diagonalization of the dual tight binding model of the Eq. (S20). Including more dual sites, there would appear many additional bands which, however, would be just replicas of the ones already shown, but shifted by a different origin in the BZ. As a consequence, they would not change appreciably the DOS and the other physical quantities, as discussed in the Ref. [52] .The parametrization of the hBN induced potential given by the Ref. [48] , Eq. (S11), accounts for a small value of the staggered potential: ∆ = 3.62meV. This is the reason why the unperturbed band structure, obtained by neglecting δV SL (r) and shown by the red lines in the Fig. 5 of the main text, is almost gapless at CN and strongly resembles