Exchange interaction, disorder, and stacking faults in rhombohedral graphene multilayers
EExchange interaction, disorder and stacking faults in rhombohedral graphenemultilayers
James H. Muten, Alex J. Copeland, and Edward McCann ∗ Physics Department, Lancaster University, Lancaster, LA1 4YB, UK
We generalize the mean-field Hartree Fock theory of gapped electronic states at charge neutral-ity in bilayer graphene to thin films of rhombohedral graphite with up to thirty layers. For theground state, which has an odd spatial parity within each spin-valley flavor and an antiferromag-netic arrangement of four flavors, the order parameter (the separation of bands at the valley center)saturates to a constant non-zero value as the layer number increases, whereas the band gap decreaseswith layer number. We take into account chiral symmetry breaking disorder in the form of randomlayer potentials and chiral preserving disorder in the form of random values of the interlayer cou-pling. The former reduces the magnitude of the mean band gap whereas the latter has a negligibleeffect, which is due to self-averaging within a film with a large number of layers. We determine theground state in the presence of an individual stacking fault embedded within a film of rhombohedralgraphite. For a Bernal stacking fault, the ground state is also an odd parity antiferromagnetic state,and the fault can be interpreted as introducing a small coupling between two independent sectionsof rhombohedral graphite. For a twin boundary stacking fault, however, the ground state is aneven parity antiferromagnetic state, and the fault introduces stronger coupling across the system.In the presence of stacking faults, each individual rhombohedral section with m layers contributesa pair of low-energy flat bands producing a peak in the Berry curvature located at a characteristic m -dependent wave vector. For both types of stacking fault, the Chern number per spin-valley flavorfor the filled valence bands in the ground state is equal in magnitude to the total number of layersdivided by two, the same value as for pristine rhombohedral graphite. I. INTRODUCTION
Recently, topological flat bands have been the sub-ject of intense research in twisted bilayer graphene [1, 2]as well as other two-dimensional systems including theLieb, honeycomb and kagome lattices [3–11]. Thereare also topological flat bands in rhombohedral multi-layer graphene (RMG) [12–18] in which alternating intra-and interlayer coupling act like staggered hopping in theSu-Schrieffer-Heeger (SSH) model [19–22], as shown inFig. 1. Now there is fresh interest in RMG due to progressin fabricating and characterizing samples with a largelayer number [23–31] culminating in the realization ofhigh-quality films with up to fifty layers [32].In high mobility samples, at charge neutrality, low tem-perature and for zero external fields, low energy bandshave been observed to be gapped in bilayer graphene [33–38], Bernal multilayers with up to N = 8 layers [39–41], rhombohedral multilayers with up to N = 4 lay-ers [42–44] and recently with N ≈
12 [32]. A num-ber of different interaction-induced broken symmetrystates have been proposed for bilayer graphene [17, 45–59] including pseudospin layer antiferromagnetic (LAF)states [17, 45, 46, 51, 55] in which electrons with differ-ent spin and valley flavors spontaneously accumulate ondifferent layers, creating an odd parity state that breaksinversion symmetry and opens a gap. Owing to the anti-ferromagnetic configuration of four flavors, there is no netcharge accumulation on summing over them and, thus, ∗ [email protected] RhombohedralgrapheneA B A B A B A B A B A B A B BernalfaultA B A B A B A B A B A B A B Twin boundaryfaultA B A B A B A B A B A B A B FIG. 1. Schematic side view of the lattice of rhombohedralgraphene with N = 7 layers showing the pristine lattice, alattice with a single Bernal stacking fault, and a lattice witha single twin boundary stacking fault. Labels indicate the A n , B n atomic sites on the n th layer, horizontal solid linesindicate intralayer hopping with parameter γ , vertical solidlines indicate interlayer hopping γ . no cost in terms of Hartree energy. The evolution ofsimilar gapped states with layer number has been dis-cussed for both Bernal [60, 61] and rhombohedral multi-layers [17, 62–64].In this paper, we generalize the Hartree Fock mean-field theory [45, 46, 60–62] of the pseudospin LAF state toRMG with a large number of layers (up to N = 30) in or-der to determine the layer dependence of the interaction-induced band gap at charge neutrality. We find thatthe strength of the interacting state, as characterized bythe separation of the bands exactly at the valley center,saturates with layer number, but that the actual bandgap decreases with layer number; this agrees with den-sity functional theory (DFT) which predicted the bandgap for RMG with up to eight layers [64].We then consider the robustness of the LAF state a r X i v : . [ c ond - m a t . m e s - h a ll ] F e b to defects by including values of tight-binding parame-ters that are constant within each layer (thus preserv-ing translational invariance within each layer) but thatvary randomly between layers. We compare disorder thatpreserves chiral symmetry (random values of the inter-layer hoppings) with disorder that breaks chiral symme-try (random layer potentials). For weak disorder, we findthat the mean band gap is diminished by chiral breakingdisorder, but it is almost insensitive to chiral preservingdisorder; this is similar to the behavior of gapless edgestates in the non-interacting SSH model [65–69].Another type of defect is a localized stacking fault [32,70, 71] within a large RMG system, namely a Bernalfault or a twin boundary fault, Fig. 1. They are partic-ularly interesting because they introduce additional flatbands into the energy spectrum. For a Bernal fault, wefind that it introduces a weak connection between twosections of RMG and that the interacting ground stateis a straightforward generalization of the LAF state withodd parity. The twin boundary fault, however, createsa stronger connection between two RMG sections: theinteracting ground state is also a LAF state, but witheven parity within each spin-valley flavor. For both ofthese types of ground state, the Chern number per spin-valley flavor has magnitude N/
2, as is the case for pristineRMG [20, 72, 73].Section II describes the methodology including thenon-interacting Hamiltonian and the Hartree Fock mean-field theory. We use the minimal model with nearest-neighbor intralayer and interlayer hopping parameterizedby γ and γ , respectively, but neglecting other tight-binding parameters. This is done for simplicity and, inparticular, it dramatically simplifies the calculation of theexchange interaction allowing us to consider large layernumber N (cid:29)
1. Section III describes the LAF statein pristine RMG. We introduce a toy two-band modelthat can be solved analytically to give simple expressionsfor the parameter dependence of the band gap that arebroadly in qualitative agreement with the full numeri-cal model. With the full numerical model, we determinethe layer and temperature dependence of the band gap.Then, our main results are described in Section IV fordisorder, Section V for the Bernal stacking fault, andSection VI for the twin boundary fault.
II. METHODOLOGYA. Effective mass model
The lattice of RMG with N layers consists of two in-equivalent sites A n , B n on each layer, n = 1 , , . . . N ,with sites B n located below A n +1 , Fig. 1. In the tight-binding model, interlayer coupling between p z orbitals onthe B n and A n +1 sites hybridizes those orbitals leadingto gapped bulk conduction and valence bands. In thesurface layers, however, the A and B N sites don’t haveneighbors in the next layers so their p z orbitals aren’t hybridized by interlayer coupling, resulting in low-energysurface states within the bulk gap in the vicinity of eachof two valleys K ± .In a basis of p z orbitals on A , B , A , B , . . . , A N , B N sites, the non-interacting Hamiltonian of RMG with N layers [12, 15, 16] may be written near each valley as H = D V · · · V † D V · · · V † D V · · · V † D · · · ... ... ... ... . . . , (1)where we use 2 × D = (cid:18) γ κ † γ κ (cid:19) , V = (cid:18) γ (cid:19) , κ = ξk x + ik y k c . Here k = ( k x , k y ) is the wave vector measured fromthe center of valley K ξ with valley index ξ = ±
1, and k c = γ / (¯ hv ). Block D describes intralayer nearest-neighbor hopping with velocity v = ( √ / aγ / ¯ h andin-plane lattice constant a , block V describes interlayerhopping with parameter γ between successive B n and A n +1 sites. For numerical diagonalization of (1) we use γ = 3 .
16 eV, γ = 0 .
381 eV [74], and a = 2 . B. Mean-field theory
Electron-electron interactions are included within amean-field Hartree-Fock approximation [45, 46, 62] and,in particular, we follow the methodology applied toBernal-stacked multilayer graphene in Refs. 60 and 61.The total Hamiltonian is ˆ H tot = ˆ H + ˆ V MF whereˆ H = (cid:88) k σXX (cid:48) H k XX (cid:48) c † k σX c k σX (cid:48) , (2)ˆ V MF = (cid:88) k σXX (cid:48) (cid:104) U (H) X δ XX (cid:48) − W k σXX (cid:48) (cid:105) c † k σX c k σX (cid:48) . (3)Here X = A , B , A , B , . . . indexes the sublattices, σ =1 , , , ↑ , ↓ ) and valley( K + , K − ) degrees of freedom, and c † k σX and c k σX arecreation and annihilation operators, respectively. Thenon-interacting term ˆ H contains H k XX (cid:48) which is a matrixelement of (1); the interaction term ˆ V MF consists of theHartree U (H) and exchange W potentials, U (H) X = lim q → (cid:88) X (cid:48) V ( q ; z X − z X (cid:48) ) n X (cid:48) , (4) W k σXX (cid:48) = 1 L (cid:88) k (cid:48) V ( k − k (cid:48) ; z X − z X (cid:48) ) (cid:104) c † k (cid:48) σX (cid:48) c k (cid:48) σX (cid:105) , (5)where n X = (1 /L ) (cid:80) k σ (cid:104) c † k σX c k σX (cid:105) − n , L is thesystem area, z X is the vertical coordinate of sublat-tice X , and V ( q ; z ) = (2 πe / [ (cid:15) r q ]) exp( − q | z | ) is thetwo-dimensional Fourier transform of the Coulomb po-tential [75], (cid:15) r is the dielectric constant and we use d = 3 . n represents the background density of positive charge and,in the charge neutral case considered here, it is deter-mined by (cid:80) X n X = 0. In this case, the Hartree term (4)simplifies as U (H) X = − (2 πe /(cid:15) r ) (cid:80) X (cid:48) | z X − z X (cid:48) | n X (cid:48) .The strength of the electronic interactions is character-ized by the effective fine structure constant for graphene, α g = e (cid:15) r ¯ hv . (6)For example, for (cid:15) r = 2 then α g ∼
1. However, theHartree-Fock approximation tends to overestimate thestrength of the exchange interaction by neglecting screen-ing effects so we treat α g as a fitting parameter in therange 0 < α g ≤ . U (H) and W intoaccount. The summations over k are performed withina circle k < k (cid:63) around the K point with a cutoff k (cid:63) ,and we choose ¯ hvk (cid:63) ≈ k (cid:63) ,the ground state is determined for different values of thesystem size L ( i.e. different densities of k points), thenthe band gap is evaluated by extrapolation to L → ∞ .The numerical precision of our results is high so that un-certainties are negligible, and error bars are only shownin Section IV when we study random disorder. Never-theless, there are many sources of systematic uncertaintyincluding the choice of cutoff k (cid:63) , the omission of tight-binding parameters in the minimal model, and the valueof the interaction parameter α g . As described in Sec-tion III B, we find very close agreement of band gap val-ues at α g = 0 . N = 3 to N = 8 layers, and so we use α g = 0 . σ = 1 , , , K + , ↑ ) , ( K + , ↓ ) , ( K − , ↑ ) , ( K − , ↓ ). We neglect interactions be-tween different valleys because they are described bythe Coulomb interaction V ( q ; z ) with large wave vec-tors q ≈ K + − K − and this approximation treats thefour flavors on an equal footing. Within a given fla-vor, the exchange potential breaks chiral symmetry andcharge density is transferred between layers. In addi-tion, the four flavors have a certain relative configura-tion. By beginning the iterative procedure with differentinitial exchange profiles, we find different self-consistentsolutions, and we evaluate the total energy of each so-lution in order to determine the ground state. For ex-ample, for the LAF state, exchange for a given flavorhas odd parity with respect to spatial inversion, and thefour flavors are arranged in an antiferromagnetic con-figuration so that there is net charge polarization and,hence, no cost in terms of Hartree energy. Thus thisstate has lower energy at charge neutrality than, say, aferrimagnetic or ferromagnetic configuration of the fourflavors. Since our approximation treats the four flavors equally, it is unable to differentiate three distinct com-binations [17, 45, 46] of spins and valleys within the an-tiferromagnetic configuration: layer-antiferromagnetic inwhich the polarization of flavors ( K + , ↑ ) , ( K − , ↑ ) is oppo-site to that of ( K + , ↓ ) , ( K − , ↓ ), quantum anomalous Hallwhen ( K + , ↑ ) , ( K + , ↓ ) are opposite to ( K − , ↑ ) , ( K − , ↓ ),or quantum spin Hall when ( K + , ↑ )( K − , ↓ ) are oppositeto ( K + , ↓ ) , ( K − , ↑ ).Using the minimal model, the energy spectrum isisotropic around each K point and the eigenstates ofthe non-interacting Hamiltonian (1) at an arbitrary anglemay be related to those at a specific angle by a stacking-dependent unitary transformation [45, 76]. We assumethat the eigenstates of the interacting mean-field the-ory (2,3) also satisfy this rotational transformation [60],allowing for the k summations to be performed in onlyone specific direction with the exchange interaction (5)being determined via an integration with respect to thepolar angle of wave vector k (cid:48) . This simplification dra-matically reduces the numerical cost of the calculationsallowing for a study of multilayers with a large numberof layers.The non-interacting Hamiltonian (1) obeys chiral sym-metry [20]: matrix elements only connect A and B sites(not A to A or B to B ) and chiral symmetry ensuresparticle-hole symmetry of the electronic spectrum. Whenthe spectrum is gapped due to interactions, which gen-erally break chiral symmetry, we use the wave functionsin k space to determine the non-Abelian Berry curva-ture [72] for the occupied valence bands. This is thenintegrated in the vicinity of a single valley in order todetermine the Chern number per spin-valley flavor [73]. III. RHOMBOHEDRAL GRAPHENE
The band structure of RMG with N = 16 layers for asingle spin-valley flavor is shown in Fig. 2(a,b) for non-interacting and interacting electrons, respectively, as de-termined by numerical solving the mean-field theory (2,3)for 2 N bands, Eq. (1). The orbitals on the surface sites, A and B N , contribute to a pair of low-energy bandsthat are flat for k < ∼ k c where k c = γ / ¯ hv . The wave vec-tor k c corresponds to the point of the phase transitionbetween non-trivial and trivial topological phases in theSSH model [19–22]. For the interacting ground state, theband structure displays flavor degeneracy, Fig. 2(b). Theband gap, E g , is at finite wave vector k between the sur-face state bands, this is the difference in energy betweenthe conduction band minima and the valence band max-ima. We also consider the separation of the surface statebands at k = 0 which we refer to as the order parameter∆. In principle, it is possible to have non-zero ∆ evenwhen E g = 0 and, generally, ∆ ≥ E g .The exchange potential W k =0 σXX at the valley center k = 0 and for each site X (cid:48) = X for RMG with N = 8layers is shown in Fig. 2(c) [77]. For two flavors, the ex-change is negative on A sites and positive on all B sites, (a) (b) E ( e V ) E ( e V ) k/k c k/k c E g � (c) A1 B1 A2 B2 A3 B3 A4 B4 A5 B5 A6 B6 A7 B7 A8 B8
Site W ( m e V ) FIG. 2. Low-energy band structure of RMG with N = 16 lay-ers with (a) non-interacting electrons described by Hamilto-nian (1) and (b) interacting electrons described by the mean-field theory (2,3) with interaction strength α g = 0 .
3. Blacklines are bulk bands, blue lines are surface bands for k < ∼ k c where ¯ hvk c = γ . E g is the band gap (at k ≈ k c ), ∆ is the or-der parameter (the separation of the surface bands at k = 0).(c) The values of the exchange potential W k =0 σXX at the val-ley center k = 0 and for each site X (cid:48) = X for RMG with N = 8 layers [77]. The black solid line shows the exchangefor two flavors, the blue dashed line is exchange for the othertwo flavors. with large magnitude on the surface sites A and B N (solid black line); for the other two flavors, the exchangehas an inverted profile (dashed blue line). We refer to thisas being the odd antiferromagnetic state because the ex-change potential has odd parity within each flavor, andthe four flavors are arranged in an antiferromagnetic con-figuration.Before discussing the layer number N and temperature T dependence of E g and ∆ arising from the numericalcalculations, we consider a very simple two-band modelthat may be solved analytically in order to develop aqualitative picture of the LAF state. A. Two band model
We consider a mean-field Hamiltonian [78, 79] for eachof the four σ flavours with two sublattices A and B N ,and order parameter ∆ = 2 | w | due to the exchange in-teraction w which breaks sublattice symmetry: H = (cid:18) E + w − γ ( − κ † ) N − γ ( − κ ) N E − w (cid:19) , (7)where N ≥
2. Note that similar two-band models havebeen considered in similar contexts previously [45, 80]. The eigenvalues and eigenstates of H may be written as E s = E + s (cid:113) w + γ ( k/k c ) N , (8) ψ s = 1 (cid:112) E s − E )( E s − E + w ) (cid:18) E s − E + w − γ ( − κ ) N (cid:19) , (9)where s = ± k = | k | ≡ (cid:113) k x + k y . We will show this is a self-consistent solution under the approximation that we onlytake into account the contribution exactly at the valleycenter ( k = 0) in the exchange (5). This means that theoff-diagonal in sublattice exchange potential W k σA BN is zero because the summation over all k (cid:48) includes a fac-tor such as exp( iN φ (cid:48) ) arising from the chiral wave func-tions (9), where φ (cid:48) is the polar angle of the wave vector k (cid:48) .Considering the diagonal in sublattice exchange poten-tial W k σA A (5), then w = − π ¯ hvα g L (cid:88) s = ± (cid:88) k (cid:48) f ( E s ) | k (cid:48) | w E s − E ) , (10)where f ( E s ) = 1 / (exp[( E s − E ) / ( k B T )]+1) is the Fermi-Dirac distribution, E is the chemical potential, k B isBoltzmann’s constant and T is absolute temperature.Then, the equation for the order parameter ∆ = 2 | w | is ∆ = 2 γ (cid:20) α g g N (cid:18) ∆ k B T (cid:19)(cid:21) N/ ( N − , (11) g N ( x ) = (cid:90) ∞ dy (cid:112) y N sinh (cid:16) x (cid:112) y N (cid:17) (cid:16) x (cid:112) y N (cid:17) . (12)At zero temperature, g N ( ∞ ) is simply a number, and theorder parameter is explicitly given by∆( T = 0) = 2 γ (cid:18) α g g N ( ∞ )2 (cid:19) N/ ( N − , (13) g N ( ∞ ) = 1 √ π Γ (cid:18) N (cid:19) Γ (cid:18) − N (cid:19) , (14)where Γ( x ) is the gamma function. For bilayer graphene, N = 2, then g ( ∞ ) = 1 .
854 and ∆( T = 0) = 1 . α γ ,whereas, for N (cid:29)
1, then g N ( ∞ ) = 1 and ∆( T = 0) = α g γ . For finite T , the temperature dependence of the or-der parameter (11) is similar to the self-consistent equa-tion for the magnetization in the Weiss mean-field ap-proximation [81] and, for N (cid:29)
1, we find that the criticaltemperature is given by k B T c = α g γ / B. Full band model
We now discuss the results of the numerical calculationto solve the mean-field theory (2,3) for 2 N bands. Theband gap, E g , and the order parameter (the band sepa-ration at k = 0), ∆, are plotted in Fig. 3 as a function (a) � g (b) � g E g ( m e V ) � ( m e V ) k g / k c FIG. 3. The band gap of rhombohedral multilayer graphene atzero temperature as a function of layer number N for differentvalues of interaction strength α g showing (a) the band gap E g and (b) the order parameter ∆. Solid lines are a guide for theeye. Red squares in (a) show the results of Ref. 64 for N = 3to N = 8 obtained using density functional theory. The insetof (a) shows the position of the band gap k g in units of k c asa function of N for α g = 0 . of layer number N for different interaction strengths α g .Red squares in Fig. 3(a) show the results of Ref. 64 for E g for N = 3 to N = 8 obtained using DFT. Our resultsare in qualitative agreement and, by choosing α g = 0 . E g , grows for small N , until it peaks around N = 6 and, then, falls for larger N . The order parame-ter, ∆, also grows for small N , but it saturates for larger N values (as in the simple two-band model).The decrease of E g at large N is largely due to thenon-interacting band structure, as described by Hamil-tonian (1), in that the position of the band gap movesfrom k ≈ N = 2 to k ≈ k c for N (cid:29) k ≈ k c for N (cid:29) k ) is plotted in Fig. 4(a) as a function of the magnitudeof the wave vector k plotted from the valley center (theBerry curvature is isotropic in the minimal model) for N = 16. The position of the maximum in Ω( k ) is givenby k max ≈ s N k c (∆ /γ ) /N , s N = (( N − / (2 N +4)) / N ,which moves from k = 0 to k = k c as N increases [17,18, 82, 83]. The integral of Ω( k ) with respect to thewave vector area also increases with N , as characterizedby the Chern number per flavor which has magnitude N/ T ) is shown in Fig. 5(a) for different N values and interaction strength α g = 0 .
3. We fit ∆( T )to a form suggested in Ref. 64,∆( T )∆(0) = (cid:20) A (cid:18) − TT c (cid:19) + (3 − A ) (cid:18) − TT c (cid:19) +( A − (cid:18) − TT c (cid:19) (cid:21) / , (15) (b) � / k c k/k c Bernalfault (8,8) oddparity(a) � / k c k/k c RMGoddparity (d) � / k c k/k c Bernalfault (12,4) evenparity(c) � / k c k/k c Bernalfault (12,4) oddparity
FIG. 4. Berry curvature Ω( k ) as a function of the magnitudeof the wave vector k plotted from the valley center (the Berrycurvature is isotropic in the minimal model) for N = 16, α g = 0 . T = 0 K. (a) is for pristine RMG in the oddparity ground state, (b) is for RMG with a Bernal stackingfault at its center in the odd parity ground state, (c) is foran off-centre Bernal stacking fault in the odd parity groundstate, (d) is for an off-centre Bernal stacking fault in an evenparity state (which is not the ground state). with A and T c as fitting parameters. As shown inFig. 5(a), the quality of this fit is generally excellent.The resulting layer dependence of the critical tempera-ture T c ( N ) is shown in Fig. 5(b). Red squares in Fig. 5(b)show the results of Ref. 64 for N = 3 to N = 8 obtainedusing density functional theory (DFT); our results are inclose agreement. As in the simple two-band model, T c saturates at a finite value for N (cid:29)
1, and, for α g = 0 . T = 0) / ( k B T c ) ∼ T = 0) / ( k B T c ) = 4 forthe two-band model. IV. INTERLAYER DISORDER
In this section, we consider the influence of interlayerdisorder on the interacting mean-field state in RMG, i.e. we preserve translational invariance within each graphenelayer, but take into account random tight-binding param-eters in the perpendicular-to-layer direction. This hasa close analogy with studies of the SSH model [19, 20]wherein the effects on the zero-energy edge states ofchiral-symmetry-preserving or breaking disorder are con-sidered [65–69]. A major difference here is that we con-sider the influence of disorder on the interacting mean-field state in which the exchange potential has already (a) (b) T c ( K ) � ( m e V ) FIG. 5. (a) the order parameter ∆( T ) of rhombohedral multi-layer graphene as a function of temperature T for interactionstrength α g = 0 . N = 2 , , , , , , T c as a function of layer number N for interaction strength α g = 0 . N = 3 to N = 8 obtained using densityfunctional theory. broken chiral symmetry and gapped the spectrum.We take into account two types of disorder. The first isdue to random layer energies which break chiral symme-try. For a given realization of disorder, the diagonal ele-ments of the non-interacting Hamiltonian (1) take values H A n A n = H B n B n = δ n for layer index n = 1 , , , . . . , N where each δ n takes a random value uniformly distributedin the range [ − δ, δ ] for disorder strength δ . We considerweak disorder up to δ = 10 meV so that δ (cid:28) { E g , γ } for typical values of the band gap E g . Figure 6(a) showsthe mean band gap E g for N = 12 layers as a function ofdisorder strength δ for the odd parity antiferromagneticstate, Fig. 6(b) shows the mean order parameter. Eachdata point (triangles) is an average over twenty differentrealizations of disorder; error bars increase with disorderin Fig. 6 because the standard deviation increases whilethe number of realizations is constant. The mean valuesof both E g and ∆ decrease with disorder δ , although theyappear to be quite robust for weak disorder. We restrictthe study to weak disorder values because the groundstate will change ( e.g. to an odd ferrimagnetic state asmodeled in bilayer graphene [46]) for certain realizationsat higher disorder.The second type of disorder is due to random val-ues of the interlayer coupling γ which preserve chi-ral symmetry. For a given realization of disorder, ele-ments of the non-interacting Hamiltonian (1) describinginterlayer coupling take random values, i.e. H B A = H A B = γ + δ , H B A = H A B = γ + δ , etc. , n = 1 , , , . . . , N −
1, where each δ n takes a randomvalue uniformly distributed in the range [ − δ, δ ] for dis-order strength δ . Figure 6 shows the dependence of themean values of E g and ∆ for N = 12 layers as a functionof disorder strength δ for the odd parity antiferromag-netic state, each data point (circles) is an average overtwenty different realizations of disorder. We find that themean value of ∆ is not affected by disorder (within the (a) (b) E g ( m e V ) � ( m e V ) �� (meV) �� (meV) FIG. 6. The mean band gap of rhombohedral multilayergraphene with N = 12 layers in the odd parity antiferromag-netic state at zero temperature and α g = 0 . δ showing (a) the mean band gap E g and (b)the mean order parameter ∆. Triangles show data for ran-dom layer energies, circles show data for random interlayercoupling. Mean values and error bars are found by averag-ing over twenty disorder realizations. Note that the scale onthe vertical axes is offset from zero, and is a factor of threedifferent between (a) and (b). error bars), and that disorder slightly reduces the meanvalue of E g for the weak disorder values we consider. Thisis in line with studies of the SSH model [65–69] whereone expects chiral-preserving disorder to have a negli-gible effect on the zero-energy edge states, although theexchange interaction has already broken chiral symmetryin the interacting mean-field state considered here.The influence of disorder may be understood by consid-ering the form of the two-band Hamiltonian (7). For ran-dom layer energies, the energies of the outer layers wouldappear directly in the two-band Hamiltonian as randomdiagonal elements for H A A and H B N B N , having a directimpact on the exchange potential and the band gap in theform of random numbers δ and δ N . Interlayer coupling,however, appears in the off-diagonal term as H A B N = − γ ( − κ † ) N = − ( − ¯ hv [ ξk x − ik y ]) N γ − N , i.e. the con-nection between the surface states involves a product ofthe N − γ + δ n for n = 1 , , , . . . , N − N (cid:29)
1, the system self-averages so that the effect ofrandom values δ n is negligible for the low-energy bands. V. BERNAL STACKING FAULT INRHOMBOHEDRAL GRAPHENE
Stacking faults have been considered previously ingraphene multilayer systems [15, 84] and in RMG in par-ticular [32, 70, 71]. Single stacking faults, e.g. a Bernalfault or a twin boundary fault, Fig. 1, within a largeRMG system are interesting because they introduce ad-ditional flat bands into the energy spectrum. Here wefocus on the properties of the interacting ground state athalf filling. (a) (b) E ( e V ) E ( e V ) k/k c k/k c (c) A1 B1 A2 B2 A3 B3 A4 B4 A5 B5 A6 B6 A7 B7 A8 B8
Site W ( m e V ) FIG. 7. Low-energy band structure of RMG with N = 16layers and a Bernal stacking fault at its center with (a) non-interacting electrons described by Hamiltonian (1) and (b)interacting electrons described by the mean-field theory (2,3)with interaction strength α g = 0 .
3. Black lines are bulkbands, blue lines are almost doubly-degenerate low-energybands for k < ∼ k c where ¯ hvk c = γ (8-fold degenerate tak-ing into account spin and valley degrees of freedom). (c) Thevalues of the exchange potential W k =0 σXX at the valley center k = 0 and for each site X (cid:48) = X for RMG with N = 8 layersand a Bernal stacking fault at its center [77]. The black solidline shows the exchange for two flavors, the blue dashed lineis exchange for the other two flavors. A. Band structure of non-interacting electrons
In order to describe the influence of a stacking fault onthe low-energy band structure, we begin by consideringnon-interacting electrons. The number of zero energystates at k = 0 is determined by the stacking structureof the multilayer, particularly the degree of hybridization of p z orbitals (one per site) caused by interlayer coupling γ with a neighboring site in an adjacent layer directlyabove or below. For n atoms coupled in a vertical line byinterlayer coupling γ , even n ( e.g. a dimer) contributes n bulk bands, but no zero energy states at k = 0. For odd n ( e.g. a monomer not directly connected to a neighborin an adjacent layer or a trimer), there are n − k = 0 [85].For RMG, all sites are part of a dimer apart from A and B N at the surfaces which are monomers: hence thereare two zero energy states as shown in Fig. 2. For aBernal-stacked multilayer with N layers, there is one N -mer contributing one zero energy state if N is odd andin addition there are N monomers. Thus, overall, thereare N zero energy states if N is even, N + 1 if N is odd.We consider RMG with a single Bernal stacking fault,as illustrated in Fig. 1 (central panel) for N = 7 layers.Specifically, we use integers ( m, n ) to denote a rhombohe-dral section with m layers and sites A , B , . . . , A m , B m connected by a Bernal stacking fault to a rhombohedralsection with n layers and sites A m +1 , B m +1 , . . . A N , B N where the total layer number is N = m + n and m ≥ n ≥
2. Thus, the example in the central panel of Fig. 1 isa (3 ,
4) Bernal stacking fault. At the stacking fault, thereare four vertically connected atomic sites (sites B , A , B , A in Fig. 1) which make up a 4-mer; this is even,so it contributes 4 bulk bands, but no zero energy statesat k = 0. Rather, the zero energy states arise from thesites not directly connected to a neighbor in an adjacentlayer, namely A , B m , A m +1 , B N (sites A , B , A , B in Fig. 1), so there are four low-energy states per spinand valley flavor.The low-energy bands of non-interacting electrons for N = 16 layers with a Bernal stacking fault at the center n = m = 8 are shown in Fig. 7(a). For m, n (cid:29) m and n layers, respectively. This may be understood byderiving an effective low-energy four band Hamiltonian,following the procedure described previously for bilayergraphene [78, 86], in a basis of orbitals on A , B m , A m +1 , B N sites. For k (cid:28) k c and E (cid:28) γ , we find H ( m,n )Bernal = γ − ( − κ † ) m − κ † ) m + n − − ( − κ ) m − c mn ( k/k c ) (cid:96) − κ − c mn ( k/k c ) (cid:96) − ( κ † ) − ( − κ † ) n ( − κ ) m + n − − ( − κ ) n , (16)where c mn = (1 + δ mn ) / (cid:96) = min( m, n ). As theHamiltonian is chiral, every matrix element between twoA sites or between two B sites is zero. For the non-zero elements (between A and B sites), we keep only theleading terms in k/k c .The diagonal 2 × m and n layers; the off-diagonal2 × − κ ) m + n − describes effective coupling betweenthe A and B N sites which are on opposite surfaces of thesample and this is very small for N (cid:29) k/k c (cid:28) − c mn ( k/k c ) (cid:96) − ( κ † ) describes effective coupling (a) � g (b) � g E g ( m e V ) � � ( m e V ) FIG. 8. The band gap of RMG with a Bernal stacking faultat its center, for even N , at zero temperature as a functionof layer number N for different values of interaction strength α g showing (a) the band gap E g and (b) the order parameter∆ . Points are data for the system with the stacking fault,solid lines are data from Fig. 3 for a single RMG section with m = N/ between the B m and A m +1 sites. Although they are onadjacent layers, this coupling is of order ( k/k c ) (cid:96) whichis also very small for (cid:96) (cid:29) k/k c (cid:28)
1. Weak effectivecoupling between the B m and A m +1 sites arises from thefact that the Bernal stacking fault consists of four ver-tically coupled sites ( B , A , B , A in Fig. 1). Theireffective coupling in the basis of Eq. (16) is described byinverting the 4 × − = −
11 0 0 00 0 0 1 − , which has an exactly zero matrix element between thesecond and third components. The reason that matrixelement − c mn ( k/k c ) (cid:96) − ( κ † ) in Eq. (16) is not alsoidentically zero is that this small contribution arises froma slight rotation of the low-energy basis states that isrequired to preserve their normalization [86].Since the four band Hamiltonian (16) is chiral, theenergy spectrum of non-interacting electrons displayselectron-hole symmetry and the band energies E maybe determined as the solution of a quadratic equation:( E/γ ) = 12 β B ( k/k c ) ± (cid:113) β ( k/k c ) − η B ( k/k c ) ,β B ( x ) = x m + x n + x m +2 n − + c mn x (cid:96) ,η B ( x ) = x m +2 n + 2 c mn x m +2 n +2 (cid:96) − + c mn x m +2 n +4 (cid:96) − . Of particular interest is when the fault lies exactly in thecenter of a long RMG system: n = m = N/ c mn = 1,and (cid:96) = m = N/ N (cid:29)
1. Then, the stacking fault(as described by the off-diagonal 2 × E ≈ ± γ (cid:2) ( k/k c ) N/ ± ( k/k c ) N − (cid:3) for k (cid:28) k c . (a) (b) E ( e V ) E ( e V ) k/k c k/k c (c) A1 B1 A2 B2 A3 B3 A4 B4 A5 B5 A6 B6 A7 B7 A8 B8
Site W ( m e V ) FIG. 9. Low-energy band structure of RMG with N = 16layers and a Bernal stacking fault between its third and fourthlayers, ( m, n ) = (3 , α g = 0 .
3. Black lines are bulk bands, blue lines aresurface bands for k < ∼ k c where ¯ hvk c = γ . (c) The values ofthe exchange potential W k =0 σXX at the valley center k = 0and for each site X (cid:48) = X for RMG with N = 8 layers and aBernal stacking fault off center [77] between the second andthird layers, ( m, n ) = (2 , B. Numerical mean-field theory
The numerical mean-field theory calculations proceedas for pristine RMG, except that the stacking fault istaken into account by a different position of the inter-layer coupling γ in the non-interacting Hamiltonian (1).We find a number of self-consistent solutions with differ-ent exchange profiles including both even and odd paritywithin each flavor, and ferromagnetic, ferrimagnetic orantiferromagnetic arrangements of flavors. However, ondetermining the total energy of each, we find that theground state is the odd antiferromagnetic state, a gener-alization of the ground state in pristine RMG (discussedin Section III). That the odd antiferromagnetic state isthe ground state agrees with Refs. 60 and 61 for the (2 , N = 4 Bernal-stackedmultilayer. For benchmarking, with α g = 0 . E g = 1 .
42 meV for the (2 ,
2) fault which compares with E g ≈ .
44 meV for the second red circle in Fig. 3(a) ofRef. 60.Since there are now four low-energy bands, we identifythe separation at k = 0 of the lowest conduction band andthe highest valence band as ∆ and the separation at k =0 of the second lowest conduction band and the second (a) (b) � ( m e V ) E g � � m � � � ( m e V ) E g � � m � � FIG. 10. (a) The band gap E g and order parameters ∆ , ∆ of RMG with N = 16 layers and a Bernal stacking fault ofstructure ( m, − m ) consisting of a rhombohedral section oflength m coupled to a rhombohedral section of length 16 − m .The interaction strength is α g = 0 .
3. Points are data forthe system with the stacking fault, solid lines are data fromFig. 3 for a single RMG section with m layers for ∆ and E g , and 16 − m layers for ∆ . (b) The band gap E g andorder parameters ∆ , ∆ of RMG with N = 15 layers and atwin boundary fault of structure ( m, − m ) consisting of arhombohedral section of length m coupled to a rhombohedralsection of length 16 − m . The interaction strength is α g = 0 . m layers for ∆ and E g , and 16 − m layers for ∆ . highest valence band as ∆ , ∆ ≥ ∆ . For m, n (cid:29)
1, wefind that the low-energy bands of the interacting systembehave as if they arise from two disconnected pieces ofRMG. This is not surprising as, due to their chirality, thewave functions of the B m and A m +1 sites have differentdependences on the polar angle in the graphene planeand, at k = 0, this suppresses the exchange interactionmatrix element between them.In the special case of an even number of layers N = 2 m with a stacking fault on the central layer ( m, m ), the low-energy bands are almost doubly degenerate (i.e. 8-folddegenerate taking into account spin and valley), as shownfor N = 16 in Fig. 7(a,b) for non-interacting and inter-acting electrons, respectively. The exchange potential W k =0 σXX at the valley center k = 0 and for each site X (cid:48) = X for RMG with N = 8 layers and a Bernal faultat its center is shown in Fig. 7(c) [77]. This is a gen-eralization of the odd antiferromagnetic state in pristineRMG, but now the exchange (and carrier density) perflavor has a substantial magnitude on sites B m , A m +1 by the stacking fault as well as the surface states. Theband gap, E g , and the order parameter ∆ are plottedas data points in Fig. 8 as a function of layer number N for different interaction strengths α g (for clarity, wedon’t plot ∆ because ∆ ≈ ∆ when the fault is at thecenter). To illustrate that the system behaves almost astwo separate RMG sections of m = N/ N/ N >
4, theagreement is very close. More generally, the stacking fault breaks spatial inver-sion symmetry. As an example, bands for a N = 16 layersystem ( m, n ) = (3 ,
13) with a three-layer section con-nected to a 13-layer section are shown in Fig. 9(a,b) fornon-interacting and interacting electrons, respectively.Although spatial inversion symmetry is absent, the bandstructure has flavor degeneracy for the antiferromagneticground state. This is not generally the case, e.g. withina ferrimagnetic configuration, the flavors with differentorientation are usually not degenerate. For the interact-ing case, Fig. 9(b) shows the bands of a single flavor forthe antiferromagnetic ground state; within a flavor, thefour low-energy bands are not degenerate. In Fig. 9(b),∆ > ∆ , where ∆ (the separation at k = 0 of the low-est conduction band and the highest valence band) is theorder parameter related to the short section m = 3, ∆ (the separation at k = 0 of the second lowest conductionband and the second highest valence band) is related tothe long section n = 13.The exchange potential W k =0 σXX at the valley center k = 0 and for each site X (cid:48) = X for RMG with N = 8 lay-ers and a Bernal fault off center, ( m, n ) = (2 , A , B , A , B ). The exchange potential(and carrier density) per flavor has a larger magnitude onsites with low-energy orbitals A , B , A , B sites, buthas a much larger magnitude on sites A , B associatedwith the longer RMG section than A , B related to thesmall section. This is reflected in the relative magnitudesof ∆ and ∆ .Figure 10(a) shows E g , ∆ and ∆ for a N = 16 sys-tem ( m, − m ) plotted as data points as a function ofthe number of layers m in the short section ( i.e. for dif-ferent positions of the stacking fault). The solid linesare not fits, but they show data taken from Fig. 3 fora single RMG section: a section of length m is used tocompare E g and ∆ , a section of length 16 − m is usedto compare ∆ . For ∆ and ∆ the agreement is close, itis slightly less close for the band gap E g ; ∆ and ∆ aredetermined at k = 0 where the approximate splitting ofthe system into two parts is clearcut whereas E g is gener-ally determined at non-zero k . Weak coupling due to theBernal stacking fault is indicated by the close agreementof the data points and lines for ∆ and ∆ , as well as thenearly equal values of ∆ and ∆ for the spatially sym-metric case m = 8 (∆ and ∆ differ by about 1 meV);the stacking fault only breaks the degeneracy slightly.For the odd parity ground state, the Berry curvatureΩ( k ) is plotted in Fig. 4(b) for N = 16 with a stack-ing fault at the center ( m, n ) = (8 , k ) for the oddparity ground state with the stacking fault off-center0 (a) (b) E ( e V ) E ( e V ) k/k c k/k c (c) A1B1 A2 B2 A3B3 A4B4 A5 B5 A6B6 A7 B7A8 B8 A9 B9
Site W ( m e V ) FIG. 11. Low-energy band structure of RMG with N = 15layers and a twin boundary stacking fault at its center with(a) non-interacting electrons described by Hamiltonian (1)and (b) interacting electrons described by the mean-field the-ory (2,3) with interaction strength α g = 0 .
3. Black linesare bulk bands, blue lines are doubly-degenerate low-energybands for k < ∼ k c where ¯ hvk c = γ (8-fold degenerate takinginto account spin and valley degrees of freedom). (c) The val-ues of the exchange potential W k =0 σXX at the valley center k = 0 and for each site X (cid:48) = X for RMG with N = 9 layersand a twin boundary stacking fault at its center [77]. Theblack solid line shows the exchange for two flavors, the bluedashed line is exchange for the other two flavors. ( m, n ) = (12 , m = 12and n = 4. For both of these cases, the Berry cur-vature Ω( k ) sums to give a Chern number of magni-tude N/ k ) for the same system in an even par-ity state [the exchange has relative signs of (-,+,+,-) onthe four low-energy orbitals ( A , B m , A m +1 , B N )] which is not the ground state. In this case, the contributionsof the two sections of length m = 12 and n = 4 appearwith an opposite sign and the Berry curvature Ω( k ) sumsto give a Chern number of magnitude | m − n | / i.e. itdepends on the position of the stacking fault. VI. TWIN BOUNDARY STACKING FAULT INRHOMBOHEDRAL GRAPHENEA. Band structure of non-interacting electrons
As a second example, we consider RMG with a twinboundary stacking fault, as illustrated in Fig. 1 (rightpanel) for N = 7 layers. Specifically, we use integers( m, n ) to denote a rhombohedral section with m lay-ers and sites A , B , . . . , A m , B m connected by a twinboundary stacking fault to a rhombohedral section with n layers and sites B m , A m , . . . B N , A N where m ≥ n ≥
2. The total layer number is N = m + n − A m , B m . Thus, the example in Fig. 1 is a (3 ,
5) twinboundary stacking fault with N = 7 total layers. Thisfault contrasts with the Bernal fault. For example, at thestacking fault there are only three vertically connectedatomic sites (sites B , A , B , in Fig. 1) which make upa trimer; this is odd so it contributes two bulk bandsand one band near zero energy at k = 0 [32] in a simi-lar way to trilayer Bernal graphene [87]. Thus, overall,there are four low energy bands per spin and valley re-lated to three sites not directly connected to a neighborin an adjacent layer, namely A , B m , A N (sites A , B , A in Fig. 1) plus an odd combination of trimer sites( B m − − B m +1 ) / √ B − B ) / √ B m and to ( B m − − B m +1 ) / √ A , ( B m − − B m +1 ) / √ A N , B m sites, is given by H ( m,n )twin = γ − ( − κ † ) m − / √ − ( − κ † ) m / − ( − κ ) m − / √ − κ ) n − / √ − κ † ) n − / √ − ( − κ † ) n / − ( − κ ) m / − ( − κ ) n / . (17)As the Hamiltonian is chiral, every matrix element be-tween two A sites or between two B sites is zero. Forthe non-zero elements (between A and B sites), we keeponly the leading terms in k/k c . The second and fourthcolumns indicate that the B orbitals at the fault are cou-pled to both of the RMG sections. Since the four band Hamiltonian (17) is chiral, theenergy spectrum of non-interacting electrons displayselectron-hole symmetry and the band energies E are1 (a) (b) � g E g ( m e V ) � � ( m e V ) � g FIG. 12. The band gap of RMG with a twin boundary stack-ing fault at its center, for odd N , at zero temperature as afunction of layer number N for different values of interactionstrength α g showing (a) the band gap E g and (b) the orderparameter ∆ . Points are data for the system with the stack-ing fault, solid lines are data from Fig. 3 for a single RMGsection with ( N − / given by( E/γ ) = 12 β t ( k/k c ) ± (cid:113) β ( k/k c ) − η t ( k/k c ) ,β t ( x ) = 12 (cid:0) x m − + x n − (cid:1) + 14 (cid:0) x m + x n (cid:1) ,η t ( x ) = 12 x m +2 n − . When the fault lies exactly in the center of a long RMGsystem: n = m = ( N + 1) / ≈ N/ N (cid:29) E ≈ ± γ ( k/k c ) N/ and E ≈ ± ( γ / √ k/k c ) N/ .The low-energy dispersion of a pair of the bands acquiresan additional factor of 1 / √ N/ B. Numerical mean-field theory
The numerical mean-field theory calculations proceedas for pristine RMG, except that the stacking fault istaken into account by a different position of the inter-layer coupling γ in the non-interacting Hamiltonian (1).In the special case of an odd number of layers N = 2 m − m, m ), the low-energy bands are almost doubly degenerate (i.e. 8-folddegenerate taking into account spin and valley), as shownfor N = 15 in Fig. 11(a,b) for non-interacting and inter-acting electrons, respectively. Although it is not clearlyvisible in Fig. 11(b), the low-energy bands actually havea small separation of a few meV (this is indicated inFig. 10(b) where there is a non-zero separation of ∆ and ∆ for m = 8).The exchange potential W k =0 σXX at the valley center k = 0 and for each site X (cid:48) = X for RMG with N = 9layers and a fault at its center is shown in Fig. 11(c) [77].This ground state has an even parity of exchange (and (a) (b) E ( e V ) E ( e V ) k/k c k/k c (c) A1 B1 A2 B2 A3 B3 A4 B4 A5 B5 A6 B6 A7 B7 A8 B8
Site W ( m e V ) FIG. 13. Low-energy band structure of RMG with N = 16layers and a twin boundary fault off center at its third layer,( m, n ) = (3 , α g = 0 .
3. Black lines are bulk bands, blue lines are surfacebands for k < ∼ k c where ¯ hvk c = γ . (c) The values of theexchange potential W k =0 σXX at the valley center k = 0 andfor each site X (cid:48) = X for RMG with N = 8 layers and atwin boundary stacking fault off center [77] at its third layer,( m, n ) = (3 , carrier density) per flavor, with the large magnitude ofexchange on the surface orbitals A , A having the samesign, low-energy orbitals at the fault B , B , B have ex-change potentials with the opposite sign. Within the fourspin-valley flavors, the ground state has an antiferromag-netic configuration, because this minimizes the Hartreeenergy as previously.The band gap, E g , and the order parameter ∆ areplotted as data points in Fig. 12 as a function of layernumber N for different interaction strengths α g and afault at the center (for clarity, we don’t plot ∆ because∆ ≈ ∆ when the fault is at the center). The solid linesare data taken from Fig. 3 for a single RMG section with( N − / N − / N/ N = 16 layer system (3 ,
14) with an offcenter fault, namely a three-layer section connected toa 14-layer section are shown in Fig. 13(a,b) for non-interacting and interacting electrons, respectively. Theexchange potential W k =0 σXX at the valley center k = 0and for each site X (cid:48) = X for RMG with N = 8 layers and2 (b) � / k c k/k c twinfault (12,4) evenparity(a) � / k c k/k c twinfault (8,8) evenparity FIG. 14. Berry curvature Ω( k ) as a function of the magnitudeof the wave vector k plotted from the valley center (the Berrycurvature is isotropic in the minimal model) for N = 15, α g = 0 . T = 0 K. (a) is for a twin boundary fault atthe center in the even parity ground state, (b) is for a twinboundary fault off center in the even parity ground state. a fault off center is shown in Fig. 13(c) [77]. This is alsothe even parity state [the exchange has relative signs of(-,+,+,+,-) on the low-energy orbitals ( A , B m − , B m , B m +1 , A N )]. In the interacting case, Fig. 13(b), thefour low-energy bands are clearly not degenerate. Fig-ure 10(b) shows E g , ∆ and ∆ for a N = 15 system( m, − m ) plotted as data points as a function of thenumber of layers m in the short section. The solid linesshow data taken from Fig. 3 for a single RMG section:a section of length m is used to compare E g and ∆ , asection of length 16 − m is used to compare ∆ . Theagreement of the data points and the solid lines is rea-sonable, although not as close as in the case of the Bernalfault Fig. 10(a), this is due to the larger coupling betweenthe two RMG sections in the twin stacking fault case. Inparticular, when the stacking fault is in the center anddoesn’t break spatial inversion symmetry, m = 8, there’sstill a significant difference between ∆ and ∆ (of about8 meV), whereas this difference is small (about 1 meV) inthe Bernal fault case, Fig. 10(a).For the even parity ground state, the Berry curva-ture Ω( k ) is plotted in Fig. 14(a) for N = 15 with astacking fault at the center ( m, n ) = (8 , k ) for theeven parity ground state with the stacking fault off-center( m, n ) = (12 , m = 12and n = 4. For both of these cases, the Berry curva-ture Ω( k ) sums to give a Chern number of magnitude N/ A and B sites within a layer for the layers above thefault (right panel of Fig. 1) compensating the change of relative sign of potential differences for the even paritystate. VII. DISCUSSION
We have generalized the mean-field Hartree Fock de-scription [45, 46, 60–62] to provide a comprehensive qual-itative description of broken symmetry ground states inRMG, including the effects of defects including randomdisorder and stacking faults. The non-trivial topology ofthe low-energy bands is reflected in large Berry curvatureand Chern numbers per spin-valley flavor. An obviousgeneralization is to a number of stacking faults separat-ing rhombohedral sections with different numbers of lay-ers m , each section contributing a pair of low-energy flatbands and a peak in the Berry curvature at a character-istic m -dependent wave vector. In RMG, each stackingfault contributes a pair of low-energy flat bands becausethey are more complicated than the domain walls usuallyconsidered in the SSH model which consist of isolatedmonomers or trimers [20]. The Bernal fault correspondsto two monomers (and a 4-mer), the twin boundary faultis a monomer plus a trimer.As described in Section III B, sources of systematic un-certainty include the choice of cutoff k (cid:63) , the omission oftight-binding parameters in the minimal model, and thevalue of the interaction parameter α g . Additional tight-binding parameters such as γ and γ will introduce trig-onal warping of the dispersion around each valley (sothe Berry curvature, Fig. 4, will be anisotropic), and γ will break particle-hole symmetry [16, 18]; this is likelyto reduce the value of the band gap. The tight-bindingparameters are usually smaller in magnitude than thetypical values of the band gap that we predict, but thereis a possibility that additional parameters will changethe qualitative nature of the ground state [61]. However,even without these parameters, for our choice of cutoffand for α g = 0 .
3, we find close agreement of band gapvalues in RMG with DFT calculations of Ref. 64 (whichconsidered N = 3 to N = 8 layers).The mean-field Hartree Fock approach neglects strongcorrelation effects, and there have been predictions ofmagnetic ordering [88, 89] and superconductivity [90–92] due to the flat bands in RMG. We speculate thatthe additional flat bands localized at stacking faults, andin close spatial proximity to each other, are more likelyto support strongly-correlated states than the widely-separated surface states in pristine RMG. ACKNOWLEDGMENTS
The authors thank Mikito Koshino for discussions.Computer time was provided by the Lancaster Univer-sity High End Computing facility.3 [1] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi,E. Kaxiras, and P. Jarillo-Herrero, Nature , 43(2018).[2] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken,J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T.Taniguchi, E. Kaxiras, R. C. Ashoori, and P. Jarillo-Herrero, Nature , 80 (2018).[3] B. Sutherland, Phys. Rev. B , 5208 (1986).[4] E. H. Lieb, Phys. Rev. Lett. , 1201 (1989).[5] C. Wu, D. Bergman, L. Balents, and S. Das Sarma, Phys.Rev. Lett. , 070401 (2007).[6] H.-M. Guo and M. Franz Phys. Rev. B , 113102 (2009).[7] K. Sun, Z. Gu, H. Katsura, and S. Das Sarma, Phys.Rev. Lett. , 236803 (2011).[8] R. Drost, T. Ojanen, A. Harju, and P. Liljeroth, Nat.Phys. , 668 (2017).[9] D. Leykam, A. Andreanov, and S. Flach, Adv. Phys. X , 677 (2018).[10] S. N. Kempkes, M. R. Slot, J. J. van den Broeke,P. Capiod, W. A. Benalcazar, D. Vanmaekelbergh, D.Bercioux, I. Swart, and C. Morais Smith, Nat. Mater. , 1292 (2019).[11] M. Kang, S. Fang, L. Ye, H. C. Po, J. Denlinger, C.Jozwiak, A. Bostwick, E. Rotenberg, E. Kaxiras, J. G.Checkelsky, and R. Comin, Nat. Commun. , 4004(2020).[12] J. W. McClure, Carbon , 425 (1969).[13] Latil and L. Henrard, Phys. Rev. Lett. , 036803 (2006).[14] M. Aoki and H. Amawashi, Solid State Commun. ,123 (2007).[15] D. P. Arovas and F. Guinea, Phys. Rev. B , 245416(2008).[16] M. Koshino and E. McCann, Phys. Rev. B , 165409(2009).[17] F. Zhang, J. Jung, G. A. Fiete, Q. Niu, and A. H. Mac-Donald, Phys. Rev. Lett. , 156801 (2011).[18] S. Slizovskiy, E. McCann, M. Koshino, and V. I. Fal’ko,Commun. Phys. , 164 (2019).[19] W. P. Su, J. R. Schrieffer, and A. J. Heeger, Phys. Rev.Lett. , 1698 (1979).[20] J. K. Asb´oth, L. Oroszl´any, and A. P´alyi, A Short Courseon Topological Insulators (Springer, Switzerland, 2016).[21] T. T. Heikkil¨a and G. E. Volovik, JETP Lett. , 59(2011).[22] R. Xiao, F. Tasn´adi, K. Koepernik, J. W. F. Vender-bos, M. Richter, and M. Taut, Phys. Rev. B , 165404(2011).[23] D. Pierucci, H. Sediri, M. Hajlaoui, J.-C. Girard, T.Brumme, M. Calandra, E. Velez-Fort, G. Patriarche, M.G. Silly, G. Ferro, V Souli`ere, M. Marangolo, F. Sirotti,F. Mauri, and A. Ouerghi, ACS Nano , 5432 (2015).[24] Y. Henni, H. P. O. Collado, K. Nogajewski, M. R. Molas,G. Usaj, C. A. Balseiro, M. Orlita, M. Potemski, and C.Faugeras, Nano Lett. , 3710 (2016).[25] H. Henck, J. Avila, Z. B. Aziza, D. Pierucci, J. Baima,B. Pamuk, J. Chaste, D. Utt, M. Bartos, K. Nogajewski,B. A. Piot, M. Orlita, M. Potemski, M. Calandra, M. C.Asensio, F. Mauri, C. Faugeras, and A. Ouerghi, Phys.Rev. B , 245421 (2018).[26] T. Latychevskaia, S.-K. Son, Y. Yang, D. Chancellor, M.Brown, S. Ozdemir, I. Madan, G. Berruto, F. Carbone, A. Mishchenko and K. S. Novoselov, Front. Phys. ,13608 (2019).[27] Y. Yang, Y.-C. Zou, C. R. Woods, Y. Shi, J. Yin, S. Xu,S. Ozdemir, T. Taniguchi, K. Watanabe, A. K. Geim,K. S. Novoselov, S. J. Haigh, and A. Mishchenko, NanoLett. , 8526 (2019).[28] F. R. Geisenhof, F. Winterer, S. Wakolbinger, T. D.Gokus, Y. C. Durmaz, D. Priesack, J. Lenz, F. Keil-mann, K. Watanabe, T. Taniguchi, R. Guerrero-Avil´es,M. Pelc, A. Ayuela and R. Thomas Weitz, ACS Appl.Nano Mater. , 6067 (2019).[29] Y. Lee, S. Che, J. Velasco Jr., D. Tran, J. Baima,F. Mauri, M. Calandra, M. Bockrath, and C. N. Lau,arXiv:1911.04450[30] C. Bouhafs, S. Pezzini, N. Mishra, V. Mi˘seikis, Y. Niu,C. Struzzi, A. A. Zakharov, S. Forti, and C. Coletti,arXiv:2006.06667[31] A. Kerelsky, C. Rubio-Verd´u, L. Xian, D. M. Kennes, D.Halbertal, N. Finney, L. Song, S. Turkel, L. Wanga, K.Watanabe, T. Taniguchi, J. Hone, C. Dean, D. N. Basov,A. Rubio, and A. N. Pasupathy, Proc. Natl. Acad. Sci.U.S.A. , 2017366118 (2021).[32] Y. Shi, S. Xu, Y. Yang, S. Slizovskiy, S. V. Morozov, S.-K. Son, S. Ozdemir, C. Mullan, J. Barrier, J. Yin, A. I.Berdyugin, B. A. Piot, T. Taniguchi, K. Watanabe, V. I.Fal’ko, K. S. Novoselov, A. K. Geim, and A. Mishchenko,Nature , 210 (2020).[33] R. T. Weitz, M. T. Allen, B. E. Feldman, J. Martin, andA. Yacoby, Science , 812 (2010).[34] F. Freitag, J. Trbovic, M. Weiss, and C. Sch¨onenberger,Phys. Rev. Lett. , 076602 (2012)[35] J. Velasco Jr, L. Jing, W. Bao, Y. Lee, P. Kratz, V.Aji, M. Bockrath, C. N. Lau, C. Varma, R. Stillwell, D.Smirnov, F. Zhang, J. Jung, and A. H. MacDonald, Nat.Nanotechnol. , 156 (2012).[36] W. Bao, J. Velasco Jr., F. Zhang, L. Jing, B. Standley,D. Smirnov, M. Bockrath, A. H. MacDonald, and C. N.Lau, Proc. Natl. Acad. Sci. , 10802 (2012).[37] A. Veligura, H. J. van Elferen, N. Tombros, J. C. Maan,U. Zeitler, and B. J. van Wees, Phys. Rev. B , 155412(2012).[38] F. Freitag, M. Weiss, R. Maurand, J. Trbovic, and C.Sch¨onenberger, Phys. Rev. B , 161402(R) (2013).[39] A. L. Grushina, D.-K. Ki, M. Koshino, A. A. L. Nico-let, C. Faugeras, E. McCann, M. Potemski, and A. F.Morpurgo, Nat. Commun. , 6419 (2015).[40] Y. Nam, D.-K. Ki, M. Koshino, E. McCann, and A. FMorpurgo, 2D Mater. , 324 (2018).[42] W. Bao, L. Jing, J. Velasco Jr, Y. Lee, G. Liu, D. Tran,B. Standley, M. Aykol, S. B. Cronin, D. Smirnov, M.Koshino, E. McCann, M. Bockrath and C. N. Lau, Nat.Phys. , 948 (2011).[43] Y. Lee, D. Tran, K. Myhro, J. Velasco, N. Gillgren, C.N. Lau, Y. Barlas, J. M. Poumirol, D. Smirnov and F.Guinea, Nat. Commun. , 5656 (2014).[44] K. Myhro, S. Che, Y. Shi, Y. Lee, K. Thilahar, K. Bleich,D. Smirnov and C. N. Lau, 2D Mater. , 045013 (2018).[45] H. Min, G. Borghi, M. Polini, and A. H. MacDonald,Phys. Rev. B , 041407(R) (2008). [46] J. Jung, F. Zhang, and A. H. MacDonald, Phys. Rev. B , 115408 (2011).[47] J. Nilsson, A. H. Castro Neto, N. M. R. Peres, and F.Guinea, Phys. Rev. B , 214418 (2006).[48] K. Sun, H. Yao, E. Fradkin, and S. A. Kivelson Phys.Rev. Lett. , 046811(2009).[49] O. Vafek and K. Yang, Phys. Rev. B , 041401(R)(2010).[50] R. Nandkishore and L. Levitov, Phys. Rev. Lett. ,156803 (2010).[51] F. Zhang, H. Min, M. Polini, and A. H. MacDonald,Phys. Rev. B , 041402(R) (2010).[52] R. Nandkishore and L. Levitov Phys. Rev. B , 115124(2010).[53] O. Vafek, Phys. Rev. B , 205106 (2010).[54] Y. Lemonik, I. L. Aleiner, C. Toke, and V. I. Fal’ko, Phys.Rev. B , 201408(R) (2010).[55] A. H. MacDonald, J. Jung and F. Zhang, Phys. Scr. T146 , 014012 (2012).[56] Y. Lemonik, I. Aleiner, and V. I. Fal’ko, Phys. Rev. B , 245451 (2012).[57] M. Kharitonov, Phys. Rev. B , 195435 (2012).[58] V. Cvetkovic, R. E. Throckmorton, and O. Vafek, Phys.Rev. B , 075467 (2012).[59] L. Zhu, V. Aji, and C. M. Varma, Phys. Rev. B ,035427 (2013).[60] C. Yoon, Y. Jang, J. Jung and H. Min, 2D Mater. ,021025 (2017).[61] M. Koshino, K. Sugisawa, and E. McCann, Phys. Rev. B , 235311 (2017).[62] J. Jung and A. H. MacDonald, Phys. Rev. B , 075408(2013).[63] J. Jia, E. V. Gorbar, and V. P. Gusynin, Phys. Rev. B , 205428 (2013).[64] B. Pamuk, J. Baima, F. Mauri, and M. Calandra, Phys.Rev. B , 075422 (2017).[65] I. Mondragon-Shem, T. L. Hughes, J. Song, and E. Pro-dan, Phys. Rev. Lett. , 046802 (2014).[66] T. Liu and H. Guo, Phys. Lett. A , 3287 (2018).[67] B. P´erez-Gonz´alez, M. Bello, ´A. G´omez-Le´on, and G.Platero, Phys. Rev. B , 035146 (2019).[68] C. J¨urß and D. Bauer, Phys. Rev. B , 195428 (2019).[69] M. Scollon and M. P. Kennett Phys. Rev. B , 144204(2020).[70] M. Taut, K. Koepernik, and M. Richter, Phys. Rev. B , 085312 (2014).[71] A. Garc´ıa-Ruiz, S. Slizovskiy, M. Mucha-Kruczy´nski, andV. I. Fal’ko, Nano Lett. , 6152 (2019). [72] T. Fukui, Y. Hatsugai, and H. Suzuki, J. Phys. Soc.Japan , 1674 (2005).[73] We calculate the Chern number per spin-valley flavor byintegrating the non-Abelian Berry curvature [72] for theoccupied valence bands in the vicinity of a single val-ley. Since this is not the actual Chern number (obtainedby integrating over the whole Brillouin zone), the Chernnumber per flavor may take half integer values.[74] A. B. Kuzmenko, I. Crassee, D. van der Marel, P. Blakeand K. S. Novoselov, Phys. Rev. B , 165406 (2009).[75] We use Gaussian units. For SI units, V ( q ; z ) should havean additional factor of 1 / (4 π(cid:15) ).[76] Y. Jang, E. H. Hwang, A. H. MacDonald and H. Min,Phys. Rev. B , 041411 (2015).[77] When plotting the exchange interaction W k =0 σXX , wesubtract a constant term so that the mean value is zero.[78] E. McCann and V. I Fal’ko, Phys. Rev. Lett. , 086805(2006).[79] K. S. Novoselov, E. McCann, S. V. Morozov, V. I. Fal’ko,M. I. Katsnelson, U. Zeitler, D. Jiang, F. Schedin, andA. K. Geim, Nat. Phys. , 177 (2006).[80] M. Koshino, Phys. Rev. B , 125304 (2010).[81] C. Kittel, Introduction to Solid State Physics (Wiley,Hoboken, N.J., 2005, 8th Ed.)[82] D. Xiao, W. Yao, and Q. Niu, Phys. Rev. Lett. , 236809(2007).[83] D. Xiao, M.-C. Chang, and Q. Niu, Rev. Mod. Phys. ,1959 (2010).[84] M. Koshino and E. McCann, Phys. Rev. B , 045420(2013).[85] H. Min and A. H. MacDonald, Phys. Rev. B , 155416(2008).[86] E. McCann and M. Koshino, Rep. Prog. Phys. , 056503(2013).[87] M. Koshino and E. McCann, Phys. Rev. B , 125443(2009).[88] M. Otani, M. Koshino, Y. Takagi, and S. Okada, Phys.Rev. B , 161403(R) (2010).[89] D.-H. Xu, J. Yuan, Z.-J. Yao, Y. Zhou, J.-H. Gao, andF.-C. Zhang, Phys. Rev. B , 201404(R) (2012).[90] N. B. Kopnin, M. Ij¨as, A. Harju, and T. T. Heikkil¨a,Phys. Rev. B , 140503(R) (2013).[91] W. A. Mu˜noz, L. Covaci, and F. M. Peeters, Phys. Rev.B , 134509 (2013).[92] T. L¨othman and A. M. Black-Schaffer, Phys. Rev. B96