Interacting Dirac fermions under spatially alternating pseudo-magnetic field: Realization of spontaneous quantum Hall effect
IInteracting Dirac fermions under spatially alternating pseudo-magnetic field:Realization of spontaneous quantum Hall effect
J¨orn W. F. Venderbos and Liang Fu
Department of Physics, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA (Dated: May 7, 2018)Both topological crystalline insulators surfaces and graphene host multi-valley massless Diracfermions which are not pinned to a high-symmetry point of the Brillouin zone. Strain couplesto the low-energy electrons as a time-reversal invariant gauge field, leading to the formation ofpseudo-Landau levels (PLL). Here we study periodic pseudo-magnetic fields originating from strainsuperlattices. We study the low-energy Dirac PLL spectrum induced by the strain superlattice andanalyze the effect of various polarized states. Through self-consistent Hartree-Fock calculations weestablish that, due to the strain superlattice and PLL electronic structure, a valley-ordered statespontaneously breaking time-reversal and realizing a quantum Hall phase is favored, while othersare suppressed. Our analysis applies to both topological crystalline insulators and graphene.
I. INTRODUCTION
The discovery of graphene and topological insula-tors has significantly boosted the ubiquity of condensedmatter realizations of Dirac fermions as emergent elec-tronic excitations at low-energy [1–3]. Dirac electronsin condensed matter systems have enjoyed an enor-mous amount of interest both from a fundamental andtechnological application perspective [4]. A key differ-ence between graphene and topological insulators is thenumber of species, or valleys, of Dirac fermions andtheir locations in momentum space. Topological insu-lators (TI) protected by time reversal symmetry host asingle-valley Dirac fermion, which is pinned to a time-reversal-invariant (TRI) momentum in the surface Bril-louin zone [5]. In contrast, graphene hosts two valleys ofDirac fermions located at non-TRI momenta [6, 7], eachvalley having an additional spin degeneracy. More re-cently, a new type of Dirac fermions was discovered on thesurface of topological crystalline insulators (TCI) SnTe,(Sn,Pb)Se and (Sn,Pb)Te [8–12], which are protected bymirror symmetry of the crystal [8, 13–15]. These Diracfermions exhibit spin-momentum locking as in TIs; how-ever, there is an even number of Dirac cones at non-TRImomenta, a feature similar to graphene.In general, when Dirac points are located at non-TRImomenta, nonmagnetic perturbations such as strain areable to move Dirac points in momentum space, therebyacting as an effective gauge field on Dirac fermions. Forexample, strain induces opposite gauge fields for the twoDirac valleys in graphene, and spatially inhomogeneousstrain gives rise to effective magnetic fields that are op-posite in two valleys, preserving time reversal symme-try [17–19]. In the presence of such a pseudo-magneticfield B , the low-energy electronic structure takes the formof pseudo-Landau levels (PLLs) with energies charac-teric of Dirac electrons in magnetic fields, i.e. ∼ √ n B ,where n is the Landau level (LL) index. Key signaturesof pseudo-magnetic fields have been experimentally ob-served in graphene [20, 21].The PLLs have a large single-particle degeneracy, which makes them susceptible to many-body instabilitiesin a manner similar to magnetic-field induced LLs [22].Electronic interactions are expected to lift the degener-acy and drive the system into various gapped states. Twoprimary examples are spin-polarized and valley-polarizedstates of PLLs in graphene [23–28].In this work we consider interacting Dirac electronsunder periodically modulated pseudo-magnetic fields,where regions of positive and negative fields alternatein space, forming a superlattice. This field profile leadsto a novel electronic structure markedly different fromuniform pseudo-magnetic fields [29, 30]. There arevarious ways in which periodic pseudo-magnetic fieldscan arise, one prominent way being a strain superlat-tice in graphene or TCI. Such spatially periodic strainfields are particularly relevant, as they were experimen-tally found to develop at interfaces of heterostructuresbuilt from TCIs (e.g., SnTe) and trivial insulators (e.g.,PbTe) [31–33]. At these interfaces, the lattice con-stant mismatch causes dislocations which self-organizeinto a periodic array and therefore produce a naturalrealization of periodic strain fields. The key charac-teristic of the corresponding periodic pseudo-magneticfields is that they can exist over macroscopic regions.In contrast, uniform pseudo-magnetic fields cannot ex-ist in the thermodynamic limit, owing to the bounded-ness of strain (=pseudo-gauge field) [30], unlike a realmagnetic field. Apart from strain superlattices, periodicpseudo-magnetic fields can arise as a result of incom-mensurate electrostatic potentials originating from, forinstance, from lattice-mismatched substrates with a twistangle [34, 35, 41].Starting from the strain-induced pseudo-magnetic su-perlattice, we address the effect of electron-electron in-teractions. Spatially alternating pseudo-magnetic fieldschange the low-energy electronic structure close to theDirac points. Most strikingly, the energy-momentum dis-persion in the vicinity of each Dirac point becomes nearlyflat, leading to a segment of flat band with a twofold de-generacy. These flat bands arise from the zeroth PLLin regions with strong pseudo-magnetic fields; and thetwofold degeneracy corresponds to Landau orbitals that a r X i v : . [ c ond - m a t . s t r- e l ] J un reside in different spatial regions of opposite fields andhave opposite Dirac spinor components [30], a novel fea-ture that is absent in the case of uniform magnetic fields.Counting the two valleys, the flat bands have fourfolddegeneracy. The presence of flat bands leads to a diverg-ing density of states, in contrast to the vanishing densityof states at the Dirac point of massless Dirac fermions.Consequently, periodic strain fields provide a feasible andeffective way of engineering density of states, i.e., elec-tronic compressibility, at zero energy.At charge neutrality, the degenerate flat bands aremainly responsible for driving the spontaneous forma-tion of ordered states. We discuss the various possibilitiesfor degeneracy lifting in the flat band and discriminatebetween energetically favorable and unfavorable states.Two prominent candidate ordered states are the charge-ordered state, where charge is redistributed from the re-gion of positive (negative) to negative (positive) pseudo-magnetic field, and the valley-ordered state, where ineach spatial region the valley degeneracy is lifted. Weshow that the valley-ordered state in graphene and TCIsspontaneously breaks time-reversal symmetry and real-izes an integer quantum Hall effect similar to the Hal-dane state [7], but with the important difference of be-ing driven by electron interactions without any externaltime-reversal-breaking field.We determine the mean-field ground state by self-consistently solving the full gap equation of interactingDirac fermions under a periodically alternating pseudo-magnetic field. The continuum Hamiltonian is micro-scopically implemented using a lattice model with astrain superlattice. Our analysis shows that the supportof the flat band wavefunctions is of great importance.Flat bands in any spatial region only have a twofold val-ley degeneracy, protected by the time reversal symmetry.Therefore lifting this degeneracy by interactions impliestime reversal symmetry breaking. For this reason, we findthe valley-ordered quantum Hall state is greatly favoredover the charged-ordered state under generic forms ofelectron interactions. We show that the order parameterscorresponding to the ordered states follow the strain pro-file, highlighting the crucial role of the pseudo-magneticfield.The present setup for a spontaneous time-reversal sym-metry breaking quantum Hall state relying on a strain-induced flat band should be contrasted with the proposedHaldane mass generation for interacting massless Diracfermions in graphene [43–45]. Exact diagonalization andDMRG studies [48–51] seem to have failed to find theinteraction-driven Haldane phase in models so far pro-posed, in contradiction to Hartree-Fock results. It is be-lieved that the absence of the Haldane phase in ED andDMRG phase diagrams stems from the vanishing densityof states at the Dirac point, and the resulting absence of aweak-coupling instability. In contrast, the quantum Hallstate in our setup already occurs spontaneously at weakcoupling, owing to the strain-induced flat band. II. DIRAC FERMIONS IN TWO DIMENSIONS
We set out to study the coupling of time-reversal in-variant pseudo-gauge fields to Dirac electrons. With twospecific realizations in mind, graphene and TCI surfacestates, we focus on Dirac electrons in two dimensions.Let us start by stating the essential features of Diracelectrons coupled to pseudo-gauge fields, independent ofspecific context.Any two dimensional system respecting time-reversalinvariance and having Dirac fermions not pinned to a par-ticular time-reversal invariant momentum, will consist oftwo species of Dirac fermions. Labeling the two speciesby + and − , the Dirac Hamiltonian describing the twospecies takes the general formˆ H ± = ± (cid:126) v F ˆΨ †± ( − iτ x ∂ x − iτ y ∂ y ) ˆΨ ± (1)where τ i is a set of Pauli matrices acting on the pseu-dospin degree of freedom of the Dirac fermions. Time-reversal symmetry relates the two species by exchangingˆΨ + ↔ ˆΨ − .Pseudo-gauge fields couple to the Dirac fermions in amanner similar to real electromagnetic gauge fields, withone crucial difference, however. In order to respect time-reversal invariance, the pseudo-gauge field must coupleto the fermions in such a way that the two species seeopposite fields. As a result, in the presence of a pseudo-gauge field given by A µ ( µ = x, y ), the Dirac Hamiltonianbecomes ˆ H ± = ± v F ˆΨ †± τ µ (ˆ p µ ± A µ ) ˆΨ ± . (2)This Hamiltonian (with ˆ p µ = − i∂ µ ) describes the genericDirac electrons coupled to pseudo-gauge fields. Pseudo-Landau level quantization will occur when the gauge field A µ acquires spatial dependence, i.e., A µ = A µ ( (cid:126)r ).The interpretation of the Dirac fermion pseudospin andvalley degrees of freedom will depend on the particularrealization of pseudo-magnetic field coupling in a givenmaterial. In this work we will discuss two examples oflow-energy Dirac electrons coupled to time-reversal in-variant gauge fields, which we introduce in the remainderof this section. First we consider the case of graphene,and then we consider surface states of TCIs. Whereasin graphene the Dirac pseudospin degree of freedom de-rives from the two sublattices [6], the pseudospin of theTCI surface state is more complicated due to intrinsicspin-orbit coupling, as we will discuss below.Importantly,in case of the latter, spin-orbit coupling leads to spin-momentum locking in the surface state Dirac theory.In both cases, graphene and TCI, the emphasis will beon strain-induced pseudo-magnetic field coupling. How-ever, we will use the case of graphene to point out thatpseudo-magnetic fields can have a physical origin differ-ent from strain, giving way to an even wider applicationof our results. A. Dirac fermions in graphene
The low-energy theory of graphene at charge neutral-ity is one of the hallmark examples of a 2 D Dirac the-ory [4, 6, 37]. The two species of nodal Dirac fermionsare located at the two inequivalent BZ corners, i.e., theDirac points or valleys, and are labeled by K + and K − corresponding to the momenta (cid:126)K + = (4 π/ ,
0) and (cid:126)K − = − (4 π/ , K and K (cid:48) in small momenta (cid:126)q relative to the Diracpoints [6]. It is given by H ( (cid:126)q ) = (cid:126) v F ν z ( q x τ x + q y τ y ) ≡ (cid:126) v F q µ Γ µ (3)(where v F = √ ta/ (2 (cid:126) )). The set of Pauli matrices τ i acts on the sublattice degree of freedom ( A / B ) and theset of matrices ν i acts on the valley degree of freedom( K + / K − ). In addition, we have defined the Dirac ma-trices Γ x = ν z τ x and Γ y = ν z τ y . The Hamiltonian actson the Dirac spinor ˆΨ( (cid:126)q ) defined byˆΨ( (cid:126)q ) = (cid:104) ˆ ψ A + ( (cid:126)q ) ˆ ψ B + ( (cid:126)q ) ˆ ψ B − ( (cid:126)q ) ˆ ψ A − ( (cid:126)q ) (cid:105) T . (4)Note that we choose the basis so that the A and B latticeare exchanged in the K − valley, meaning that we areworking in the chiral representation (i.e., H ± = ± (cid:126)q · (cid:126)τ forvalley K ± ).Starting from the low-energy Dirac Hamiltonian ofEq. (3), we introduce a generalized time-reversal invari-ant pseudo-gauge field by coupling the Dirac fermions tothe field (cid:126) A i = ( A ix , A iy ) , (5)which consists of three components i = 1 , ,
3. The cou-pling to the fermions has the same form as ordinary min-imal coupling, but with different gauge charges Ω i ex-pressed as H ( (cid:126)q ) = (cid:126) v F Γ µ ( q µ + A iµ Ω i ) . (6)The gauge charge matrices Ω i encode the distinct natureof the pseudo-gauge field as compared to the ordinarygauge field, and are given by Ω i = ( ν x τ z , ν y τ z , ν z ). Thethird gauge charge matrix Ω = ν z is diagonal in valleyspace and assigns opposite sign to the two valleys. There-fore, the component A µ Ω realizes the general Hamilto-nian of Eq. (2) in graphene. In graphene, this is thepseudo-gauge field component that arises in the presenceof strain and plays a central role in this work. The pres-ence of a field (cid:126) A coupling to Ω leads to a moving of theDirac points away from K + and K − , in opposite direc-tions. This is shown in Fig. 1 (left), where the bold bluedots denote the Dirac points moving towards the zonecenter.The following properties of the gauge field charges willbe important for our analysis. The charges Ω i real-ize a pseudospin SU (2) algebra, expressed as [Ω i , Ω j ] = 2 i(cid:15) ijk Ω k . The matrices Ω i commute with the Hamilto-nian in the absence of fields, and as a consequence gen-erate a continuous SU (2) symmetry of the low-energygraphene Hamiltonian. This symmetry is broken whenmass terms are introduced to the Hamiltonian, i.e., whenthe Dirac electrons are gapped out. In particular, theset of mass matrices (cid:126) Γ = (Γ , Γ , Γ ) ≡ ( ν x , ν y , ν z τ z )describes masses that anti-commute with the Hamilto-nian and between themselves. They constitute a set ofcompatible masses, the physical nature of which is well-known. Specifically, the mass Γ = ν z τ z correspondsto an electrostatic potential making the two honeycombsublattices inequivalent and breaking inversion symme-try. Such term is diagonal in valley-space, i.e., it does notcouple the two Dirac points. The other two masses, Γ and Γ , which are off-diagonal in valley space, are knownas Kekul´e masses and correspond to modulations of thetight-binding nearest-neighbor hopping parameter t withtripled unit cell [38]. The breaking of translational in-variance and the modulations over small distances (largemomenta) couple the Dirac points.The gauge charges Ω i act as generators of rotationswithin the space of masses, which follows from the com-mutation relation [Ω i , Γ j ] = 2 i(cid:15) ijk Γ k . In addition to themass terms (cid:126) Γ, there is a mass term τ z , the time-reversalodd Haldane mass [7], which anti-commutes with theHamiltonian (3), but commutes with both the Γ i andthe Ω i . Hence, whereas (cid:126) Γ is a vector under the transfor-mations generated by Ω i , τ z is a scalar.The two remaining gauge charges Ω = ν x τ z andΩ = ν x τ z are off-diagonal in valley space but diagonal insublattice space. The former implies translational sym-metry breaking and the latter implies that these termsarise due to charge density modulations. Consequently,charge density waves (CDWs) with a six-site unit cell,which we will refer to as valley-coupling CDWs, lead toa pseudo-gauge coupling in the same way as strain [34].The SU (2) structure of the gauge charges Ω i implies thatwithin the low-energy theory, the pseudo-gauge field com-ponents are unitarily equivalent to each other. B. TCI surface state Dirac fermions
Topological insulator materials are bulk insulatorshosting gapless Dirac fermions at their surfaces [1, 2].The spin-momentum locked surface Dirac fermions areprotected by time-reversal symmetry, and as a resultthey are pinned at the time-reversal invariant momenta(TRIM). Due to this symmetry-protected pinning, thesurface states of topological insulators do not allow fortime-reversal invariant pseudo-gauge field coupling. Inparticular, strain is not able to move the Kramer’s dou-blet away from the TRIM.In contrast, the TCIs are topological materials pro-tected by crystalline symmetries [8, 16], which host sur-face Dirac fermions not pinned to the TRIM [15, 36, 39].As a result, strain can couple to the low-energy Dirac X X k x k y k k K + K FIG. 1. (Left) Hexagonal Brillouin zone of the honeycomblattice. The two Dirac points K and K are marked by boldblue dots. The blue arrows indicate possible Dirac pointsmoving towards the zone center due to strain, i.e. ∼ u xx − u yy .(Right) Square surface Brillouin of TCI surface state withtwo sets of Dirac points located at X (blue) and X (red).Arrows indicate the moving of Dirac points towards the zonecenter due to symmetric strain ∼ u xx + u yy . fermions as a pseudo-gauge field and can move the Diracpoints in momentum space, in a way that depends on thesymmetry of the strain tensor [15, 30].In this work we specfically focus on the SnTe mate-rial class [8] and its mirror symmetry-protected surfaceDirac fermions appearing on the (001) surface. The sur-face Brillouin zone of the (001) surface is shown in Fig. 1.Two species of low-energy Dirac fermions related by time-reversal symmetry exist in the vicinity of the surfacetime-reversal invariant momenta X and X , representedas blue and red dots in Fig. 1. The surface state DiracHamiltonian at X , given by the terms that respect thecrystal symmetries leaving X invariant, reads H X ( (cid:126)q ) = v q σ y − v q σ x + mν x + δσ x ν y , (7)and a similar expression can be derived for X . Here σ i isa set of Pauli matrices that represents a Kramers doublet,and ν i is a valley degree of freedom corresponding to thetwo inequivalent bulk L -points mapped onto X . Themomentum (cid:126)q is measured with respect to X ; the spin-momentum locking shown in Hamiltonian (7) (i.e., firsttwo terms) comes from spin-orbit coupling. For m = δ = 0 there are two degenerate Kramer’s doublets at X , which are split in energy by finite m and δ . Mostimportantly, finite m and δ leads to the appearance oftwo species low-energy Dirac points, which are locatedat (cid:126) Λ ± = (0 , ±√ m + δ /v ), measured from X .The Hamiltonian of Eq. (8) can be projected into thesubspace corresponding to (cid:126) Λ ± to obtain the effective low-energy Dirac theory. This yields [15] H (cid:126) Λ ± ( (cid:126)q ) = − v (cid:48) q ˜ τ x + v q ˜ τ z , (8)where ˜ τ i is the effective pseudospin degree of freedom, (cid:126)q is now measured with respect to (cid:126) Λ ± , and v (cid:48) = v δ/ √ m + δ . Note that in the chosen basis the Hamil- tonian is valley-isotropic, taking ˜ ν z = ± τ z and ˜ ν z ˜ τ z . The former is a time-reversal even mass andcorresponds to a ferroelectric distortion of the crystal. Itderives from the Dirac bilinear ν z in the basis of (7). Themass term ˜ ν z ˜ τ z breaks time-reversal symmetry and orig-inates from the terms σ z and σ y ν z in the basis of Eq. (7).The mass gap originating from ˜ ν z ˜ τ z is equivalent to thegraphene Haldane gap, and consequently corresponds toa QAH phase [15].Similarly, by using symmetry arguments, the time-reversal invariant pseudo-gauge field couplings can beidentified. As a consequence of the low symmetry ofthe X point, there are no two dimensional representa-tions which directly imply pseudo-gauge coupling. How-ever, since the symmetric terms ν x and σ x ν y displacethe Dirac points in momentum space, any perturbationcoupling to them, will have the effect of a pseudo-gaugefield. Looking for other terms both even under time-reversal and inversion (as expected for strain), one findsanother Dirac bilinear given by σ y ν y . We will show in thenext section, when we discuss strain and strain superlat-tices, that components of the strain tensor couple to theseterms. The effect of these terms is shown schematicallyin Fig. 1 (right), where bold blue and red dots denote theDirac points in the vicinity of X and X , respectively,shifting towards the zone center as a result of strain. III. PERIODIC STRAIN SUPERLATTICES
In the previous section we introduced pseudo-gaugefield coupling in the context of graphene and TCI sur-face states, and argued that strain realizes such coupling.We now turn to a more detailed discussion of strain, andmore specifically periodic pseudo-magnetic fields arisingdue to periodically varying strain, i.e., a strain superlat-tice.Elastic deformations of the crystal lattice, i.e., strainfields, are described by the strain tensor u ij given by u ij = 12 ( ∂ i u j + ∂ j u i ) , (9)where u i ( i = x, y ) is the displacement field. Given thesymmetry of the crystal lattice, the strain tensor canbe decomposed into components transforming as distinctrepresentations of the symmetry group. From this de-composition one can read off which lattice deformationscouple as (pseudo-)gauge fields to the Dirac fermions.In case of hexagonal symmetry, applicable to graphene,there are two d -wave strain components ( u xx − u yy , − u xy ) ∼ ( d x − y , d xy ), which couple to the Diracfermions as the valley-diagonal field (cid:126) A of Eq. (6) (let usomit the label 3, i.e. (cid:126) A = (cid:126) A , for the moment) [40]. Notethat assuming full hexagonal symmetry implies a singlecoupling constant, (cid:32) A x A y (cid:33) ∼ α (cid:32) u xx − u yy − u xy (cid:33) . (10)For illustration purposes, the effect of finite and constant A x is shown graphically in Fig. 1 (left), where the Diracpoints K + and K − move along the k x axis. In case ofthe square symmetry, which applies to the (001) surfacestates of TCIs, the d -wave components d x − x an d x x are not degenerate, and one finds at X [30], (cid:32) A x A y (cid:33) ∼ (cid:32) α ( u xx − u yy ) α u xy (cid:33) . (11)In addition, in the previous section we observed thata perturbation respecting all symmetries can move theDirac points in momentum space, implying that u xx + u yy enters the expression for A x as well, with an independentcoupling. It is the latter strain component, u xx + u yy , theeffect of which is shown in 1 (right).Let us now come to the case of periodic strain fieldswith wavelength λ . More specifically, we consider u ij → u ij ( (cid:126)r ) and consequently (cid:126) A → (cid:126) A ( (cid:126)r ). The periodicity of (cid:126) A ( (cid:126)r ) is directly reflected in the periodicity of the pseudo-magnetic field B = (cid:126) ∇ × (cid:126) A ( (cid:126)r ) = B ( (cid:126)r ), which shouldbe compared to and contrasted with uniform pseudo-magnetic field B . In order to implement the strain super-lattice in a tight-binding setting, we take the graphenelattice as a simple regularization of the continuum the-ory. To solve the superlattice Hamiltonian, we estab-lish a connection between the strain components and thechange in overlap intergals δt n , where n = 1 , , { (cid:126)δ n } . The overlap in-tegral change is expressed in terms of the strain tensor u ij as δt n = (cid:80) n δ in δ jn u ij , which becomes (cid:32) u xx − u yy − u xy (cid:33) ∼ (cid:32) δt − δt − δt √ δt − δt ) (cid:33) . (12)This expresses the pseudo-gauge field in terms of themodulation of the hopping t n → t + δt n .We proceed to consider a single-propagation vectorpseudo-gauge field modulation and obtain the electronicspectrum. A particularly convenient choice is the prop-agation vector δ (cid:126)G ≡ (cid:126)G/λ , where (cid:126)G = (0 , π/ √
3) is areciprocal lattice vector. Then λ is the superlattice wave length, given in terms of graphene unit cells, i.e. λ = 700leads to a superlattice unit cell containing 700 grapheneunit cells. The pseudo-gauge field (cid:126) A and correspondingpseudo-magnetic field are given by A x ( (cid:126)r ) = A cos( δ (cid:126)M · (cid:126)r ) B ( (cid:126)r ) = − ∂ y A x = δM y A sin( δ (cid:126)M · (cid:126)r ) , (13)while A y = 0, since we are only interested in the trans-verse component. A denotes the amplitude of the strainfield, i.e., the maximal change in overlap δt .The spectra in the presence of the strain superlatticefor some values of A are shown in Fig. 2. We observe thatfor increasing A , implying increasing pseudo-magneticfield strength, a flat zero-energy band forms at the Diracpoints, in addition to higher energy dispersive but doublydegenerate bands. This specific reorganization of the low-energy electronic spectrum resembles the Landau levelstructure of external magnetic fields. We will establisha detailed connection between Landau level physics andperiodic strain in the next section. A key feature we wishto stress here, is that the formation of the zero-energy flatband, the degeneracy of which is related to the strengthof the pseudo-magnetic field, leads to a finite and consid-erable density of states (DOS) at the charge neutralitypoint. In stark constrast, in the unstrained case Diracelectrons have linearly vanishing DOS at the charge neu-trality point.Instead of the propagation vector δ (cid:126)G , we can take thepropagation vector δ (cid:126)K = (cid:126)K/λ , where (cid:126)K is the Diracpoint vector defined above. This may be viewed asa simple rotation of δ (cid:126)G , which will result in modifiedMoir´e pattern. We then have for the spatially dependentpseudo-gauge field A x ( (cid:126)r ) = A cos( δ (cid:126)K · (cid:126)r ) B ( (cid:126)r ) = − ∂ y A x = δK y A sin( δ (cid:126)K · (cid:126)r ) (14)where is it important to choose (cid:126)K such that the fieldhas a nonzero transverse component. For given λ thestrain superlattice unit cell constains 3 λ graphene unitcells, which has the benefit that it is commensurate withany perturbation modulated by (cid:126)K coupling the Diracpoints. This allows to treat valley-diagonal and valley-off diagonal perturbations on the same footing.We recall that the low-energy Dirac Hamiltonian in thepresence of strain reads H = (cid:126) v F Γ µ ( − i∂ µ + A µ ( (cid:126)r )Ω ). Asnoted above, a unitary matrix U can be used to rotate toanother gauge field component, U † H U = (cid:126) v F Γ µ ( − i∂ µ + A µ ( (cid:126)r )Ω ). Clearly, this does not change the spectrumand therefore electrostatic potential superlattices, whichwould couple to the gauge field components Ω andΩ [41], are equivalent to strain superlattices. Therefore,even though we focus on strain in this work, we high-light that in the context of graphene spatially modulatedvalley-coupling CDWs induce periodic time-reversal in-variant pseudo-magnetic fields in the same way as strain. -0.10-0.05 0 0.05 0.10 K K K K K K K ( a ) ( b ) ( c ) E / t FIG. 2. Spectra of graphene in the presence of periodically modulated strain for different values of the amplitude of modulation(given by A ∼ δt − δt − δt ; see main text). The modulation wavelength λ is chosen to be 700 graphene unit cells and thepropagation vector is δ (cid:126)G = (0 , π/ √ /λ . The amplitudes are (a) A = 0 . t , (b) A = 0 . t , (c) A = 0 . t . Note that the plotsshow K (cid:48) folded onto the k x -axis. IV. FLAT BAND PSEUDO-LANDAU LEVELSINDUCED BY STRAIN
In this section we address the spectral properties ofDirac electrons in the presence of a time-reversal invari-ant pseudo-magnetic field. As in the previous section,we particularize to the case of graphene (i.e., our lat-tice regularization of the continuum theory), and startby considering a spatially uniform field induced by strain.Uniform fields are fundamentally different from periodi-cally modulated fields induced by the strain superlattice,but we can use the results for the former to develop anintuition for the case of periodic pseudo-magnetic fields.In particular, we may, for the sake of argument, think ofthe periodic field as alternating regions of positive andnegative constant fields, which is schematically shown inFig. 4 (top).
A. Uniform pseudo-magnetic fields and PLLs
The Hamiltonian is given by Eq. (5), for which we onlyretain the strain component (cid:126)
A ≡ (cid:126) A , H ( (cid:126)q ) = (cid:126) v F Γ x ( q x + A x Ω ) + (cid:126) v F Γ y ( q y + A y Ω ) , (15)and (cid:126) A = (cid:126) A ( (cid:126)x ) is taken so as to describe a constant field.The mathematical structure of the Hamiltonian is equiv-alent to that of a real eletromagnetic field and we canuse known techniques to solve it. For completeness webriefly review the essentials here, leaving more details toAppendix A. We first introduce dynamical momenta ˆΠ ± x and ˆΠ ± y for each valley ν = ± , reflecting the fact the signof the pseudo-magnetic field is opposite for the valleys(recall that Ω = ν z ). The Hamiltonian for each of thevalleys reads H ± ( (cid:126)q ) = ± v F (cid:32) ˆΠ ± x + i ˆΠ ± y ˆΠ ± x − i ˆΠ ± y (cid:33) . (16) and the dynamical momenta obey the commutation re-lations [ ˆΠ ± x , ˆΠ ± y ] = ∓ i (cid:126) l B . (17)These commutation relations can be used to define rais-ing and lowering operators in the standard way (seeApp. A), in terms of which the Hamiltonian takes theform H ± = ± (cid:32) √ ξ ˆ a ± √ ξ ˆ a †± (cid:33) . (18)Here we have defined ξ = v F (cid:126) /l b , and the the rais-ing and lowering operators obey the commutation rela-tion [ˆ a ± , ˆ a †± ] = ±
1. This commutation relation is a keyfeature of time-reversal invariant pseudo-magnetic fields,since it reflects anti-parallel field alignment in the twovalleys. The operation of raising and lowering is inter-changed for the two valleys, which has important conse-quences for the structure of the eigenstates. In particular,the eigenstates of the PLL zero modes are localized onthe same sublattice, instead of on opposite sublattices.More specifically one finds | Ψ +0 (cid:105) = (cid:32) | ϕ ,k (cid:105) (cid:33) , | Ψ − (cid:105) = (cid:32) | ϕ ∗ ,k (cid:105) (cid:33) . (19)We stress that this implies localization on the same sub-lattice, given the choice of basis in Eq. (4). Time-reversalsymmetry is preserved by counterpropagation of Landauorbitals in the two valleys (i.e., | ϕ ∗ ,k (cid:105) = | ϕ , − k (cid:105) ). The n = 0 PLL has energy E = 0. Eigenstates correspondingto n (cid:54) = 0 PLLs take the form | Ψ + n ± (cid:105) = 1 √ (cid:32) | ϕ n − ,k (cid:105)±| ϕ n,k (cid:105) (cid:33) , | Ψ − n ± (cid:105) = 1 √ (cid:32) | ϕ ∗ n,k (cid:105)∓| ϕ ∗ n − ,k (cid:105) (cid:33) (20) -0.10-0.05 0 0.05 0.10 200 400 600 800 1000 M/ M/ t / t . . . graphene unit cell graphene unit cell E / t | Bn =0 ,k | | An =0 ,k | FIG. 3. Upper left: Spectrum in the presence of astrain-induced periodic pseudo-magnetic field, shown withboth Dirac points folded onto Γ. Periodicity of the pseudo-magnetic field is 1200 graphene unit cells, and we used A =0 . t . Black arrows indicate for which k the zero mode eigen-states are plotted in the upper and lower right panels. In thesetwo panels we show the wave function distribution ∼ | ψ A,Bn =0 ,k | of the full zero mode (i.e. n = 0) subspace over the grapheneunit cells for the A -sublattice (lower right, red) and B -sublattice (upper right, black). Black arrows indicate which k they correspond to in the upper left plot; the blue arrow in-dicates in which order. Lower left: plot of the periodic strainmodulation A x = A cos(2 πy/λ ) ∼ δt ( y ) − δt ( y ) − δt ( y )(red) and corresponding pseudo-magnetic field. Note that forclarity we have rescaled the amplitude of the pseudo-magneticfield to A . and they have energies E ± ( n ) = ± (cid:112) ξ n for each valley.In order to gain insight into the effect of perturba-tions on periodic strain superlattices, we review the ef-fect of such perturbations on the PLL spectra for uniformpseudo-magnetic field we just presented. Consider the n = 0 PLL, the eigenstates of which are given in Eq. (19).Let us first comment on the mass gap terms ν z τ z and τ z . The effect of the sublattice CDW ν z τ z and the TRSbreaking Haldane term τ z is reversed as compared to realmagnetic fields [23]. The sublattce polarized CDW sim-ply shifts the energy of zero modes but does not breaktheir degeneracy, whereas the Haldane mass τ z energet-ically splits the zero modes in a symmetric way. Thisis an immediate consequence of the sublattice structureof the PLL zero modes. The two Kekul´e masses ν x and ν y do not affect the zero modes at all, they are neithersplit nor shifted, since they are off-diagonal is sublatticespace.Perturbations that do split the zero modes in a fash-ion similar to the Haldane term are charge density waveswith tripled unit cell, i.e., charge density waves that cou-ple the valleys K + and K − [34]. These charge densitywaves couple to the Dirac mastrices ν τ , ν τ , ν τ and ν τ . Projecting these into the PLL zero mode subspace,one finds effective Pauli matrices ˜ τ x and ˜ τ y , which anti-commute with the Haldane mass projected into the zero mode space, τ z → ˜ τ z . This leads to the counterintuitivesituation of anticommuting masses only one of which isTRS breaking and nontrivial [42]. We note that thesecharge density waves with tripled unit cell correspondto the other gauge field components of the SU (2) gaugefield, i.e. they enter as A x , A y , A x , and A y in Eq. (6).Yet another perturbation that splits the zeroth PLL isthe valley mass, given by ν z , making the valleys inequiv-alent, but acting as the identity in sublattice space. Itsspectral effect is equivalent to that of the Haldane term,meaning a symmetric splitting of the zero modes.Understanding the spectral effect of these Diracfermion bilinears on the zeroth PLL gives a first ideaof the ways in which their spontaneous formation canlower the energy for charge neutral systems. For a morerefined understanding of the energetics it is necessary toconsider the effect of perturbations on higher PLLs (i.e., n (cid:54) = 0). For both the CDW term ν z τ z and the Hal-dane term τ z all PLLs with n (cid:54) = 0 get pushed up ordown in energy depending on whether they have energies ± (cid:112) ξ n , i.e. positive (negative) solutions get pushed up(down). This is different for the valley mass ν z , whichpushes all PLLs of valley K + up and of valley K − down,effectively splitting all PLLs, even the n (cid:54) = 0 levels, leav-ing do degeneracies behind. The charge density waveswith enlarged unit cell (i.e., coupling the valleys) bothsplit and shift the higher order PLLs, which may be seenstraightforwardly by using perturbation theory up to sec-ond order. Based on these considerations we obtain anintuition for the spontaneous generation of Dirac fermionbilinears due to interactions, depending on the locationof the Fermi level. B. Alternating pseudo-magnetic fields andsuperlattice PLLs
The assumption of uniform pseudo-magnetic field is auseful first step towards understanding the physics of su-perlattice PLLs. As a next step we consider the caseof alternating pseudo-magnetic field, which for simplic-ity we will take as a periodic arrangement of regions ofpositive and negative constant field. This will providevaluable insight in the case of periodic harmonic pseudo-magnetic fields. Schematically this field arrangement isshown in the top row of Fig. 4. Figure 4 shows how wecan think of the a periodically alternating field as a strainsuperlattice with effective “two-site” unit cell (i.e., posi-tive and negative field), reminiscent of an antiferromag-net, leading to a doubling of the PLL degeneracies. Forinstance, the space of zero mode PLLs is doubled, sincewe have the spatial degeneracy in addition to valley de-generacy. For each valley there is a zero mode localizedin the positive field region, meaning on the A sublattice,and a zero mode localized in the negative field region, onthe B sublattice.The additional degree of freedom originating from theperiodicity of the pseudo-magnetic field gives rise to a + B B + B B + B B K K AA B K K B K K AA B K K B K K AA B K K B K K AA B K K B K K AA B K K B K K AA B K K B K K AA K K AA K K AA B K B K B K K B K B K B FIG. 4. Schematic representation of the PLL structure incase of spatially alternating pseudo-magnetic field. The al-ternating regions of positive and negative field are depictedin the upper row. In the row below we depict the degener-ate n = 0 PLLs in the two regions, the structure of whichperiodically alternates in accordance with the field. Here A and B label a general pseudo-spin degree of freedom, whichcorresponds to the A/B sublattice in graphene. The bottomtwo rows show the two prime polarized state candidates. Inthe charge ordered state with ferroelectric polarization, theenergy levels in the positive field ( A ) have higher energy thanthe negative field ( B ) states, the latter being fully occupied.In the valley-ordered quantum Hall state (very bottom) lev-els are split in each region, occupying a single valley in eachregion. richer structure of polarized or ordered states. Focus-ing on the PLL zero mode subspace, relevant at chargeneutrality, there are multiple ordered states that lift thedegeneracy of the zero mode subspace. Two of them areshown in Fig. 4. The first is a charge ordered state, wherethe zero mode PLLs in one of the two spatial regions areboth occupied, leaving the zero modes in the other re-gion unoccupied. This leads to a redistribution of chargebetween the two regions and an associated ferroelectricpolarization along the propagation direction of the su-perlattice wave vector. In case of graphene this stateis realized by the sublattice CDW, which energeticallydiscriminates the sublattices. The other state shown inFig. 4 is the valley-ordered quantum Hall (or Haldane)state. In such state the zero modes corresponding toan “up” pseudo-magnetic field are occupied. Note thatthis implies an alternating occupation of valleys K ± , asshown in Fig. 4. Therefore, this state may be called anti-ferro-valley -ordered. In graphene the valley-orderedquantum Hall state is realized by the time-reversal sym-metry breaking Haldane term. A third PLL polarizedstate is obtained by occupying the same valley in eachspatial region. This state also breaks time-reversal, butcontrary to the valley-ordered quantum Hall state thepseudo-magnetic field seen by the occupied PLLs alter- nates. The inversion of PLL occupation in one of the tworegions with respect to the anti-ferro-valley-ordered statesuggest the name ferro-valley -ordered. In graphene suchstate is realized by the valley mass term.We now establish a connection between the simplifieddescription of alternating pseudo-magnetic fields in termsof continuum PLLs, and the periodic strain (super)latticemodel introduced in the previous section. In order todo so we take the unidirectional periodic strain profilecompatible with an enlarged six-site unit cell (i.e., Diracvalleys folded onto Γ) defined in Eq. (14). Solving thetight-binding Hamiltonian in the presence of the strainsuperlattice yields the spectrum shown in Fig. 3 (upperleft). The connection is made by interpreting this spec-trum in terms of PLLs.Let us first study the wave function support of thezero energy solutions and compare that to the zerothPLL. The right column of Fig. 3 shows the wave func-tion support | ψ A,Bn =0 ,k | of the wave functions ψ n =0 ,k cor-responding to zero energy solutions ( n = 0) labeled by k ( k should be identified with k y ). Black arrows ex-plicitly indicate which k corresponds to which | ψ A,Bn =0 ,k | profile. Since there are two valleys and the strain su-perlattice unit cell consists of two distinct regions withopposite pseudo-magnetic field, there is a fourfold de-generacy at each k . This is reflected in Fig. 3 where twopseudo-magnetic Landau-like orbitals are localized on the A -sublattice (lower right) and two localized on the B -sublattice. The wave function support clearly shows thespatial separation of solutions living on distinct sublat-tices. In the region of positive field (see Fig. 3 lower left)“zero modes” are localized on the A -sublattice, and onthe B -sublattice in the region of negative field. In addi-tion, we observe that for k moving away from Γ (followingthe blue arrow) the two Landau orbitals in each regionmove away from the position of maximum field towardsthe position of vanishing field. In particular, they moveaway in opposite directions, which is a direct consequenceof their different valley index. As the two valleys effectivesee opposite fields, and the spatial position of a Landauorbital is given by x ∝ sgn( B ) k y , the Landau orbitals areexpected to spatially move in opposite directions. Withincreasing k ( ∼ k y ) Landau orbitals of different spatialregions and the same valley start to overlap, eventuallyleading to the splitting observed in spectrum indicatedby the most right black arrow in the upper left panel ofFig. 3. Hence, at the junctions between regions of pos-itive and negative field, the Landau orbitals acquire adispersion and form a series of snake states [30].Next, we ask what the spectral effect is of the perturba-tions that are expected to split degeneracies, in particularthe flat band degeneracies as depicted in Fig. 4. To thisend, we solve the strain superlattice Hamiltonian in thepresence of various perturbations, which for the momentwe take to be spatially uniform, i.e., not follow the super-lattice envelope, in accordance with the schematic pictureof alternating regions of constant field (Fig. 4). Figure 5shows the low-energy spectrum in the vicinity of Γ. The -0.10-0.05 0 0.05 0.10-0.10-0.05 0 0.05 0.10 ⇡ ⇡ ak x ⇡ ⇡ ak x ak x ( a ) ( b ) ( c )( d ) ( e ) ( f ) E / t E / t FIG. 5. Spectra of graphene in the presence of periodic strain and in the presence of additional perturbations splitting or liftingPLL energies. Spectra are obtained for strain superlattice unit cells containing λ = 600 graphene unit cells and A = 0 . t . (a)Free graphene (b) CDW (c) Haldane term (d) valley mass term (e) CDW with tripled unit cell (f) CDW with tripled unitcell. unperturbed case corresponding to ˆ H + ˆ H strain is shownin Fig. 5(a) for reference. The sublattice charge orderedstate and (anti-ferro-))valley-ordered quantum Hall stateare shown in Fig. 5(b) and Fig. 5(c), respectively. Bothopen up a full gap, splitting the zero mode subspace andshifting the higher PLLs, as expected. Figure 5(d) showsthe effect of a valley mass term, i.e., the ferro-valley-ordered, on the low-energy spectrum. Whereas in eachvalley degeneracies are preserved, the valleys are split, asexpected. As a consequence, the spectrum is not gappedand the Fermi level crosses the propagating snake statesassociated with the flat band PLLs [30].Figure 5(e) and 5(f) show the spectrum obtained inthe presence of valley-coupling CDWs, which we havediscussed split the continuum n = 0 PLL in a way sim-ilar to the Haldane term. We observe that in case of astrain superlattice, the CDWs do not lead to a full gap,but only lift the degeneracy of the flat band states local-ized at the position where the pseudo-magnetic field hasits extrema. The degeneracy is not lifted in the vicin-ity of the nodes of the periodic pseudo-magnetic field.The absence of a full gap in case of periodic strain canbe understood by considering the spectral effect of thevalley-coupling CDWs in case of zero pseudo-magneticfield. The valley-coupling CDWs do not open up a gapin that case, but only shift the Dirac points. Hence, atthe nodes of the pseudo-magnetic field one expects theabsence of a gap. Note also that the spectrum numeri-cally obtained numerically shows both split and shiftedhigher ( n (cid:54) = 0) PLLs. From this analysis we conclude that in the presence ofthe strain superlattice, the low-energy electronic struc-ture can be approximated by sets of PLLs for the twodistinct regions of the superlattice unit cell. In addi-tion, based on their effect on the degenerate low-energyPLLs, we expect the charge-ordered and anti-ferro-valley-ordered states to be the dominant instabilities in the pres-ence of the pseudo-magnetic field superlattice.We end this section with two remarks. First, we notethat the analysis presented here is based on the the as-sumption of a strain-induced pseudo-magnetic field, i.e., (cid:126) A = (cid:126) A in Eq. (6). As mentioned in Section II, a unitarymatrix can rotate to another component, e.g. (cid:126) A or (cid:126) A .This is does not change the (low-energy) spectrum, butit does change the nature of the eigenstates. In addition,it also changes the nature of the perturbations, since theunitary matrix rotates within the space of masses rep-resented by (cid:126) Γ as well. In particular, this implies that asublattice polarized term ν z τ z will be rotated into oneof the Kekule terms. Interestingly, the Haldane term isa scalar under these unitary rotations and hence is in-variant. To summarize, the analysis of this section stillapplies, but in a rotated basis.The second remark concerns the applicability of ouranalysis to TCI surface states. The arguments put for-ward in the present section build on the specific exampleof graphene PLL physics. They remain valid in the con-text of TCI surface states. Most importantly, in the pres-ence of a uniform pseudo-magnetic field, the TCI surface n = 0 PLLs are localized on the TCI pseudospin de-0gree of freedom in such a way that the time-reversal in-variant ferroelectric distortion of the crystal lattice onlyshifts them, whereas a time-reversal breaking Zeeman-type spin-coupling splits them in energy. As a result,there is a one-to-one correspondence between the valley-ordered quantum Hall state and charge-ordered state ingraphene, and Zeeman term and ferroelectric distortionof TCI surface states. V. INTERACTING ELECTRONS IN A STRAINSUPERLATTICE
In order to systematically investigate the patterns ofsymmetry breaking and PLL splitting, as a consequenceof interacting flat band electron, we have studied an in-teracting Hamiltonian on the graphene honeycomb lat-tice and performed extensive self-consistent Hartree Fockcalculations. We report the results in this section.Based on the intuitive picture of the PLLs in the twospatially separated regions we anticipate both the forma-tion of a charge-ordered state with ferroelectric polariza-tion, corresponding to a redistribution of charge betweenthe regions of positive and negative pseudo-magneticfield, and the formation of a valley-ordered quantum Hallground state. Interactions which may lead to the forma-tion of these states in a graphene lattice model are thenearest neighbor (NN) and next-nearest neighbor (NNN)density-density interactions [43], respectively, as will bedemonstrated below. We therefore consider the interact-ing Hamiltonian ˆ H = ˆ H + ˆ H (cid:126) A + ˆ H int where the inter-acting part ˆ H int is given byˆ H int = V (cid:88) (cid:104) rr (cid:48) (cid:105) ˆ n r ˆ n r (cid:48) + V (cid:88) (cid:104)(cid:104) rr (cid:48) (cid:105)(cid:105) ˆ n r ˆ n r (cid:48) . (21)Here V is the NN interaction strength and V the NNNinteraction strength, and ˆ n r is the number operator atsite r . The sums over (cid:104) rr (cid:48) (cid:105) and (cid:104)(cid:104) rr (cid:48) (cid:105)(cid:105) are over NN andNNN, respectively.The spatially modulated strain is implemented in theway described above, meaning that we take ˆ H (cid:126) A to con-tain (cid:126) A = ( A x , A x component of the strain-induced gauge field originates from hopping amplitudemodulations A x ∼ δt − δt − δt , which are given aspatial dependence δt n → δt n ( (cid:126)r ). A schematic repre-sentation of the way we set up the calculations is givenin Fig. 6, which is a generalization of the approach ofRef. 44 to the case of strain superlattices (more detailsin Appendix B). In order to allow for charge densitywaves that couple the Dirac points K + and K − we workwith a tripled (6-site) unit cell (red dashed hexagonsin Fig. 6) [44]. The corresponding lattice vectors are (cid:126)b = 2 (cid:126)a + (cid:126)a and (cid:126)b = (cid:126)a − (cid:126)a in terms of the graphenelattice vectors (cid:126)a i . We take the strain-induced gauge fieldto be A x ( (cid:126)r ) = A cos( (cid:126)K (cid:48) · (cid:126)rλ/ k x k y M K K M A B B B A A n nn + 1 n + 1 t t t ~a ~a ~b ~b K + K FIG. 6. Upper left: Schematic representation of the grapheneBrillouin zone (outer black hexagon) and the folded Bril-louin zone corresponding to the tripled unit cell (inner blackhexagon). Dirac points of pristine graphene are located at K and K (cid:48) , which are folded onto Γ when tripling the unitcell. The vectors connecting Γ to K and K (cid:48) , indicated byblue arrows, are reciprocal lattice vectors of the enlarged lat-tice vectors. Lower left: Wigner-Seitz sells containing threeelementary graphene unit cells, with lattice vectors (cid:126)a and (cid:126)a . Right: graphene lattice (black hexagons) and tripled unitcells (red dashed hexagons), including labeling of sites. Ai and Bi label the six sites within the enlarged unit cell; n and n +1 label the position in the strain-induced superlattice. δt n are the hopping amplitude changes due to strain. where (cid:126)r denotes the position of an elementary grapheneunit cell, (cid:126)K is the K − wave vector as shown in Fig. 6and λ equals the number of graphene unit cells in theperiodic strain superlattice unit cell. The number of 6-site unit cells in the superlattice unit cell is therefore λ/ (cid:126)b and λ(cid:126)b .We stress that we consider a fully periodic system withwavelength ∼ λ set by the periodic strain superlattice.Spinless electrons on the graphene honeycomb latticeserve as a model system for strain-induced PLLs. Wetherefore first present Hartree-Fock results for spinlesselectrons and then briefly comment on the spin degree offreedom, recalling that the case of TCI surface states issimilar to spinless graphene electrons. A. Spinless electrons
Within the framework of standard Hartree-Fock theorywe decouple the interactions both in diagonal (Hartree)and off-diagonal (Fock) channels. The diagonal (orcharge density) order parameter at a given site in thesuperlattice unit cell is labeled by ( α, i, l ), and is definedas ρ αi ( l ) = 1 N (cid:88) (cid:126)k (cid:104) ˆ ψ † αi ( l, (cid:126)k ) ˆ ψ αi ( l, (cid:126)k ) (cid:105) , (23)1where α = A, B ; i = 1 , , l labels the 6-site cell inthe superlattice cell, and N is the total number of super-lattice unit cells. We define ρ α ( l ) as the average over i within cell l , ρ α ( l ) = 13 (cid:88) i =1 ρ αi ( l ) . (24)Based on previous considerations we are interested inpossible local sublattice imbalances ∆ ρ − ( l ) expressed as∆ ρ − ( l ) = ρ A ( l ) − ρ B ( l ) , (25)and possible ferroelectric redistribution of charges be-tween regions of positive and negative pseudomagneticfield, which can be expressed as∆ ρ + ( l ) = ρ A ( l ) + ρ B ( l ) − . (26)In addition to the charge density order parameter,in order to detect the valley-ordered Haldane state, westudy the quantum Hall (QH) order parameter originat-ing from spontaneously generated next-nearest neighbortunneling with complex amplitude. Decoupling in theoff-diagonal Fock channel leads to χ αil,jl (cid:48) = 1 N (cid:88) (cid:126)k Q α ∗ il,jl (cid:48) ( (cid:126)k ) (cid:104) ˆ ψ † αi ( l, (cid:126)k ) ˆ ψ αj ( l (cid:48) , (cid:126)k ) (cid:105) , (27)where the phase factors Q αil,jl (cid:48) ( (cid:126)k ) are defined in Ap-pendix B, and from these we extract the QH order pa-rameter ∆ QH ( l ) as∆ QH ( l ) = (cid:88) ijl (cid:48) ∈ NNN( l ) η α Im χ αil,jl (cid:48) . (28)The sum is over all NNN pairs which belong to cell l and η A = − η B = 1 since in the QH state fluxes on the A and B sublattices have opposite sign.The results of self-consistent HF calculations wepresent here were performed for spinless electrons on lat-tices of size N = 16 ×
16 and λ = 3 ×
40 = 120. In thefollowing we map out the phase diagram as function of V and V by discussing the results for three specific regimesseparately. First we will focus on ( V = 0 , V (cid:54) = 0) toshow that the QH state is stabilized for the smallest val-ues of V as a result of the strain induced flat PLL-likebands. Second, we will look at the case ( V (cid:54) = 0 , V = 0)to show that a NN interaction will induce the ferroelec-tic charge ordered state. Third and last we will focuson selected cases of ( V (cid:54) = 0 , V (cid:54) = 0) to show that theferroelectric charge density wave is strongly suppressedas compared to the valley-ordered QH state, preciselydue to the localization of the low-energy states flat bandstates on a single sublattice.Figure 7 shows the HF results for various values ofNNN interaction V while keeping V = 0. Panels (b)-(c)show the QH order parameter ∆ QH ( l ) defined in Eq. (28)as function of then tripled cell index l . Black and red -0.25-0.20-0.15-0.15-0.05 0 -0.25-0.20-0.15-0.15-0.058 16 24 32-0.25-0.20-0.15-0.15-0.05 0 8 16 24 32 -0.30-0.25-0.20-0.15-0.10-0.05 0 0.2 0.4 0.6 0.8 V /t tripled unit cell tripled unit cell ( a ) ( b )( c )( d ) Q H Q H l l FIG. 7. Panel (a) shows the maxima of QH order param-eter ∆ QH as function of NNN interaction V for various val-ues of the pseudo-magnetic field amplitude: A = 0 .
10 (red), A = 0 .
15 (blue), A = 0 .
20 (black). Panels (b)-(c) show theQH order parameter as function of tripled unit cell ∆ QH ( l )( l labeling tripled unit cell) for the same values of A [ (b) A = 0 .
10; (c) A = 0 .
15; (d) A = 0 . V = 0 . , . , . , . , . , . , . , . , . , .
0, in descend-ing order, i.e. bottom most curves V = 1 . V = 0 .
1. Black and red curves correspond to A and B sublattice, respectively. V = 0 in all cases. curves correspond to A and B sublattice, respectively,and the strength of the tunneling amplitude variation A is A = 0 .
10 (b), A = 0 .
15 (c), and A = 0 .
20 (c). It isapparent from these figures that the QH order param-eter follows the profile of the effective pseudo-magneticfield. In the spatial region where the flat band statesare localized on the A sublattice, the QH order param-eter develops predominantly on the A sublattice, andvice versa for the B sublattice region. Ordinarily, in thehoneycomb lattice QH state, the spontaneously inducedmagnetic fluxes are opposite on the A and B sublattice,averaging to zero over an elementary unit cell [7]. Inthe present case the localization of flat band states on asingle sublattice leads to finite fluxes in regions of pos-itive and negative pseudo-magnetic field, which averageto zero only over the larger strain superlattice unit cell.In addition to the locking of the QH order parameter∆ QH ( l ) to the sublattice structure of the flat band zeromodes, we observe that the strain-induced reorganizationof the low-energy electronic structure into PLLs, funda-mentally changes the impact of interactions. Figure 7(a)shows the dependence of the QH order ∆ QH ( l max ), where l max is the cell index where the QH order is strongest, onthe strength of the interaction V for various values ofthe pseudo-magnetic field strength (i.e., A ). Whereasfor the unstrained honeycomb lattice a finite interaction V ∼ . tripled unit cell tripled unit cell ⇢ A + ⇢ B ⇢ A ⇢ B / ( a ) ( b ) l l FIG. 8. (a) Plots of the charge density wave order pa-rameter ρ A ( l ) − ρ B ( l ) (red) and total charge redistribution ρ A ( l ) + ρ B ( l ) − V = 0 . , . , . , . , . , .
8, in both cases plotted in ascend-ing order (top curve V = 0 .
8, bottom V = 0 . A = 0 . A = 0 . bands such as PLLs one expects the interaction-inducedorder to scale linearly with interaction for weak coupling[46, 47]. This is reflected in Fig. 7(a) which shows that forincreasing pseudo-magnetic field, setting the PLL bandflatness, the QH is more robust and its dependence onthe interaction is approximately linear, with deviationsat very small values of V .Next, we turn to the case of finite NN interaction V while keeping V = 0. In the absence of strain the unfrus-trated NN interaction will favor a CDW characterizedby translational symmetry preserving sublattice chargeimbalance [43–45]. On the contrary, in the presence ofstrain, adopting the PLL picture for the low-energy elec-trons and focusing on the zeroth PLL, the effect of theNN interaction is expected to be suppressed, as the ze-roth PLL states live exclusively on one sublattice. Nev-ertheless, since the NN has the potential to cause chargeasymmetry between the sublattices (i.e., the regions ofpositive and negative pseudo-magnetic field), and higherPLLs may be relevant to the energetics, we anticipatethe system to develop a charge density wave of ferroelec-tric type (i.e., ∆ ρ + ), with excess charge in regions wherethe pseudo-magnetic field is positive, and defect chargewhere it is negatice, or vice versa. In addition, in theprevious section we observed that under the assumptionof a uniform CDW the strain-induced PLL spectrum isgapped out.In Fig. 8 we show results for both ∆ ρ + ( l ) and ∆ ρ − ( l ),defined in Eqs. (26) and (25), for different values of the in-teraction V and strain A . As expected, we find the CDWorder parameter ∆ ρ − ( l ) (shown in red) to become finitein regions where the pseudo-magnetic field is strongest,but to have the same sign in both positive and negativeregions. The sign of concomittant ferroelectric polariza-tion depends on the sign of the pseudo-magnetic field. Inthe region where the field is positive the flat band statesare localized on the A sublattice and pushing them downin energy, signaled by positive ∆ ρ − ( l ), leads to excesscharge in that region at the expense of charge in the re-gion of negative field. At the same time we observe thatthe FP is more pronounced for stronger strains, and that -0.14-0.12-0.10-0.08-0.06-0.04-0.02 0 0.02 8 16 24 32 tripled unit cell ⇢ A + ⇢ B / V = 1 . V = 1 . V = 0 . V = 0 . V = 0 . V = 1 . V = 1 . l Q H FIG. 9. Plot of both the QH order parameter ∆ QH ( l ) (only A -sublattice shown; red dot curves) and total charge redis-tribution ρ A ( l ) + ρ B ( l ) − n for A = 0 .
20 and V = 0 .
20. Differentcurves correpond to different values of V , explicitly labeledfor clarity. it is very weak for small interaction. The latter may beattributed to the fact that NN interactions have no effectin the zeroth PLL.We proceed to consider the case of both finite V and V . We have seen that both of these interaction individ-ually favor different gapped ground states. At the sametime we argued that as a consequence of the different na-ture of these interactions, i.e. V being inter-sublatticeand V being intra-sublattice, they have different impacton the low-energy flat band electrons. Due to the spatialseparation of states localized on different sublattices, it isexpected that the effect of V is suppressed. One there-fore expects the QH state to survive even for V < V ,which corresponds to the physically relevant regime.In Fig. 9 we present results for various V with finite V = 0 . A = 0 .
20. We plot both the ferroelec-tric order parameter ∆ ρ + ( l ) and the QH order param-eter ∆ QH ( l ) for the A sublattice. The key observationis that even for small V = 0 . V ∼ .
0. The therefore conclude thatthe effect of interactions on periodic strain induced flatbands follows from their PLL character. In particular,the sublattice structure of the low-energy flat bands isthe decisive factor in determining which order is sponta-neously generated by interactions.
B. Remarks on spinful electrons
We close this section with a number of remarks on theelectron spin. The numerical calculations have been per-formed for spinless electrons in graphene. Taking spininto account gives rise to a richer structure of polarizedstates, specifically the ferromagnetically (FM) and anti-ferromagnetically (AFM) polarized states should then beconsidered, in addition to the Quantum Spin Hall polar-ized state. Moreover, the argument for the suppression ofthe NN interactions does not apply to the onsite Hubbard3interaction, ˆ H U = U (cid:80) r ˆ n ↑ r ˆ n ↓ r , which must be includedin the interacting Hamiltonian.The results for spinless electrons in graphene do, how-ever, directly apply to TCI surface, which do not have anadditional degenerate spin degree of freedom. Instead,as a consequence of spin-orbit coupling, spin is alreadypart of the low-energy Dirac structure. In particular,our results imply that on the surface of a TCI and in thepresence of periodic strain, interactions will lead to theformation of the QH state.As a first step towards understanding the polarizationof spinful strain superlattice-induced PLLs in graphene,we have performed numerical Hartree-Fock calculationswith an interacting Hamiltonian given by H U . Wefind that the mean-field ground is a superlattice anti-ferromagnet . The superlattice anti-ferromagnet exhibitsanti-ferromagnetic order, as expected on a bipartite hon-eycomb lattice, however, since the flat band states arelocalized on one sublattice only, in each of the two re-gions of the strain superlattice an effective magnetiza-tion develops. Hence, as a consequence of the particularstructure of the zeroth PLL states, the anti-ferromagneticorder is transferred to the superlattice.A similar result was reported in Ref. 28, which foundanti-ferromagnetic order induced by non-periodic strain,where the bulk and the sample boundary have an effectivebut opposite magnetization. VI. DISCUSSION AND CONCLUSION
We have shown that in the presence of a strain su-perlattice, a periodic modulation of elastic lattice defor-mations, a system of low-energy Dirac electrons exhibitsa fourfold degenerate zero energy flat band, reminiscentof a zeroth PLL. The PLL structure originates from thepseudo-magnetic field, generated by nonuniform strain.The strain superlattice unit cell consists of two spatiallydistinct regions, one in which electrons see a positivepseudo-magnetic field and one in which they see a nega-tive field. The single-particle states of the degenerate flatband have a special and important localization property:in each of the two regions their wavefunction has supportonly on one of the Dirac pseudospin species.Periodic pseudo-magnetic fields can occur both ingraphene and on the surface of a TCI, which hostspairs of Dirac fermions at opposite momenta relatedby time-reversal symmetry. The important fact thatDirac fermions are unpinned to time-reversal-invariantmomenta in the BZ allow for pseudo-gauge field undertime-reversal invariant perturbations such as strain.Interactions between electrons in continuum PLLs areexpected to lead to the formation of polarized states split-ting the degeneracies of PLLs. We have investigated PLLpolarization for the case of lattice PLLs correspondingto periodic pseudo-magnetic fields. Two polarized stateswere shown to fully lift the zero energy flat band de-generacy while at the same time pushing all occupied (unoccupied) lattice PLLs down (up). The first state isthe sublattice polarized charge-ordered state, which canbe pictured as a spatially polarized state with all PLLsof the positive (or negative) field region occupied. Thesecond is the anti-ferro-valley ordered state, or sponta-neous quantum Hall state, for which all PLL states effec-tively seeing a positive (or negative) field are occupied,implying time-reversal symmetry breaking. We foundthat other polarized states, even though they have simi-lar characteristics in the continuum, do not fully lift thelattice flat band PLL degeneracies.Using self-consistent Hartree-Fock calculations we havestudied an interacting honeycom lattice model with pe-riodic strain-induced lattice PLLs. Our results demon-strate that the strain-induced reconstruction of low-energy elecronic structure, in particular the presence of azero energy flat band, determines the impact of interac-tions. Three key results highlight this conclusion. First,the mean-field order parameters clearly reflect the pe-riodicity of the pseudo-magnetic field, showing that theamplitude of the order parameter is tied to the strengthof pseudo-magnetic field.Second, as a consequence of the characteristic wave-function support of the flat band single-particle states,the effect of interactions that favor time-reversal invari-ant pseudo-spin order is suppressed. This effectively en-hances interactions that favor time-reversal symmetrybreaking valley order with associated spontaneous quan-tum Hall effect. We have established this result in thecontext of the graphene lattice model, where the pseudo-spin corresponds to the sublattice degree of freedom.The NN interactions are inter-sublattice interactions, andtherefore supressed, as the flat band single-particle stateshave support on one of the sublattices only, in each pos-itive or negative field region. The result, however, isgeneral and applies equally to TCI surface states. Thephysical interpretation of pseudo-spin and valley are dif-ferent, as explained in Sec. II, yet the effect of interac-tions favoring time-reversal symmetry breaking remainsstrongly enhanced.Third, the valley-ordered spontaneous quantum Hallstate already occurs for small interactions when the PLLsare well-developed. The PLL structure is a way to signif-icantly enhance density of states near the charge neutralpoint. We conclude that in the presence of periodic strainand interactions, a system of unpinned Dirac electronshas a generic instability towards a spontaneous quantumHall phase.The emphasis of the numerical calculations we report,has been on spinless electrons in graphene. In TCI sur-face states, however, no additional spin degeneracy ispresent. Due to spin-orbit coupling the electron spin isan intrisinc part of the low-energy Dirac structure. Inparticular, this implies that there are no purely spin-polarized phases, such as the global anti-ferromagnet,competing with the spontaneous quantum Hall phase,favoring the latter as ground state.Whereas uniform strain-induced pseudo-magnetic4fields suffer from implementation limitations, particu-larly beyond the nano-scale, periodic strain can poten-tially be realized in macroscopic sample sizes. In fact,such periodic strain fields and induced pseudo-magneticfields were demonstrated in TCI heterostructure [32, 33],making TCI surface states the prime candidate to exhibitspontaneous formation of nontrivial electronic states.
ACKNOWLEDGMENTS
We thank Evelyn Tang, Vlad Kozii and SarangGopalakrishnan for interesting discussions. J.V. ac-knowledges support from the Netherlands Organizationfor Scientific Research (NWO). L.F. is supported byDavid and Lucile Packard Foundation.
Appendix A: Landau levels for Dirac electrons
For the purpose of being selfcontained we collect somestandard results of Dirac fermions in a constant magneticfield in this Appendix. These may be directly applied tothe case of time-reversal invariant pseudo-magnetic fields,bearing in mind the key characteristic of opposite sign ofthe pseudo-magnetic field in the two valleys.In the presence of a magnetic field we define the dy-namical momenta using the Peierls substitutionˆ p α → ˆΠ α = ˆ p α − eA α (ˆ r ) = − i (cid:126) ∂ α + | e | A α (ˆ r ) , where ˆ p α is the momentum operator and α = x, y (werestrict the description to two dimensions). ˆ r α are theposition operators obeying [ˆ r α , ˆ p β ] = i (cid:126) and A α (ˆ r ) isthe electromagnetic gauge field. In a magnetic field themomentum operators ˆΠ α do not commute but insteadobey the canonical commutation relation[ ˆΠ α , ˆΠ β ] = | e | [ˆ p α , A β (ˆ r )] + | e | [ A α (ˆ r ) , ˆ p β ]= − i (cid:126) | e | F αβ where F αβ = ∂ α A β − ∂ β A α is the field strength. Assum-ing a uniform field strength in the ˆ z direction the mag-netic field is given by B λ = (cid:15) λµν F µν /
2, implying that F µν = (cid:15) µνλ B λ . In particular we have for a uniform field B ≡ B z in the ˆ z direction[ ˆΠ α , ˆΠ β ] = − i (cid:126) | e | sgn( B ) B(cid:15) αβz . (A1)In this expression we have explicitly separated the signof the magnetic field from its strength B = | B | so as tomake dependencies on the sign of the field transparent.Defining the fundamental characteristic length scale inthe system, the magnetic length, as l B = (cid:112) (cid:126) / ( | e | B ), wecan write [ ˆΠ α , ˆΠ β ] = − i (cid:126) sgn( B ) (cid:15) αβz /l B . This implies acanonical commutation relation between the dynamicalmomenta and inspires to define creation and annihilation operators in the usual way asˆ a † = l b √ (cid:126) ( ˆΠ x + i sgn( B ) ˆΠ y ) , ˆ a = l b √ (cid:126) ( ˆΠ x − i sgn( B ) ˆΠ y ) , (A2)which obey [ˆ a, ˆ a † ] = 1. Note that the definition of theseoperators depends on the sign of the B -field, which is adirect consequence of Eq. (A1). We note in passing thatall of the above did not require specifying a gauge for A α .The Landau level spectrum of a Dirac Hamiltonian ofthe form H = (cid:126) v F (Γ x q x + Γ y q y )is then straightforwardly obtained by making the substi-tution (cid:126) q µ → ˆΠ µ . Squaring the Hamiltonian yields H = v F ( ˆΠ x + ˆΠ y ) + v F [ ˆΠ x , ˆΠ y ]Γ x Γ y = v F (cid:126) l b (2ˆ a † ˆ a − B ) τ z )We use that ˆ a † ˆ a = n for standard oscillator wave func-tions ϕ n , i.e. ˆ a † ˆ aϕ n = nϕ n , ˆ aϕ n = √ nϕ n − andˆ a † ϕ n = √ n + 1 ϕ n +1 . One therefore obtains the Landaulevel energies E ± ( n ) = ± (cid:112) ξ n, n = 1 , , . . . ,ξ ≡ v F (cid:126) l b (A3)Each of the E ± ( n ) is two-fold degenerate because of thevalley degree of freedom, in addition to an N φ = A/ πl b degeneracy where A is the area of the system.We find the corresponding eigenstates by taking acloser look at the explicit expression for the Hamiltonianin each valley. Writing H ν for the Hamiltonian in valley ν = ± we have H ν = ν (cid:32) √ ξ ˆ a √ ξ ˆ a † (cid:33) . The eigenstates belonging to the eigenvalues E ± ( n ) ofEq. (A) are easily obtained as | Ψ nν ± (cid:105) = 1 √ (cid:32) | ϕ n − ,k (cid:105)± ν | ϕ n,k (cid:105) (cid:33) , (A4)In addition to the states | Ψ nν ± (cid:105) ( n = 1 , , . . . ) thereare also zero mode states | Ψ ν (cid:105) , one for each valley, whichhave zero energy ( E = 0). Inspecting of Eq. (18) revealsthat these states are given in each of the valleys as | Ψ ν (cid:105) = (cid:32) | ϕ ,k (cid:105) (cid:33) . (A5)5We stress that this implies the zero modes are localizedon opposite sublattices for the two valleys, as we hadexchanged sublattices for the (cid:126)K − valley.We proceed to consider the effect of symmetry breakingterms on the Landau level spectrum. Specifically, weconsider first the set of time-reversal symmetry invariantmasses (cid:126)m which enter the Hamiltonian as H (cid:126)m = (cid:126)m · (cid:126) Γ = m ν x + m ν y + m ν z τ z For small masses we may use perturbation theory tostudy the splitting or shifting of Landau levels. It turnsout however that the exact energies can be obtained inthe presence of H (cid:126)m . The energies are found by squaringthe Hamiltonian H = H + H (cid:126)m , which gives H = v F ( ˆΠ x + ˆΠ y ) + v F [ ˆΠ x , ˆΠ y ]Γ x Γ y + m = v F (cid:126) l b (2ˆ a † ˆ a − B ) τ z ) + m , where we use the anticommutation relations of the Γ-matrices. We directly find the energies E ± ( n ) = ± (cid:112) ξ n + m , n = 1 , , . . . (A6)Expanding the square root (cid:112) ξ n (cid:112) m / ξ n insmall m / ξ n yields the same result as second orderperturbation theory.The Landau level spectrum in the presence of massescan alternatively obtained by using the relation [Ω i , Γ j ] =2 i(cid:15) ijk Γ k to construct a unitary matrix U which rotatesthe vector (cid:126)m so that one has U † (cid:126)m · (cid:126) Γ U = m Γ , with m = | (cid:126)m | . Such transformation block diagonalizes theHamiltonian and we obtain the energies in a more directway. In particular, we can employ the unitary rotation tofind the energies of the zero modes for the case of massiveDirac fermions. In the rotated basis we can constructzero mode states in the same way as before, which willhave energies E ν = − νm (A7)where ν = ± represents the valley degree of freedom.We conclude by taking into account the time-reversalsymmetry breaking but chiral symmetry preserving mass η entering as H η = ητ z . Since such this term is a scalarunder chiral rotations generated by Ω i we may consider H = H + H (cid:126)m + H η and use the chiral rotation U toblock diagonalize the Hamiltonian. The total mass termthen is mν z τ z + ητ z , which directly leads to the energies E ν ± ( n ) = ± (cid:112) ξ n + ( νm + η ) (A8)for the Landau levels n = 1 , , . . . . The presense of bothof these masses leads to a splitting of Landau levels,whereas the presence of either only shifts the energies.The presence of a time-reversal breaking mass breaksparticle-hole symmetry in the n = 0 Landau level. Specif-ically, the zero mode energies become E ν = − νm − η. (A9) Appendix B: Setup of Hartree-Fock calculations
Our Hartree-Fock calculations in the presence of thestrain superlattice and with interacting Hamiltonian (21)follow the scheme of Ref. 44. In particular, we choose theunit cell of the unstrained lattice such that it containssix honeycomb lattice sites, as show in Fig. 6 of the maintext. This allows intra-sublattice mean-field structuresto form, corresponding to modulation vectors K ± , whichconnect the Dirac points of the honeycomb lattice. Interms of the elementary graphene lattice vectors (cid:126)a = a (1 , √ / , (cid:126)a = a (1 , −√ / , the lattice vectors of the six-site unit cell are given by (cid:126)b = 2 (cid:126)a + (cid:126)a , (cid:126)b = (cid:126)a − (cid:126)a . They are shown in Fig. 6, together with the folded BZ.Figure 6 also shows how, in the presence of the periodicstrain superlattice, the superlattice unit cell is defined.The superlattice vectors are (cid:126)b and λ(cid:126)b , in terms of thesuperlattive wave length λ .The electronic Hamiltonian is expressed in terms ofthe fermion annihilation (and corresponding creation)operators ˆ ψ αi ( l, (cid:126)x ). Here α labels th sublattice ( A/B ), i = 1 , , l = 1 , . . . , λ lables the (six-site) unit cells in the superlattice unit cell, and (cid:126)x is aposition index for the superlattice unit cell. The Fouriertransforms is defined asˆ ψ αi ( l, (cid:126)k ) = 1 √ N (cid:88) (cid:126)x ˆ ψ αi ( l, (cid:126)x ) e − i(cid:126)x · (cid:126)k . The interacting Hamiltonian of Eq. (21) consists of twoterms: the NN interaction and the NNN interaction. Inmomentum space the NN interaction ˆ H V is given byˆ H V = V N (cid:88) (cid:126)k(cid:126)k (cid:48) (cid:126)q ˆ ψ † Ai ( l, (cid:126)k ) ˆ ψ Ai ( l, (cid:126)k − (cid:126)q ) X il,jl (cid:48) ( (cid:126)q ) × ˆ ψ † Bj ( l (cid:48) , (cid:126)k (cid:48) ) ˆ ψ Bj ( l (cid:48) , (cid:126)k (cid:48) + (cid:126)q ) , (B1)where repeated indices are summed. The matrix function X ( (cid:126)q ) connects NNs, and it is convenient to decompose itin the following way X il,jl (cid:48) ( (cid:126)q ) = X − ij ( (cid:126)q ) δ l,l (cid:48) − + X ij ( (cid:126)q ) δ l,l (cid:48) + X + ij ( (cid:126)q ) δ l,l (cid:48) +1 . Clearly, X ij ( (cid:126)q ) connects sites within the same (six-site)unit cell, and it is explicitly given by X ij ( (cid:126)q ) = e i(cid:126)b · (cid:126)q . (B2)The functions X ± ij ( (cid:126)q ) connect sites in different (six-site)unit cells and each only have a single nonzero entry. Theyare given by X +12 ( (cid:126)q ) = e − i ( (cid:126)b + (cid:126)b ) · (cid:126)q and X − ( (cid:126)q ) = e i(cid:126)b · (cid:126)q .6Similarly, the NNN interaction Hamiltonian is givenby ˆ H V = V N (cid:88) (cid:126)k(cid:126)k (cid:48) (cid:126)q ˆ ψ † αi ( l, (cid:126)k ) ˆ ψ αi ( l, (cid:126)k − (cid:126)q ) Y αil,jl (cid:48) ( (cid:126)q ) × ˆ ψ † αj ( l (cid:48) , (cid:126)k (cid:48) ) ˆ ψ αj ( l (cid:48) , (cid:126)k (cid:48) + (cid:126)q ) (B3)The functions Y α ( (cid:126)q ) connect NNNs on each sublattice α . They are directly obtained from Ref. 44, taking intoaccount the additional superlattice index l (in the sameway as for X ( (cid:126)q )).The quartic interactions of the interacting Hamiltoni-ans, schematically written as ˆ ψ † i ˆ ψ i ˆ ψ † j ˆ ψ j , are decoupled inthe standard mean-field way as (written schematically) → ˆ ψ † i ˆ ψ i (cid:104) ˆ ψ † j ˆ ψ j (cid:105) + (cid:104) ˆ ψ † i ˆ ψ i (cid:105) ˆ ψ † j ˆ ψ j − (cid:104) ˆ ψ † i ˆ ψ i (cid:105)(cid:104) ˆ ψ † j ˆ ψ j (cid:105) , → − ˆ ψ † i ˆ ψ j (cid:104) ˆ ψ † j ˆ ψ i (cid:105) − (cid:104) ˆ ψ † i ˆ ψ j (cid:105) ˆ ψ † j ˆ ψ i + (cid:104) ˆ ψ † i ˆ ψ j (cid:105)(cid:104) ˆ ψ † j ˆ ψ j (cid:105) , the first line representing charge density order and thesecond bond density order.In terms of the actual superlattice electrons, the charge density order parameter is defined as ρ αi ( l ) = 1 N (cid:88) (cid:126)k (cid:104) ˆ ψ † αi ( l, (cid:126)k ) ˆ ψ αi ( l, (cid:126)k ) (cid:105) . (B4)Bond order parameters are defined as straightforwardgeneralizations of Ref. 44. Of particular interest in ourcase is the QH order parameter, constructed from NNNbond order and defined by Eqs. (27) and (28). To obtainEq. (27) we decompose Y αil,jl (cid:48) ( (cid:126)k − (cid:126)k (cid:48) ), wich arises due tothe reordering of (B3), as Y αil,jl (cid:48) ( (cid:126)k − (cid:126)k (cid:48) ) = (cid:88) µ =1 Q αµil,jl (cid:48) ( (cid:126)k ) Q αµil,jl (cid:48) ( − (cid:126)k (cid:48) ) . (B5)The sum over µ follows from the three ways in whichNNN sites may be connected. Note that Q αµil,jl (cid:48) ( − (cid:126)k ) = Q αµ ∗ il,jl (cid:48) ( (cid:126)k ). We can now NNN bond order parameters usedto construct the QH order parameter of Eq. (28). TheNNN bond order mean-fields are defined by χ αil,jl (cid:48) = 1 N (cid:88) (cid:126)k Q α ∗ il,jl (cid:48) ( (cid:126)k ) (cid:104) ˆ ψ † αi ( l, (cid:126)k ) ˆ ψ αj ( l (cid:48) , (cid:126)k ) (cid:105) . (B6)The QH order paramater ∆ QH ( l ) is defined so that eachunit cell labeled by l is associated with 2 × [1] M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. , 3045(2010).[2] X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. , 1057(2011).[3] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang,Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A.Firsov, Science , 666 (2004).[4] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S.Novoselov, and A. K. Geim, Rev. Mod. Phys. , 109(2009).[5] L. Fu, C. L. Kane, and E. J. Mele, Phys. Rev. Lett. ,106803 (2007).[6] G. W. Semenoff, Phys. Rev. Lett. , 2449 (1984).[7] F. D. M. Haldane, Phys. Rev. Lett. , 2015 (1988).[8] T. H. Hsieh, H. Lin, J. Liu, W. Duan, A. Bansil, and L.Fu, Nat. Commun. , 982 (2012).[9] Y. Tanaka, Z. Ren, T. Sato, K. Nakayama, S. Souma, T.Takahashi, K. Segawa, and Y. Ando, Nat. Phys. , 800(2012).[10] P. Dziawa, B.J. Kowalski, K. Dybko, R. Buczko, A.Szczerbakow, M. Szot, E. (cid:32)Lusakowska, T. Balasubrama-nian, B.M. Wojek, M.H. Berntsen, O. Tjernberg, and T.Story, Nat. Matter. , 1023 (2012).[11] S.-Y. Xu, C. Liu, N. Alidoust, M. Neupane, D. Qian,I. Belopolski, J.D. Denlinger, Y.J. Wang, H. Lin, L.A.Wray, G. Landolt, B. Slomski, J.H. Dil, A. Marcinkova,E. Morosan, Q. Gibson, R. Sankar, F.C. Chou, R.J.Cava, A. Bansil, and M.Z. Hasan, Nat. Commun. , 1192(2012). [12] Y. Ando and L. Fu, Ann. Rev. Cond. Matt., , 361(2015).[13] Y. Okada, M. Serbyn, H. Lin, D. Walkup, W. Zhou,C. Dhital, M. Neupane, S. Xu, Y. Wang, R. Sankar, F.Chou, A. Bansil, M. Z Hasan, S.D. Wilson, L. Fu, andV. Madhavan, Science, , 1496 (2013).[14] I. Zeljkovic, Y. Okada, M. Serbyn, R. Sankar, D. Walkup,W. Zhou, J. Liu, G. Chang, Y. J. Wang, M. Z. Hasan,F. Chou, H. Lin, A. Bansil, L. Fu, V. Madhavan, Nat.Mater. , 318 (2015).[15] M. Serbyn and L. Fu, Phys. Rev. B , 035402 (2014).[16] L. Fu, Phys. Rev. Lett. , 106802 (2011).[17] F. Guinea, M. I. Katsnelson, and A. K. Geim, Nat. Phys. , 30 (2010).[18] M. A. H. Vozmediano, M. I. Katsnelson, and F. Guinea,Phys. Reps.-Rev. Sec. Phys. Lett. , 109 (2010).[19] B. Amorim, A. Cortijo, F. de Juan, A. G. Grushin, F.Guinea, A. Guti´errez-Rubio, H. Ochoa, V. Parente, R.Roldan, P. San-Jose, J. Schiefele, M. Sturla, and M. A. H.Vozmediano, arXiv:1503.00747v2 (2015).[20] K. K. Gomes, W. Mar, W. Ko, F. Guinea, and H. C.Manoharan, Nature , 306 (2012).[21] N. Levy, S. A. Burke, K. L. Meaker, M. Panlasigui, A.Zettl, F. Guinea, A. H. Castro Neto, and M. F. Crommie,Science , 544 (2010).[22] V. P. Gusynin, V. A. Miransky, and I. A. Shovkovy, Phys.Rev. Lett. , 3499 (1994).[23] I. F. Herbut, Phys. Rev. B , 205433 (2008).[24] P. Ghaemi, J. Cayssol, D. N. Sheng, and A. Vishwanath, Phys. Rev. Lett. , 266801 (2012).[25] D. Abanin and D. A. Pesin, Phys. Rev. Lett. , 066802(2012).[26] B. Roy and I. F. Herbut, Phys. Rev. B , 045425 (2013).[27] B. Roy and J. D. Sau, Phys. Rev. B , 075427 (2014).[28] B. Roy, F. F. Assaad, and I. F. Herbut, Phys. Rev. X ,021042 (2014).[29] M. A. H. Vozmediano, F. Guinea, and M. I. Katsnelson,Phys. Rev. B , 075422 (2008).[30] E. Tang and L. Fu, Nat. Phys. , 964 (2014).[31] A. Y. Sipatov, Funct. Mater. , 374 (2009).[32] L. S. Palatnik and A. I. Fedorenko, J. Cryst. Growth ,917 (1981).[33] G. Springholz and K. Wiesauer, Phys. Rev. Lett. ,015507 (2001).[34] S. Gopalakrishnan, P. Ghaemi, and S. Ryu, Phys. Rev.B , 081403 (2012).[35] F. de Juan, Phys. Rev. B , 125419 (2013).[36] J. Liu, W. Duan, and L. Fu, Phys. Rev. B , 241303(2013).[37] C. W. J. Beenakker, Rev. Mod. Phys. , 1337 (2008).[38] C.-Y. Hou, C. D. Chamon, and C. Mudry, Phys. Rev.Lett. , 186809 (2007).[39] C. Fang and L. Fu, Phys. Rev. B , 161105 (2015). [40] J. L. Manes, Phys. Rev. B , 045430 (2007).[41] J. R. Wallbank, M. Mucha-Kruczynski, and V. I. Falko,Phys. Rev. B , 155415 (2013).[42] P. Ghaemi, S. Gopalakrishnan, and S. Ryu, Phys. Rev.B , 155422 (2013).[43] S. Raghu, X.-L. Qi, C. Honerkamp, and S.-C. Zhang,Phys. Rev. Lett. , 156401 (2008).[44] A. G. Grushin, E. V. Castro, A. Cortijo, F. de Juan,M. A. H. Vozmediano, and B. Valenzuela, Phys. Rev. B , 085136 (2013).[45] C. Weeks and M. Franz, Physical Review B 81 (2010).[46] N. B. Kopnin, T. T. Heikkil¨a, and G. E. Volovik, Phys.Rev. B , 220503 (2011).[47] B. Uchoa and Y. Barlas, Phys. Rev. Lett. , 046604(2013).[48] M. Daghofer and M. Hohenadler, Phys. Rev. B ,035103 (2014).[49] N. A. Garcia-Martinez, A. G. Grushin, T. Neupert, B.Valenzuela, and E. V. Castro, Phys. Rev. B , 245123(2013).[50] T. Duric, N. Chancellor, and I. F. Herbut, Phys. Rev. B89