Open Momentum Space Method for Hofstadter Butterfly and the Quantized Lorentz Susceptibility
OOpen Momentum Space Method for Hofstadter Butterfly and the Quantized LorentzSusceptibility
Biao Lian, Fang Xie, and B. Andrei Bernevig Department of Physics, Princeton University, Princeton, New Jersey 08544, USA (Dated: February 10, 2021)We develop a generic k · p open momentum space method for calculating the Hofstadter butterflyof both continuum (Moir´e) models and tight-binding models, where the quasimomentum is directlysubstituted by the Landau level (LL) operators. By taking a LL cutoff (and a reciprocal latticecutoff for continuum models), one obtains the Hofstadter butterfly with in-gap spectral flows. Forcontinuum models such as the Moir´e model for twisted bilayer graphene, our method gives a sparseHamiltonian, making it much more efficient than existing methods. The spectral flows in theHofstadter gaps can be understood as edge states on a momentum space boundary, from whichone can determine the two integers ( t ν , s ν ) of a gap ν satisfying the Diophantine equation. Thespectral flows can also be removed to obtain a clear Hofstadter butterfly. While t ν is known as theChern number, our theory identifies s ν as a dual Chern number for the momentum space, whichcorresponds to a quantized Lorentz susceptibility γ xy = eBs ν . Two-dimensional (2D) lattice electrons in large mag-netic fields are known to exhibit Hofstadter butterflyspectra [1]. Conventionally, the Hofstadter butterfly iscalculated at rational fluxes per unit cell ϕ = 2 πp/q in abasis with translation symmetry of q unit cells, where p and q are coprime integers. The calculation ofteninvolves a complicated construction of the matrix ele-ments. In particular, for continuum k · p models obtainedfrom plane wave expansions such as the Moir´e modelfor twisted bilayer graphene (TBG) [2–4], the HofstadterHamiltonian matrix is infinite dimensional and dense [5–10], which requires a large cutoff for the spectrum toconverge.In contrast, the Landau levels (LLs) of a k · p Hamil-tonian at small magnetic fields can be calculated by sim-ply substituting the quasimomentum k = ( k x , k y ) with( a + a † √ (cid:96) , a − a † i √ (cid:96) ), where a and a † are the LL lowering andraising operators, and (cid:96) is the magnetic length [11]. Inthis letter, we demonstrate that such a substitution witha LL cutoff (and a reciprocal lattice cutoff for contin-uum models) provides an efficient method for calculatingthe Hofstadter butterfly in large magnetic fields, whichgreatly simplifies the Hamiltonian matrix elements [12].In particular, for continuum models, this method yieldsa sparse Hamiltonian, whose spectrum can be efficientlycalculated by the shift-and-invert Lanczos method.The method can be understood as an open momentumspace calculation, where the smaller of the momentum-space LL wavefunction radius cutoff and reciprocal latticeradius cutoff plays the role of a momentum space bound-ary. As a result, the spectrum contains not only theHofstadter butterfly, but also in-gap spectral flow levels[12, 13] which can be understood as “momentum spaceedge states”. We show that the spectral flows of theseedges allow us to determine the two integers ( t ν , s ν ) ina Hofstadter gap ν satisfying the Diophantine equation[14–16], where t ν is the Chern number of the gap. More-over, we show that s ν can be interpreted as a dual Chernnumber for the momentum space, which yields a quan- tized Lorentz susceptibility (Eq. (15)). Furthermore,by identifying and removing the momentum space edgestates, one can obtain the Hofstadter butterfly withoutspectral flows. We demonstrate our method for both con-tinuum models and tight binding models in a 2D periodiclattice. We shall denote the lattice Bravais vectors as d and d , and the reciprocal vectors as g and g , whichsatisfy g i · d j = 2 πδ ij ( i, j = 1 , Continuum models.
At zero magnetic field, a contin-uum model can be written in the real space basis | r , α (cid:105) as [3–7, 17] H αβ ( r ) = (cid:15) αβ ( − i ∇ ) + (cid:88) j V αβj e i q αβj · r , (1)where r = ( x, y ) is the real space position, − i ∇ = − i ( ∂ x , ∂ y ) is the canonical momentum, and we assumethere are M intrinsic orbitals labeled by α, β . (cid:15) αβ ( − i ∇ )and V αβj e i q αβj · r are the electron kinetic term in free spaceand the periodic lattice potential between orbitals β and α , respectively. If one denote Q ∈ g Z + g Z as the recip-rocal lattice, and choose the momentum origin of orbital α at p α , one can define a momentum lattice Q α = p α + Q for orbital α , and q αβj in Eq. (1) must be the difference Q (cid:48) α − Q β between some sites Q (cid:48) α and Q β (see supplemen-tary material (SM) [18] Sec. S2A). Generically, one canalways fix all p α = ; however, in certain models (e.g.,the TBG model [2]) nonzero p α choices are preferred.One can transform the zero-magnetic-field Hamilto-nian (1) into the momentum eigenbasis | k , Q α , α (cid:105) = (cid:82) d r e i ( k + Q α ) · r | r , α (cid:105) , where k is in the first Brillouinzone (BZ). The momentum space Hamiltonian under ba-sis | k , Q α , α (cid:105) then takes the form [2, 18] H αβ Q (cid:48) α Q β ( k ) = (cid:15) αβ ( k + Q β ) δ Q (cid:48) α Q β + (cid:88) j V αβj δ Q (cid:48) α , Q β + q αβj . (2)When a uniform out-of-plane magnetic field B = B ˆ z is added, − i ∇ in Eq. (1) is replaced by the kinematic a r X i v : . [ c ond - m a t . m e s - h a ll ] F e b momentum Π = − i ∇ − A ( r ), where A ( r ) is the vec-tor potential satisfying ∂ x A y − ∂ y A x = B . The kineticmomentum satisfies [Π x , Π y ] = i/(cid:96) , where (cid:96) = 1 / √ B isthe magnetic length. We also define the guiding center R = r − (cid:96) ˆ z × Π , which satisfies [ R x , R y ] = − i(cid:96) , and[ R , Π ] = 0.The usual Hofstadter method for continuum modelsemploys the Landau basis defined by eigenstates of R x and Π , which has complicated matrix elements [5–10].Here we shall take a different basis, under which we provethe nonzero magnetic field Hamiltonian can be simplyobtained by the zero-field momentum-space Hamiltonian(2) with the substitution of Eq. (5).We define R ˆ τ = R · ˆ τ as the guiding center alongunit vector ˆ τ , where we choose ˆ τ · (ˆ z × g )ˆ τ · (ˆ z × g ) irrational. Wealso define a set of (linearly dependent) LL operators a Q α = (cid:96) √ [Π x − Q α,x − k ,x + i (Π y − Q α,y − k ,y )] andtheir conjugates a † Q α associated with momentum sites Q α , where k = ( k ,x , k ,y ) is a freely chosen real vectorwhich we call the center momentum . We then constructan orthonormal basis | λ, Q α , n, α (cid:105) for orbital α and re-ciprocal site Q α by requiring R ˆ τ | λ, Q α , n, α (cid:105) = [ λ − (cid:96) ˆ τ · (ˆ z × Q α )] | λ, Q α , n, α (cid:105) ,a † Q α a Q α | λ, Q α , n, α (cid:105) = n | λ, Q α , n, α (cid:105) . (3)Here n ≥ λ is a realnumber chosen in the set λ + (cid:96) ˆ τ · (ˆ z × g ) Z + (cid:96) ˆ τ · (ˆ z × g ) Z representing the set, or abstractly, λ ∈ R / [ (cid:96) ˆ τ · (ˆ z × g ) Z + (cid:96) ˆ τ · (ˆ z × g ) Z ] (see SM [18] Sec. S2B).It can then be proved that all the states | λ, Q α , n, α (cid:105) form a complete basis for the continuum model satisfying (cid:104) λ, Q α , n, α | λ (cid:48) , Q (cid:48) β , n (cid:48) , β (cid:105) = δ λλ (cid:48) δ Q α , Q (cid:48) β δ nn (cid:48) δ αβ .The above basis | λ, Q α , n, α (cid:105) is advantageous becausethe nonzero-magnetic-field Hamiltonian is diagonal in λ and independent of λ . In SM [18] Sec. S2B, we show theHamiltonian in a fixed λ subspace is H λ,αβ Q (cid:48) α Q β = (cid:15) αβ (ˆ κ Q β + k + Q β ) δ Q (cid:48) α Q β + (cid:88) j V j δ Q (cid:48) α , Q β + q αβj , (4)where we have defined ˆ κ Q α = √ (cid:96) ( a Q α + a † Q α , − ia Q α + ia † Q α ). Without ambiguity, we can drop the subindex Q α and simplify ( a Q α , a † Q α ) as ( a, a † ), which acts as a | λ, Q α , n, α (cid:105) = √ n | λ, Q α , n − , α (cid:105) and a † | λ, Q α , n, α (cid:105) = √ n + 1 | λ, Q α , n + 1 , α (cid:105) . The Hamiltonian (4) is then ex-actly the zero-field Hamiltonian H αβ Q (cid:48) α Q β ( k ) in Eq. (2)with the substitution k x → a + a † √ (cid:96) + k ,x , k y → a − a † i √ (cid:96) + k ,y , (5)as we claimed earlier. One then only need calculate thespectrum for a fixed λ . Different λ and λ (cid:48) subspaceshave identical spectra, but have eigenstates differing bydisplacement λ − λ (cid:48) in the ˆ τ direction ( R ˆ τ eigenvalue). q q q K M K M ’ Γ M Q k k Q | |= L / l κ (a) (b) A k = π N L / l A k = N Q Ω BZ FIG. 1. (a) When ϕ/ π < N Q /N L , the momentum space(the shaded area) has a circular boundary of radius √ N L /(cid:96) .(b) When ϕ/ π > N Q /N L , the momentum space boundary isthe reciprocal lattice boundary enclosing N Q BZs (the shadedarea).
To numerically calculate the spectrum of Hamiltonian(4), one can fix a center momentum k , take a LL cutoff n ≤ N L , and take a cutoff of reciprocal lattice Q α at aboundary enclosing N Q BZs. This yields a Hamiltonianof size
M N L N Q for M intrinsic orbitals. If (cid:15) ( k ) only con-tains polynomials up to ∆-th power of k , and the numberof q αβj is finite, (cid:104) λ, Q (cid:48) α , m, α | H | λ, Q β , n, β (cid:105) will be zerofor | m − n | > ∆ or | Q (cid:48) α − Q β | > max( | q αβj | ), so the Hamil-tonian H is a sparse matrix. The low-energy eigenstatesand spectrum can then be efficiently calculated by theLanczos algorithm.The cutoffs N Q and N L , however, lead to spectral flowsin the Hofstadter gaps due to the absence of periodicboundary conditions [12, 13]. As an example, we cal-culate the Hofstadter butterfly of the TBG continuummodel defined on a honeycomb momentum lattice Q α [2], which has a Dirac kinetic term (cid:15) ( k ) = v F σ ∗ · k , and2 × V j between the nearest momentumsites, where σ ∗ = ( σ x , − σ y ) are the Pauli matrices (SM[18] Sec. S3). Fig. 2(a) shows the TBG spectrum at twistangle θ = 2 . ◦ versus the flux per unit cell ϕ = B | d × d | ,where we take N Q = 37 and N L = 60. Besides the Hof-stadter butterfly, one can see numerous in-gap spectralflow levels.The in-gap spectral flows are generically due to thepresence of boundaries which host edge states [12, 13].Here, as illustrated in Fig. 1(a) and (b), the cutoff N Q sets a boundary of momentum radius (cid:113) N Q Ω BZ π enclosing N Q BZs, where Ω BZ = 4 π / | d × d | is the BZ area; whilethe LL cutoff N L yields a boundary (cid:113) (cid:104) ˆ κ Q α (cid:105) ≤ √ N L (cid:96) for ˆ κ Q α in the Hamiltonian (4). The smaller value of (cid:113) N Q Ω BZ π and √ N L (cid:96) then serves as a momentum spaceboundary radius for Hamiltonian (4) (SM [18] Sec. S2C),which gives rise to edge states.The momentum space edge state levels then generatethe spectral flows versus the magnetic field B . This canbe understood from the Diophantine equation [14–16, 18]satisfied by the ν -th Hofstadter gap ( ν ∈ Z ) at ϕ = 2 πp/q flux per unit cell: t ν p + s ν q = ν , (6)where ( t ν , s ν ) are two integer quantum numbers charac-terizing the gap. In particular, t ν is the Chern numberof the gap. It is often rewritten as t ν ϕ π + s ν = ρ , (7)where ρ = ν/q is the number of occupied bulk states perunit cell in the gap [18–20]. Here it is more useful torewrite it in a dual form t ν + s ν πϕ = ρ K , (8)where ρ K = 2 πρ/ϕ = ν/p . In SM [18] Sec. S2F, weshow that ρ K gives the number of occupied bulk states per BZ in the gap for the Hamiltonian (4) at a fixed λ .Furthermore, in Eq. (16), we show that s ν plays the roleof a dual Chern number for the momentum space. Eq.(8) then determines the in-gap spectral flows (Fig. 2(a))in two different regimes as follows.In the first regime ϕ/ π < N Q /N L , the momentumspace boundary is a circle enclosing a ϕ -dependent area A K = 2 πN L /(cid:96) centered at k (Fig. 1(a)). In a gap,the total number of occupied states in the momentumarea A K is N K = ρ K A K / Ω BZ = N L ρ K ϕ/ π = N L ρ .Therefore, by Eq. (8) we have N K = N L ( t ν ϕ/ π + s ν ) . (9)In a bulk gap, the number of occupied states N K canonly change by pumping edge states into (out of) thebulk. Therefore, the edge states necessarily produce in-gap spectral flows, where the rate of flowing levels isd N K / d( ϕ/ π ) = N L t ν by Eq. (9). Fig. 2(c) showsa gap in this regime, where the midgap line (dashedline) crosses 16 levels as ϕ/ π increases from 0 .
25 to 0 . N L = 60. The flow rate is then d N K / d( ϕ/ π ) =16 / (0 . − .
25) = 64, so we can identify the Chern numberof the gap as the integer closest to N − L d N K / d( ϕ/ π ) =1 .
07, namely, t ν = 1. Further, at ϕ/ π = 0 .
5, we countedthere are N K = 32 levels from zero energy to the midgapenergy (for TBG, N K = 0 is at zero energy (SM [18] Sec.S2G)), so we find the gap has s ν = N K /N L − t ν ϕ/ π = 0from Eq. (9).In the second regime ϕ/ π > N Q /N L , the momentumspace boundary is given by cutoff N Q , which encloses a ϕ independent area A K = N Q Ω BZ (Fig. 1(b)). The num-ber of occupied states N K = ρ K A K / Ω BZ in momentumarea A K in a gap is then N K = N Q ( t ν + 2 πs ν /ϕ ) . (10)This yields a spectral flow rate d N K / d(2 π/ϕ ) = N Q s ν .Besides, for TBG which has a Dirac kinetic term, thereare 2 N Q horizontal levels at ϕ/ π > N Q /N L in Fig. E ( e V ) π / ϕ spuriousmodes E ( e V ) ϕ /2 π E ( e V ) ϕ /2 π E ( e V ) ϕ /2 π (1,0) ( , ) (0,1)(0,1) (1,1) ( , - ) (-1,0)(-1,-1) (- , ) (0,-1)(0,-1) ( , - ) ( , - ) (- , ) (- , ) (- , ) (a) (b)(c) (d) FIG. 2. (a) Hofstadter butterfly and spectral flow of θ = 2 . ◦ TBG with N Q = 37 and N L = 60, and k at Γ point. Thehorizontal axis ϕ/ π is linearly plotted in [0 , − π/ϕ in [1 , ∞ ]. (b) The Hofstadter butterfly afterdeleting the edge states (with w = min { (cid:96) − , . √ Ω BZ } , P c =0 . t ν , s ν ) in the gaps. (c) Zoom-in plot in the regime ϕ/ π < N Q /N L . (d) Zoom-in plot in the regime ϕ/ π >N Q /N L . N L (see SM [18] Sec. S3). These spurious levels shouldbe excluded when counting N K . Fig. 2(d) shows agap in this regime, where the midgap line crosses 6 lev-els (excluding the spurious modes) as 2 π/ϕ decreasesfrom 2 / /
2, and N Q = 37. The flow rate is thend N K / d(2 π/ϕ ) = 6 / (2 / − /
2) = 36, thus s ν can be iden-tified as the integer closest to N − Q d N K / d(2 π/ϕ ) = 0 . s ν = 1. Further, we counted there are N K = 55levels (excluding the spurious modes) between midgapand zero energy at 2 π/ϕ = 1 /
2, thus the gap has a Chernnumber t ν = N K /N Q − πs ν /ϕ = 1 from Eq. (10).We note that models with a Dirac kinetic term (cid:15) ( k ) = v F σ ∗ · k would have N K = 0 defined at half filling (zeroenergy for TBG), while models with a lower-bounded ki-netic term (e.g., (cid:15) ( k ) = k / m ) would have N K = 0below the lowest band (SM [18] Sec. S2G). More generi-cally, if a gap persists below and above ϕ/ π = N Q /N L ,one can identify t ν and s ν separately from the spectralflow rates at small and large ϕ , after which one can obtain N K of the gap from Eq. (9) or (10).The edge states and spurious modes can be easily re-moved from the spectrum. We define a boundary projec-tor P κ b ,w onto basis | λ, Q α , n, α (cid:105) with n > ( κ b − w ) (cid:96) / w >
0, where κ b ≈ min { √ N L (cid:96) , (cid:113) N Q Ω BZ π } isthe radius of momentum space boundary. We can thenidentify the eigenstates with (cid:104) P κ b ,w (cid:105) > P c above certainvalue P c ∈ [0 ,
1] as momentum space edge states withindistance w to the boundary, and delete them to obtain abulk Hofstadter spectrum. For example, Fig. 2(b) is ob-tained by setting w = min { (cid:96) − , . √ Ω BZ } and P c = 0 . Tight-binding models.
The substitution (5) can also beemployed to calculate the Hofstadter butterfly of tight-binding models. Given the position u α of each Wannierorbital α in a unit cell in the continuum space, we denoteorbital α at position D + u α as | D , α (cid:105) , where D ∈ d Z + d Z is the lattice vector. The Hamiltonian under Peierlssubstitution [21–23] then takes the form H = (cid:88) j,α,β t αβj T D j + u α − u β , (11)where D j ∈ d Z + d Z , t αβj is the hopping from | D , β (cid:105) to | D + D j , α (cid:105) , and T D j + u α − u β = (cid:88) D e i (cid:82) cαβ A ( r ) · d r | D + D j , α (cid:105)(cid:104) D , β | (12)is the translation operator, with c αβ being the straightline segment from D + u β to D + D j + u α . At zero mag-netic field, the Hamiltonian can be transformed into themomentum space basis | k , α (cid:105) = (cid:80) D e i k · ( D + u α ) | D , α (cid:105) as H αβ ( k ) = (cid:88) j t αβj e − i k · ( D j + u α − u β ) . (13)At nonzero magnetic field, we define a basis as | λ, n, α (cid:105) = (cid:80) D | D , α (cid:105)(cid:104) D + u α , α | λ, , n, α (cid:105) , where | D + u α , α (cid:105) is the continuum space position eigenstate at po-sition r = D + u α , | λ, , n, α (cid:105) is the state defined in Eq.(3) in the continuum space at reciprocal site , and λ ∈ R / [ (cid:96) ˆ τ · (ˆ z × g ) Z + (cid:96) ˆ τ · (ˆ z × g ) Z ]. One can then show that | λ, n, α (cid:105) forms a complete orthonormal basis of Hamilto-nian (11) satisfying (cid:104) λ (cid:48) , n (cid:48) , β | λ, n, α (cid:105) = δ λλ (cid:48) δ n (cid:48) n δ βα (SM[18] Sec. S4A). Furthermore, T D j + u α − u β is diagonal in λ and takes the λ independent form T λ,αβ D j + u α − u β = e − i (ˆ κ + k ) · ( D j + u α − u β ) (14)in a fixed λ subspace between basis | λ, n, β (cid:105) and | λ, n (cid:48) , α (cid:105) ,where ˆ κ = √ (cid:96) ( a + a † , − ia + ia † ), with a | λ, n, α (cid:105) = √ n | λ, n − , α (cid:105) and a † | λ, n, α (cid:105) = √ n + 1 | λ, n + 1 , α (cid:105) (SM[18] Sec. S4A). Therefore, the nonzero magnetic fieldtight-binding Hamiltonian (11) in a fixed λ is given bythe zero-field momentum space Hamiltonian (13) withsubstitution (5). For nonstandard Peierls substitutionsalong nonstraight c αβ paths, e − i (ˆ κ + k ) · ( D j + u α − u β ) in Eq.(14) becomes the path-ordered integral P e − i (cid:82) cαβ (ˆ κ + k ) · d r (SM [18] Sec. S4B).The Hofstadter butterfly can then be numerically cal-culated with a LL cutoff, namely, n ≤ N L . Fig. 3(a) and(b) show the spectrum of the square lattice tight-bindingmodel H ( k ) = − cos k x − cos k y [1] with cutoffs N L = 100and N L = 500, respectively, where we set k = . The E ϕ /2 π bulk state n| ψ | edge state n| ψ | E ϕ /2 π E ϕ /2 π (a)(c)(d) (b)(e) FIG. 3. The Hofstadter spectrum for tight-binding model H ( k ) = − cos k x − cos k y with k = and LL cutoff (a) N L = 100 and (b) N L = 500. (c) Probability distribution ofa typical momentum space edge state in (b) versus LL number n . (d) Probability distribution of a typical bulk state in (b).(e) The Hofstadter butterfly obtained by deleting the edgestates in (b) (with w = 3 . (cid:96) − √ ϕ + 0 . P c = 0 . spectrum exhibits both the Hofstadter butterfly and thespectral flows, which can again be understood as momen-tum space edge states. Since tight-binding models haveno cutoff in the reciprocal lattice, the momentum spaceboundary is always at radius κ b = √ N L (cid:96) given by N L ,and the spectral flows always satisfy Eq. (9).One can define a boundary projector P κ b ,w onto basis | λ, n, α (cid:105) with n > ( κ b − w ) (cid:96) / w >
0, andidentify the eigenstates with (cid:104) P κ b ,w (cid:105) > P c for some P c ∈ [0 ,
1] as momentum space edge states. Fig. 3(c) and (d)show the probability of a typical edge state and bulk stateversus LL number n , respectively. By deleting the edgestates, one can obtain a high-quality Hofstadter butterflywithout spectral flows, as shown in Fig. 3(e) (where w =3 . (cid:96) − √ ϕ + 0 . P c = 0 . Quantized Lorentz susceptibility . The Chern number t ν is known to give a quantized Hall conductance via theKubo formula σ xy = i∂ ω (cid:82) d ω (cid:48) (cid:104) G ω + ω (cid:48) ˆ j x G ω (cid:48) ˆ j y (cid:105)| ω → = t ν e h [10], where G ω is the Green’s function at energy ω ,and ˆ j = (ˆ j x , ˆ j y ) is the uniform current operator. The du-ality between Eqs. (7) and (8) suggests that s ν behavesas a dual Chern number for the momentum space, thus s ν should also give a quantized response. Indeed, by not-ing that the natural momentum-space dual of the currentoperator ˆ j is the force operator ˆ F = e B × d R d t = ( ˆ F x , ˆ F y ),we find s ν leads to a quantized Lorentz susceptibility (SM[18] Sec. S6) γ xy = − i ∂∂ω (cid:90) d ω (cid:48) (cid:104) G ω + ω (cid:48) ˆ F x G ω (cid:48) ˆ F y (cid:105) (cid:12)(cid:12)(cid:12) ω → = eBs ν . (15)It yields a Lorentz force per unit cell F x = γ xy v y on thesystem when the lattice is moving at velocity v y . Fur-thermore, a formula similar to the Thouless-Kohmoto-Nightingale-den Nijs formula [10] at flux per unit cell ϕ = 2 πp/q can be derived for s ν (SM [18] Sec. S6B2): s ν = − i (cid:88) n ∈ occ (cid:90) d ∈ Ω M d d π ˆ z · (cid:104) ∂ d w n, d | × | ∂ d w n, d (cid:105) , (16)where Ω M is a torus with periods d and d /p serving asa “dual magnetic BZ”, | w n, d (cid:105) = e i(cid:96) − (ˆ z × R ) · d | ψ n,(cid:96) − ˆ z × d (cid:105) (see the explicit form in SM [18] Sec. S6C) is definedusing the Bloch eigenstates | ψ n, k (cid:105) of band n , and n runsover all occupied bands. Discussion . It is worth noting that the cutoffs in ourmethod affect the resolution but not the shape of theHofstadter butterfly. Our method greatly simplifies thematrix element construction compared to usual methods[1, 5], and require neither rational flux per unit cell norlarge magnetic unit cells. Moreover, it leads to a sparseHamiltonian for continuum models. At small magnetic fields, our method reduces to the LL calculations of k · p Hamiltonians expanded at center momentum k . Thelarge magnetic field spectrum is insensitive to the choiceof k . ACKNOWLEDGMENTS
Acknowledgments . We thank Michael Zaletel, HoiChun Po and Junyi Zhang for helpful discussions. B.L.acknowledge the support of Princeton Center for Theo-retical Science at Princeton University at the early stageof this work. B.A.B was supported by the DOE GrantNo. DE-SC0016239, the Schmidt Fund for Innovative Re-search, Simons Investigator Grant No. 404513, and thePackard Foundation. Further support was provided bythe NSF-EAGER No. DMR 1643312, NSF-MRSEC No.DMR-1420541 and DMR-2011750, ONR No. N00014-20-1-2303, Gordon and Betty Moore Foundation throughGrant GBMF8685 towards the Princeton theory pro-gram, BSF Israel US foundation No. 2018226, and thePrinceton Global Network Funds. [1] Douglas R. Hofstadter, “Energy levels and wave func-tions of bloch electrons in rational and irrational mag-netic fields,” Phys. Rev. B , 2239–2249 (1976).[2] Rafi Bistritzer and Allan H. MacDonald, “Moir´e bandsin twisted double-layer graphene,” Proceedings of theNational Academy of Sciences , 256802 (2007).[4] E. J. Mele, “Commensuration and interlayer coherencein twisted bilayer graphene,” Phys. Rev. B , 161405(2010).[5] R. Bistritzer and A. H. MacDonald, “Moir´e butterfliesin twisted bilayer graphene,” Phys. Rev. B , 035440(2011).[6] S. Janecek, M. Aichinger, and E. R. Hern´andez, “Two-dimensional bloch electrons in perpendicular magneticfields: An exact calculation of the hofstadter butterflyspectrum,” Phys. Rev. B , 235429 (2013).[7] Godfrey Gumbs, Andrii Iurov, Danhong Huang, andLiubov Zhemchuzhna, “Revealing hofstadter spectrumfor graphene in a periodic potential,” Phys. Rev. B ,241407 (2014).[8] Kasra Hejazi, Chunxiao Liu, and Leon Balents, “Lan-dau levels in twisted bilayer graphene and semiclassicalorbits,” Phys. Rev. B , 035115 (2019).[9] Ya-Hui Zhang, Hoi Chun Po, and T. Senthil, “LandauLevel Degeneracy in Twisted Bilayer Graphene: Role ofSymmetry Breaking,” arXiv e-prints , arXiv:1904.10452(2019), arXiv:1904.10452 [cond-mat.str-el].[10] D. J. Thouless, M. Kohmoto, M. P. Nightingale, andM. den Nijs, “Quantized hall conductance in a two-dimensional periodic potential,” Phys. Rev. Lett. ,405–408 (1982). [11] Roland Winkler, Spin-orbit coupling effects in two-dimensional electron and hole systems , Springer tractsin modern physics (Springer, Berlin, 2003).[12] Biao Lian, Fang Xie, and B. Andrei Bernevig, “TheLandau Level of Fragile Topology,” arXiv e-prints, arXiv:1811.11786 (2018), arXiv:1811.11786 [cond-mat.mes-hall].[13] J´anos K. Asb´oth and Andrea Alberti, “Spectral flow andglobal topology of the hofstadter butterfly,” Phys. Rev.Lett. , 216801 (2017).[14] F. H. Claro and G. H. Wannier, “Magnetic subbandstructure of electrons in hexagonal lattices,” Phys. Rev.B , 6068–6074 (1979).[15] I Dana, Y Avron, and J Zak, “Quantised hall conduc-tance in a perfect crystal,” Journal of Physics C: SolidState Physics , L679–L683 (1985).[16] Indubala I Satija, Butterfly in the Quantum World , 2053-2571 (Morgan & Claypool Publishers, 2016).[17] Toshikaze Kariyado and Ashvin Vishwanath, “Flatband in twisted bilayer Bravais lattices,” arXiv e-prints, arXiv:1905.02206 (2019), arXiv:1905.02206 [cond-mat.str-el].[18] See Supplemental Material for details.[19] G. H. Wannier, “A result not dependent on ra-tionality for bloch electrons in a magnetic field,”physica status solidi (b) , 757–765 (1978),https://onlinelibrary.wiley.com/doi/pdf/10.1002/pssb.2220880243.[20] P Streda, “Theory of quantised hall conductivity in twodimensions,” Journal of Physics C: Solid State Physics , L717 (1982).[21] J. M. Luttinger, “The effect of a magnetic field on elec-trons in a periodic potential,” Phys. Rev. , 814–817(1951).[22] Takafumi Kita and Masao Arai, “Theory of interact-ing bloch electrons in a magnetic field,” Journal of the Physical Society of Japan , 2813–2830 (2005),https://doi.org/10.1143/JPSJ.74.2813.[23] A. Alexandradinata and Leonid Glazman, “Semiclassi-cal theory of landau levels and magnetic breakdown intopological metals,” Phys. Rev. B , 144422 (2018).[24] Aaron L. Sharpe, Eli J. Fox, Arthur W. Barnard, JoeFinney, Kenji Watanabe, Takashi Taniguchi, M. A.Kastner, and David Goldhaber-Gordon, “Emer-gent ferromagnetism near three-quarters filling intwisted bilayer graphene,” Science , 605–608 (2019),https://science.sciencemag.org/content/365/6453/605.full.pdf.[25] M. Serlin, C. L. Tschirhart, H. Polshyn, Y. Zhang, J. Zhu, K. Watanabe, T. Taniguchi, L. Balents, andA. F. Young, “Intrinsic quantized anomalous hall effectin a moir´e heterostructure,” Science , 900–903 (2020),https://science.sciencemag.org/content/367/6480/900.full.pdf.[26] Z. Song, Z. Wang, W. Shi, G. Li, C. Fang, and B. A.Bernevig, “All “Magic Angles” Are “Stable” Topolog-ical,” ArXiv e-prints (2018), arXiv:1807.10676 [cond-mat.mes-hall].[27] B. ANDREI BERNEVIG and Taylor L. Hughes, Topo-logical Insulators and Topological Superconductors , stu -student edition ed. (Princeton University Press, 2013).[28] Robert Karplus and J. M. Luttinger, “Hall effect in fer-romagnetics,” Phys. Rev. , 1154–1160 (1954). SUPPLEMENTARY MATERIALCONTENTS
Acknowledgments 5References 5Supplementary Material 6I. The algebra in a magnetic field 7II. Basis Completeness and Matrix Elements for Continuum models 7A. Continuum model in real space 7B. Continuum model in real space with orbital-dependent momentum origin shifts 8C. Transforming the model at zero magnetic field into momentum space 9D. The Basis and Hamiltonian in nonzero magnetic field 9E. Numerical Hofstadter calculations: the cutoffs and the momentum space boundary 13F. Number of states per Brillouin zone in a fixed λ sector 14G. Numerical determination of number of occupied states 151. Determination of N K = 0 16III. The example of the TBG continuum model 17A. Description of the model 17B. Spurious zero modes 18IV. Basis Completeness and Matrix Elements for Tight-binding models 19A. Standard Peierls substitution 19B. Nonstandard Peierls substitution 23V. Review of the Diophantine equation 24VI. Quantized Lorentz Susceptibility from s ν I. THE ALGEBRA IN A MAGNETIC FIELD
We first review the algebra obeyed by 2-dimensional (2D) electrons in a uniform static magnetic field B in thecontinuum real space (in the first quantized language). We denote r = ( x, y ) as the position operator, and − i ∇ =( − i (cid:126) ∂ x , − i (cid:126) ∂ y ) as the canonical momentum operator.Assume the magnetic field B corresponds to a gauge field A ( r ) = ( A x ( r ) , A y ( r )), which satisfies ∂ x A y − ∂ y A x = B .We can then define the magnetic length (cid:96) = (cid:112) (cid:126) /eB , where e is the electron charge, and (cid:126) is the Planck constant.The kinematic momentum operator of an electron in a magnetic field is given by Π = (Π x , Π y ), which satisfiesΠ x = − i∂ x − eA x , Π y = − i∂ y − eA y , [Π x , Π y ] = i (cid:126) (cid:96) , (17)or in vector form Π = − i ∇ − A ( r ). In the absence of the gauge field, Π = − i ∇ is the same as the canonicalmomentum. We also define the real space guiding center coordinates R = ( R x , R y ), which satisfy R x = x + (cid:96) (cid:126) Π y , R y = y − (cid:96) (cid:126) Π x , [ R x , R y ] = − i(cid:96) . (18)It can be written in the vector form as R = r − (cid:96) (cid:126) ˆ z × Π . Semiclassically, the guiding center is the central position ofthe cyclotron motion of an electron in magnetic field B . The kinematic momentum operator Π commutes with theguiding center operator R , namely, [Π x , R x ] = [Π x , R y ] = [Π y , R x ] = [Π y , R y ] = 0.For convenience, hereafter we set e = (cid:126) = 1, unless recovery of the original units is needed. Besides, we alwaysunderstand r as the position operator instead of a vector parameter, except that | r (cid:105) stands for a state at position r (where r is a parameter). II. BASIS COMPLETENESS AND MATRIX ELEMENTS FOR CONTINUUM MODELS
In this section, we give the detailed derivation of the basis we choose and the Hamiltonian matrix elements forcontinuum models at zero magnetic field and nonzero magnetic field.
A. Continuum model in real space
We consider a continuum model with M intrinsic orbitals per zero-magnetic-field unit cell. For example, in the one-valley one-spin twisted bilayer graphene (TBG) continuum model in Ref. [2] (see also Sec. III), the intrinsic orbitalsare graphene sublattice and layer indices. In more generic examples with spin-orbit coupling, intrinsic orbitals alsoinclude spin, etc. We denote the lattice Bravais vectors as d and d , and the reciprocal vectors as g and g , whichsatisfy g i · d j = 2 πδ ij , ( i, j = 1 , . (19)We denote the reciprocal lattice as Q = m g + m g , ( m , m ∈ Z ) . (20) In the absence of magnetic field , the continuum model Hamiltonian in a continuum space with a lattice potential isof the generic form H = (cid:90) d r c † α ( r ) (cid:101) H αβ ( r ) c β ( r ) , (21)where (cid:101) H αβ ( r ) = (cid:101) (cid:15) αβ ( − i ∇ ) + (cid:88) j (cid:101) V αβj e i q j · r . (22)Here α, β denote the M intrinsic orbitals, c α ( r ), c † α ( r ) are the electron annihilation and creation operators of orbital α at position r , (cid:101) (cid:15) αβ ( − i ∇ ) is the kinetic term in free space, (cid:101) V αβj is the momentum q j component of the lattice potential.The momentum q j of the lattice potential component (cid:101) V αβj satisfy q j ∈ Q , (23)where Q is the set of reciprocal lattice sites in Eq. (20); thus the lattice potential is periodic with lattice Bravaisvectors d , d . Besides, the Hamiltonian is Hermitian, namely, (cid:101) (cid:15) αβ ( − i ∇ ) = (cid:101) (cid:15) βα ( − i ∇ ) † , and V αβj = V βα ∗ j for momenta q j = − q j .Here we note that H denotes the second quantized Hamiltonian, and (cid:101) H αβ ( r ) denotes the first quantized single-particle Hamiltonian. The basis of the first quantized Hamiltonian (cid:101) H αβ ( r ) is given by c † α ( r ) | (cid:105) , (24)with | (cid:105) being the vacuum state. Since we do not consider interactions, we can work in the first quantized single-particle Hamiltonian hereafter. B. Continuum model in real space with orbital-dependent momentum origin shifts
In writing the Hamiltonian (22) above, the momentum origins of all orbitals are chosen at the Γ point of the BZ ofthe lattice. In some continuum models, it is convenient to shift the momentum origins of different orbitals α to somedesired momenta p α . This is done by transforming the single-particle Hamiltonian (22) from the real space basis inEq. (24) into a new real space basis defined by | r , α (cid:105) = e − i p α · r c † α ( r ) | (cid:105) , (25)where p α is an orbital α dependent momentum vector which can be chosen freely. Here we shall restrict the choicesof p α so that p α = p β if (cid:101) (cid:15) αβ ( − i ∇ ) (cid:54) = 0 , (26)which ensures the kinetic term under the new basis (25) to be a function of − i ∇ only and is independent of r (seeEq. (27)). An example of models with such orbital-dependent momentum origin shifts is the TBG continuum modeloriginally written down in Ref. [2], where the orbitals α of the upper layer have p α = k θ (0 , β of the lower layer have p β = k θ (0 , −
1) (see Sec. III for definition of k θ and more details). Condition (26) is alsosatisfied for the TBG continuum model, since the kinetic term (cid:101) (cid:15) αβ ( − i ∇ ) between an orbital α in the upper layer andan orbital β in the lower layer is zero.Under the assumption (26), we find the first quantized single-particle Hamiltonian transforms under the new basis(25) into H αβ ( r ) = (cid:104) r , α | H | r , β (cid:105) = e i p α · r (cid:101) H αβ ( r ) e − i p β · r = (cid:15) αβ ( − i ∇ ) + (cid:88) j V αβj e i q αβj · r , (27)where we have defined (cid:15) αβ ( − i ∇ ) = (cid:101) (cid:15) αβ ( − i ∇ − p α ) , V αβj = (cid:101) V αβj , q αβj = q j + p α − p β ∈ Q α − Q β ( q j ∈ Q ) , (28)and the momentum lattice Q α for orbital α is defined as Q α = p α + Q = p α + m g + m g , ( m , m ∈ Z ) , (29)where Q is the reciprocal lattice sites in Eq. (20). In particular, we note that with the constraint (26), the transformedkinetic term (cid:15) αβ ( − i ∇ ) in Eq. (28) under the new real space basis (25) is still only a function of − i ∇ , and does notdepend on r . In contrast, if we choose vectors p α which do not satisfy the constraint (26), after the transformation ofEq. (27), we would have a kinetic term (cid:15) αβ ( − i ∇ , r ) = (cid:101) (cid:15) αβ ( − i ∇ − p α ) e i ( p α − p β ) · r that is position r dependent, whichbrings unnecessary complications. Therefore, we will impose the constraint (26).We note that the single-particle Hamiltonian in Eq. (27) is generically invariant up to a unitary transformationunder the translation of a Bravais lattice vector d i . Under the real space basis | r , α (cid:105) in Eq. (25), if we define thetranslation operator of distance d which satisfies T d | r , α (cid:105) = | r + d , α (cid:105) , its real space representation in the first quantizedlanguage is given by T d = e − d ·∇ , and we have H αβ ( r + d i ) = T − d i H αβ ( r ) T d i = e d ·∇ (cid:15) αβ ( − i ∇ ) + (cid:88) j V αβj e i q αβj · r e − d ·∇ = e i p α · d i H αβ ( r ) e − i p β · d i , (30)where d i ( i = 1 ,
2) is a Bravais lattice vector, and we have used the assumption (26), and the relation (19). It is clearthat if we set all p α = 0, we would have the usual translational invariance H αβ ( r + d i ) = H αβ ( r ) which does notinvolve a unitary transformation.For generality, we shall use the first quantized continuum model single-particle Hamiltonian in Eq. (27), which hasorbital-dependent momentum origin shifts p α , and the real space basis is defined in Eq. (25). We note that one canalways choose all p α = 0, in which case the form of the continuum model Hamiltonian reduces back to Eq. (22). C. Transforming the model at zero magnetic field into momentum space
Eq. (27) gives the single-particle Hamiltonian in the absence of magnetic field. It can be written in the momentumspace by Fourier transformation. To do this, we define the momentum space basis | k , Q α , α (cid:105) = 1 √ Ω tot (cid:90) d r e i ( k + Q α ) · r | r , α (cid:105) , (31)where k is the quasimomentum in the first Brillouin zone (BZ), Q α for intrinsic orbital α is defined in Eq. (29), andΩ tot is the total area of the system. Under this momentum space basis | k , Q α , α (cid:105) , the single-particle Hamiltonian H αβ ( r ) in Eq. (27) transforms into the following form diagonal in k (here k , k (cid:48) are in the first BZ): H αβ Q (cid:48) α Q β ( k (cid:48) , k ) = (cid:104) k (cid:48) , Q (cid:48) α , α | H | k , Q β , β (cid:105) = 1Ω tot (cid:90) d r d r (cid:48) e − i ( k (cid:48) + Q (cid:48) α ) · r (cid:104) r , α | H | r (cid:48) , β (cid:105) e i ( k + Q β ) · r (cid:48) = 1Ω tot (cid:90) d r e − i ( k (cid:48) + Q (cid:48) α ) · r H αβ ( r ) e i ( k + Q β ) · r = 1Ω tot (cid:90) d r e − i ( k (cid:48) + Q (cid:48) α ) · r (cid:15) αβ ( − i ∇ ) + (cid:88) j V αβj e i q αβj · r e i ( k + Q β ) · r = δ k , k (cid:48) (cid:15) αβ ( k + Q β ) δ Q (cid:48) α , Q β + (cid:88) j V αβj δ Q (cid:48) α , Q β + q αβj = δ k , k (cid:48) H αβ Q (cid:48) α Q β ( k ) , (32)where we have used the definition of q αβj in Eq. (28), and the condition (26). The Hamiltonian diagonal in k H αβ Q (cid:48) α Q β ( k ) = (cid:15) αβ ( k + Q β ) δ Q (cid:48) α , Q β + (cid:88) j V αβj δ Q (cid:48) α , Q β + q αβj (33)then gives Eq. (2) in the main text. D. The Basis and Hamiltonian in nonzero magnetic field
Under a uniform magnetic field B = B ˆ z , the canonical momentum − i ∇ in Eq. (27) will be replaced by the kineticmomentum Π = − i ∇ − A ( r ), where ∇ × A ( r ) = B . The (first quantized) single-particle Hamiltonian then reads H αβ ( r ) = (cid:15) αβ ( Π ) + (cid:88) j V αβj e i q αβj · r . (34)In the below, we will construct a basis under magnetic field, in which we show that the Hamiltonian in magnetic fieldis block diagonalized into blocks with identical matrix elements (different blocks differ by guiding center translationsin the real space, see explanations below Eq. (55)), and each block is simply given by Eq. (33) (the momentum spaceHamiltonian at zero magnetic field) with the substitution k → ˆ κ + k , ˆ κ = 1 √ (cid:96) ( a + a † , − ia + ia † ) , (35)where k is an arbitrarily chosen momentum vector, a and a † are lowering and raising operators satisfying [ a, a † ] = 1which we will define below, and (cid:96) is the magnetic length.To begin, we would like to find a set of mutually commuting operators to define the quantum numbers of our basis.First, we define a set of Landau level (LL) lowering and raising operators associated with the momentum lattice sites Q α = ( Q α,x , Q α,y ) of orbital α as a Q α = (cid:96) √ x − Q α,x − k ,x + i (Π y − Q α,y − k ,y )] , a † Q α = (cid:96) √ x − Q α,x − k ,x − i (Π y − Q α,y − k ,y )] , (36)0where k = ( k ,x , k ,y ) is a freely chosen fixed momentum which we call the center momentum. They satisfy thecommutation relation [ a Q α , a † Q α ] = 1. We note that for two different sites Q α and Q (cid:48) α , the operators a Q α and a Q (cid:48) α only differ by a constant shift, thus are linearly dependent. We also note that [ a † Q α a Q α , a † Q (cid:48) α a Q (cid:48) α ] (cid:54) = 0 for Q α (cid:54) = Q (cid:48) α .The eigenvalue of a † Q α a Q α for a given Q runs over all the nonnegative integers.Secondly, we define (recall that R = r − (cid:96) (cid:126) ˆ z × Π ) R ˆ τ = R · ˆ τ = R x ˆ τ x + R y ˆ τ y (37)as the guiding center along the ˆ τ direction, where ˆ τ = (ˆ τ x , ˆ τ y ) is a unit vector so chosen that ˆ τ · (ˆ z × g )ˆ τ · (ˆ z × g ) is an irrationalnumber (the reason will be explained below Eq. (41)). In an infinite real space (which we assume is the case here),the eigenvalue of R ˆ τ runs over all real numbers R . This can be seen from the fact that one can shift the eigenvalueof R ˆ τ by any value b ∈ R using the operator e ib R · (ˆ τ × ˆ z ) /(cid:96) , namely, e ib R · (ˆ τ × ˆ z ) /(cid:96) R ˆ τ e − ib R · (ˆ τ × ˆ z ) /(cid:96) = R ˆ τ + b , ( b ∈ R ) (38)We note that once R ˆ τ is diagonalized for a given direction ˆ τ , one cannot further diagonalize the guiding centercoordinate R ˆ τ (cid:48) = R · ˆ τ (cid:48) along any other direction ˆ τ (cid:48) (cid:54) = ± ˆ τ , since [ R ˆ τ (cid:48) , R ˆ τ ] = − i(cid:96) ˆ z · ( ˆ τ (cid:48) × ˆ τ ) (cid:54) = 0 unless ˆ τ (cid:48) = ± ˆ τ (Note that R ˆ τ = − R − ˆ τ , so diagonalization of R ˆ τ is sufficient).Since [ R , Π ] = 0, it is easy to see that[ R ˆ τ , a Q α ] = [ R ˆ τ , a † Q α ] = [ R ˆ τ , a † Q α a Q α ] = 0 (39)for any Q α . Therefore, R ˆ τ and a † Q α a Q α (for a given Q α ) can be diagonalized simultaneously. For a given Q α , thereis no other functions of R , Π independent of R ˆ τ and a † Q α a Q α which commute with both R ˆ τ and a † Q α a Q α , so R ˆ τ and a † Q α a Q α form the maximal commuting set of operators R , Π . We shall therefore use R ˆ τ and a † Q α a Q α to define thequantum numbers of our basis.In the previous Hofstadter method for continuum models [5], the chosen basis makes the lattice potential term e i q αβj · r in Eq. (34) a complicated dense matrix. We hope to find a different basis where the operator e i q αβj · r havesimple matrix elements, so that the Hamiltonian takes a much simpler form. Now we describe the construction ofsuch a basis (in which the operator e i q αβj · r has very simple matrix elements, see Eq. (48)). For each momentum site Q α , we define a set of basis | λ, Q α , n, α (cid:105) satisfying R ˆ τ | λ, Q α , n, α (cid:105) = [ λ − (cid:96) ˆ τ · (ˆ z × Q α )] | λ, Q α , n, α (cid:105) , a † Q α a Q α | λ, Q α , n, α (cid:105) = n | λ, Q α , n, α (cid:105) , (40)where n ≥ λ is a real number, and α is the intrinsic orbital index.Note that we have defined the eigenvalue of R ˆ τ in a Q α dependent way (we could have defined the R ˆ τ eigenvaluewithout (cid:96) ˆ τ · (ˆ z × Q α ), which would be equivalent to the basis used in usual method [5]). One advantage of sucha definition is that, the operator e i q αβj · r will not change the number λ when acting on the basis | λ, Q α , n, α (cid:105) , aswe will prove in Eq. (48). Another advantage of such a definition in Eq. (40) is that, for a fixed number λ , thebasis | λ, Q α , n, α (cid:105) of all quantum numbers α , Q α and n are orthonormal. To see this, consider two momentum sites Q α , Q (cid:48) α ∈ p α + g Z + g Z (as defined in Eq. (29)), and assume Q α − Q (cid:48) α = m g + m g . If two states | λ, Q α , n, α (cid:105) and | λ, Q (cid:48) α , n (cid:48) , α (cid:105) have equal R ˆ τ eigenvalues, namely, λ − (cid:96) ˆ τ · (ˆ z × Q α ) = λ − (cid:96) ˆ τ · (ˆ z × Q (cid:48) α ) → m ˆ τ · (ˆ z × g ) + m ˆ τ · (ˆ z × g ) = 0 , (41)we must have m = m = 0, namely, Q (cid:48) α = Q α , since ˆ τ · (ˆ z × g )ˆ τ · (ˆ z × g ) is an irrational number. Vice versa, if Q (cid:48) α (cid:54) = Q α ,the two states | λ, Q α , n, α (cid:105) and | λ, Q (cid:48) α , n (cid:48) , α (cid:105) will have different R ˆ τ eigenvalues and will be orthogonal to each other.Therefore, the basis with a fixed number λ satisfies the orthonormal relation: (cid:104) λ, Q (cid:48) β , n (cid:48) , β | λ, Q α , n, α (cid:105) = δ βα δ Q (cid:48) β Q α δ n (cid:48) n . (42)We note that Eq. (42) would not hold if we define R ˆ τ | λ, Q α , n, α (cid:105) = λ | λ, Q α , n, α (cid:105) instead, as two different Q ’s wouldhave the same eigenvalue λ .As we will show later below Eq. (48), the way of defining the basis in Eq. (40) will greatly simplify the matrixelements of Hamiltonian (34). However, before we move on, we note that the set of basis with a fixed number λ (satisfying Eq. (42)) is not a complete basis for the Hamiltonian (34), since for a fixed λ , the R ˆ τ eigenvalue1 λ − (cid:96) ˆ τ · (ˆ z × Q α ) (of all momentum sites Q α ) does not run over the entire real number set R (the complete set of R ˆ τ eigenvalues is R , see the argument above Eq. (38)). Therefore, we need to allow λ to run over a certain set tomake the basis | λ, Q α , n, α (cid:105) a complete basis. We now prove this can be done by assuming λ runs over the followingquotient set: λ ∈ Λ ˆ τ = R / [ (cid:96) ˆ τ · (ˆ z × g ) Z + (cid:96) ˆ τ · (ˆ z × g ) Z ] , (43)namely, each number λ labels a coset { λ + (cid:96) ˆ τ · (ˆ z × g ) Z + (cid:96) ˆ τ · (ˆ z × g ) Z } of the subgroup (cid:96) ˆ τ · (ˆ z × g ) Z + (cid:96) ˆ τ · (ˆ z × g ) Z inthe real number group R . Two numbers λ and λ (cid:48) label the same coset if they belongs to the same coset. For example,assume we have (cid:96) ˆ τ · (ˆ z × g ) = 1 and (cid:96) ˆ τ · (ˆ z × g ) = √
2, then the quantum number λ = 0 will label the coset { Z + √ Z } . Accordingly, λ = 0 is identical to λ = 1 and λ = √
2, etc. In contrast, λ = 0 and λ = √ λ in each coset to represent the coset,so each coset is represented by a definite number λ . For example, we can choose to use the definite number λ = 0(instead of 1, √
2, etc) to represent the coset { Z + √ Z } . In this way, we can represent each element (coset) in theset Λ ˆ τ in Eq. (43) by a definite number λ .With the set of representative numbers λ given in Eq. (43), we now show that the basis | λ, Q α , n, α (cid:105) forms acomplete orthonormal basis. Note that for each λ in Eq. (43), the R ˆ τ eigenvalue λ − (cid:96) ˆ τ · (ˆ z × Q α ) of all Q α runs overall the numbers in the coset { λ − (cid:96) ˆ τ · (ˆ z × p α ) + (cid:96) ˆ τ · (ˆ z × g ) Z + (cid:96) ˆ τ · (ˆ z × g ) Z } , where p α is the momentum originshift of orbital α defined in Eq. (25), and recall that Q α = p α + g Z + g Z as defined in Eq. (29). Therefore, the R ˆ τ eigenvalue of | λ, Q α , n, α (cid:105) of all λ in Eq. (43) and all Q α runs over all the real numbers R (since λ is a continuousvariable in the quotient set (43)). In particular, two orbital α states | λ, Q α , n, α (cid:105) and | λ (cid:48) , Q (cid:48) α , n (cid:48) , α (cid:105) will have their R ˆ τ eigenvalues in different cosets if λ (cid:54) = λ (cid:48) ( λ, λ (cid:48) ∈ Λ ˆ τ in Eq. (43), note that we have chosen a unique definite number λ in each coset to represent the coset in Eq. (43)). Therefore, we have the orthonormal relation (cid:104) λ (cid:48) , Q (cid:48) β , n (cid:48) , β | λ, Q α , n, α (cid:105) = δ βα δ λλ (cid:48) δ Q (cid:48) β Q α δ n (cid:48) n . (44)Thus, the basis | λ, Q α , n, α (cid:105) with λ defined in Eq. (43) forms a complete orthonormal basis of Hamiltonian (40).We now show that the most important advantage of the basis | λ, Q α , n, α (cid:105) satisfying Eq. (44) is, the Hamiltonian(34) is diagonal in λ , and its matrix elements are independent of λ . To see this, we first note that from the commutationrelations in Sec. (I), we have e − i q · r R ˆ τ e i q · r = R ˆ τ − (cid:96) ˆ τ · (ˆ z × q ) , e − i q · r a Q α e i q · r = a Q α + (cid:96) √ q x + iq y ) = a Q α − q , (45)where r is the position operator, and q is an arbitrary momentum vector. Therefore, we have R ˆ τ e i q · r | λ, Q α , n, α (cid:105) = e i q · r ( R ˆ τ − (cid:96) ˆ τ · (ˆ z × q )) | λ, Q α , n, α (cid:105) = [ λ − (cid:96) ˆ τ · (ˆ z × ( Q α + q ))] e i q · r | λ, Q α , n, α (cid:105) (46) a † Q α + q a Q α + q e i q · r | λ, Q α , n, α (cid:105) = e i q · r a † Q α a Q α | λ, Q α , n, α (cid:105) = ne i q · r | λ, Q α , n, α (cid:105) , (47)which indicates e i q · r | λ, Q α , n, α (cid:105) = | λ, Q α + q , n, α (cid:105) , (48)where r is the position operator (understood as (cid:80) α (cid:82) r | r , α (cid:105)(cid:104) r , α | d r when expanded in position basis), and it alsoplays the role of the generator of momentum translations.We then examine the matrix elements of Hamiltonian (34). Using Eq. (48), we find the lattice potential term V αβj e i q αβj · r in the Hamiltonian (34) has matrix elements (cid:104) λ (cid:48) , Q (cid:48) α , m, α | V αβj e i q αβj · r | λ, Q β , n, β (cid:105) = V αβj (cid:104) λ (cid:48) , Q (cid:48) α , m, α | λ, Q β + q αβj , n, β (cid:105) = V αβj δ λ (cid:48) λ δ Q (cid:48) α , Q β + q αβj δ mn , (49)where r on the left hand side is understood as the position operator (instead of a number). Therefore, the latticehopping potential term V αβj e i q αβj · r is diagonal in quantum numbers λ and n when acting on basis | λ, Q α , n, α (cid:105) , whilechanges the orbital from β to α , and shifts the reciprocal momentum Q β to Q β + q αβj .As for the kinetic term (cid:15) ( Π ), we first note that the kinetic momentum Π commutes with R ˆ τ and does not changethe R ˆ τ eigenvalue (which has one-to-one correspondence with the pair of quantum numbers λ and Q α ). Therefore,2the matrix elements of (cid:15) ( Π ) has to be diagonal in both λ and Q α . Besides, for each orbital α , we can rewrite thekinetic momentum as Π = ˆ κ Q α + k + Q α , where we have defined a Q α dependent operatorˆ κ Q α = 1 √ (cid:96) ( a Q α + a † Q α , − ia Q α + ia † Q α ) . (50)Note that ˆ κ Q α only acts on the quantum number n of basis | λ, Q α , n, α (cid:105) , which obeys the following rules: a Q α | λ, Q α , n, α (cid:105) = √ n | λ, Q α , n − , α (cid:105) , a † Q α | λ, Q α , n, α (cid:105) = √ n + 1 | λ, Q α , n + 1 , α (cid:105) . (51)Therefore, it is easy to find that (cid:104) λ (cid:48) , Q (cid:48) β , m, β | (cid:15) ( Π ) | λ, Q α , n, α (cid:105) = δ λ (cid:48) λ δ Q (cid:48) β , Q α (cid:104) λ (cid:48) , Q (cid:48) β , m, β | (cid:15) (ˆ κ Q α + k + Q α ) | λ, Q α , n, α (cid:105) = δ λ (cid:48) λ δ Q (cid:48) β , Q α [ (cid:15) βα (ˆ κ Q α + k + Q α )] mn . (52)Mathematically, we can define an operator ˆ κ = 1 √ (cid:96) ( a + a † , − ia + ia † ) (53)without the Q α subindex, where a and a † are some lowering and raising operators satisfying [ a, a † ] = 1 (which neednot have any relation with a Q α and a † Q α ). It is then easy to see that the matrix element in Eq. (52) is mathematicallyequal to [ (cid:15) βα (ˆ κ Q α + k + Q α )] mn = [ (cid:15) βα (ˆ κ + k + Q α )] mn = (cid:104) m | (cid:15) βα (ˆ κ + k + Q α ) | n (cid:105) , (54)where | n (cid:105) is the basis of a and a † defined by a | n (cid:105) = √ n | n − (cid:105) and a † | n (cid:105) = √ n + 1 | n + 1 (cid:105) . Therefore, mathematicallyone could replace ˆ κ Q α by ˆ κ without ambiguity, provided that one remembers that ˆ κ acts on the quantum number n .From Eqs. (49), (52) and (54), it is easy to see that the matrix elements of both the kinetic term and the latticepotential term of Hamiltonian (34) are diagonal in λ and independent of λ . Therefore, we can divide the entire Hilbertspace into subspaces with different (continuous) quantum number λ ; the energy spectra of all the λ subspaces are thesame. Within the subspace of a fixed λ , the Hamiltonian matrix element from basis | λ, Q β , n, β (cid:105) to | λ, Q (cid:48) α , m, α (cid:105) canbe written as H λ,αβ Q (cid:48) α Q β ,mn = (cid:2) (cid:15) αβ (ˆ κ + k + Q β ) (cid:3) mn δ Q (cid:48) α , Q β + (cid:88) j V αβj δ Q (cid:48) α , Q β + q αβj δ mn , (55)where [ (cid:15) αβ (ˆ κ + k + Q β )] mn is defined in Eq. (54). This is nothing but the zero magnetic field momentum spaceHamiltonian (33) with the substitution k → ˆ κ + k (given also in Eq. (5) of the main text). It is then sufficient tocompute the spectrum within just one fixed λ subspace.We note that different λ sectors have different eigenstate wave functions, although they have identical Hamiltonianmatrix elements independent of λ as given by Eq. (55). To see this explicitly, if we take a Landau gauge perpendicularto the ˆ τ direction, A = B ( r · ˆ τ )(ˆ z × ˆ τ ), we have (cid:104) r , α | λ, Q α , n, α (cid:105) = e i [ λ − (cid:96) ˆ τ · (ˆ z × Q α )][ r · (ˆ z × ˆ τ )] /(cid:96) − [ r · ˆ τ − λ + (cid:96) ˆ τ · (ˆ z × Q α )] / (cid:96) h n (cid:16) (cid:96) − [ r · ˆ τ − λ + (cid:96) ˆ τ · (ˆ z × Q α )] (cid:17) , (56)where h n ( x ) is the n -th Hermite polynomial. Therefore, one can see explicitly that the wave functions | λ, Q α , n, α (cid:105) atdifferent λ have different guiding centers. Assume the subspace Hamiltonian H λ,αβ Q (cid:48) α Q β ,mn in sector λ has an eigenstate | λ, u (cid:105) = (cid:80) α, Q α ,n u α, Q α ,n | λ, Q α , n, α (cid:105) at energy energy E u , where the coefficients u α, Q α ,n are independent of λ . Thenthe subspace Hamiltonian H λ (cid:48) ,αβ Q (cid:48) α Q β ,mn in sector λ (cid:48) will have an eigenstate | λ (cid:48) , u (cid:105) = (cid:80) α, Q α ,n u α, Q α ,n | λ (cid:48) , Q α , n, α (cid:105) atthe same energy E u . However, the ˆ τ direction guiding center coordinate R ˆ τ of the two wave functions | λ, u (cid:105) and | λ (cid:48) , u (cid:105) will differ by λ (cid:48) − λ , namely, the central position of the two wave functions are different.We also note that, if there is a disorder potential that breaks the periodicity of the lattice model, e.g., a potentialterm δ (cid:101) V αβ e i (cid:101) q · r from orbitals β to α with (cid:101) q (cid:54) = Q + p α − p β for any reciprocal vector Q , this term will couple different λ sectors, and the Hamiltonian will no longer be diagonal in λ and independent of λ . By Eq. (48) and Eq. (40),such a term δ (cid:101) V αβ e i (cid:101) q · r will couple the sector of coset labeled by λ (defined in Eq. (43)) with the sector of coset of λ − (cid:96) ˆ τ · (ˆ z × [ (cid:101) q + p β − p α )], which does not live in the same coset. In this paper, we shall not consider any disorderpotential breaking the periodicity of the lattice model.In particular, if the kinetic energy (cid:15) ( k ) is a polynomial of k up to power ∆ (∆ ∈ Z + ), the matrix element in Eq.(55) has to be zero for | m − n | > ∆. In this case, the Hamiltonian (55) is sparse.3 E. Numerical Hofstadter calculations: the cutoffs and the momentum space boundary
In numerical calculations, one needs to take a cutoff in the reciprocal lattice at a boundary enclosing an area of N Q BZs, and a cutoff in the LL quantum number n ≤ N L . The Hamiltonian (55) is then a matrix of size M N L N Q . Wenow explain how the cutoffs N Q and N L set a momentum space boundary for Hamiltonian (55).For concreteness, assume the cutoff of the reciprocal lattice encloses a circular momentum space area N Q Ω BZ centered at the center momentum k , where Ω BZ = | g × g | = 4 π / | d × d | is the BZ area. This restricts thereciprocal momentum sites within a momentum radius | k + Q α | ≤ (cid:113) N Q Ω BZ /π . (57)In addition, by Eq. (53), with the LL cutoff N L , for any states we have | ˆ κ | = (cid:104) ˆ κ (cid:105) = 1 (cid:96) (cid:104) (2 a † a + 1) (cid:105) ≤ N L + 1 (cid:96) ≈ N L (cid:96) , (58)where we have defined | ˆ κ | = (cid:112) (cid:104) ˆ κ (cid:105) as the norm of the operator ˆ κ . Therefore, ˆ κ is restricted within a radius | ˆ κ | ≤ √ N L /(cid:96) .The Hamiltonian (55) can be viewed as a lattice model in the momentum space with “hoppings” V j between nearbyreciprocal sites, and an “on-site potential energy” (cid:15) (ˆ κ + k + Q α ) on each site Q α , and ˆ κ plays the role of the“position” operator in the momentum space (although its x and y components are noncommuting). The momentumspace probability (norm square of amplitude) of the basis wave function | λ, Q α , n, α (cid:105) is concentrated circularly neara ring of radius (cid:112) (cid:104) ˆ κ (cid:105) ≈ √ n/(cid:96) .The radius cutoff of k + Q α and the radius cutoff of ˆ κ become equal when (cid:112) N L /(cid:96) = (cid:113) N Q Ω BZ /π , → ϕ π = B | d × d | π = 2 π(cid:96) Ω BZ = N Q N L , (59)where ϕ = B | d × d | is the magnetic flux per unit cell, and we have used the Brillouin zone area Ω BZ = 4 π / | d × d | and B = 1 /(cid:96) . Given N Q , N L which we pick in our calculation, this flux ϕ/ π = N Q /N L then separates the systeminto a small B regime and a large B regime as follows.If √ N L /(cid:96) < (cid:112) N Q Ω BZ /π , or equivalently ϕ/ π < N Q /N L (the small B field regime), the momentum space“position” ˆ κ has a hard cutoff | ˆ κ | ≤ √ N L /(cid:96) which serves as the momentum space boundary.On the contrary, if √ N L /(cid:96) > (cid:112) N Q Ω BZ /π , or equivalently ϕ/ π > N Q /N L (the large B field regime), we havethe expectation value of | ˆ κ + k + Q α | > √ N L /(cid:96) − (cid:112) N Q Ω BZ /π for all sites Q α if a state has expectation value | ˆ κ | > (cid:112) N Q Ω BZ /π . In general, the kinetic energy, or “on-site potential energy” in the momentum space (cid:15) (ˆ κ + k + Q α ),is an increasing function of | ˆ κ + k + Q α | . If an eigenstate has a large expectation value of | ˆ κ + k + Q α | for any Q α (because it has expectation value | ˆ κ | > (cid:112) N Q Ω BZ /π ), it will also have a large expectation value of kinetic energy (cid:15) (ˆ κ + k + Q α ) for any Q α , thus its eigenenergy is expected to be large, and cannot be a reliable eigenstate of thelow energy Hofstadter bulk bands.Therefore, from the reasoning the above, we can define a momentum boundary radius given by the cutoffs N L and N Q as κ b ≈ min { √ N L (cid:96) , (cid:114) N Q Ω BZ π } , (60)and a trustable low-energy eigenstate in the Hofstadter bulk bands should have its expectation value | ˆ κ | < κ b . Anystate with expectation value | ˆ κ | (cid:38) κ b are effectively localized on the momentum space boundary at radius κ b (if weview ˆ κ as a momentum space coordinate), which we will call the momentum space edge states . Accordingly, we callthe states with | ˆ κ | < κ b the momentum space bulk states . The momentum space bulk states give the Hofstadter bulkband spectra we want to calculate.As we discussed in the main text (below main text Eq. (10)), the momentum space edge states can be identifiedby a boundary projection operator P κ b ,w , the matrix elements of which are defined by (cid:104) λ (cid:48) , Q (cid:48) α , m, β | P κ b ,w | λ, Q β , n, α (cid:105) = Θ (cid:16) n − ( κ b − w ) (cid:96) / (cid:17) δ λ (cid:48) λ δ mn δ αβ δ Q (cid:48) α , Q β , (61)where Θ( x ) is the Heaviside unit step function, κ b ≈ min { √ N L (cid:96) , (cid:113) N Q Ω BZ π } as defined in Eq. (60), and w > (cid:104) P κ b ,w (cid:105) will be mainly concentrated at | ˆ κ | > κ b − w in the momentum space, namely, within distance w inside theboundary radius κ b , thus are momentum space edge states. These edge eigenstates with large (cid:104) P κ b ,w (cid:105) > P c abovea certain chosen threshold P c ∈ [0 ,
1] can then be deleted in the energy spectrum, so that only the bulk spectrumHofstadter butterfly is kept (see main text Fig. 2(b), which is calculated with cutoffs N Q = 37 and N L = 60). Inpractice, one may properly adjust w and (cid:104) P κ b ,w (cid:105) for different magnetic fields B and different energy range to reach acleaner Hofstadter butterfly. Generically, the momentum space edge states in a smaller (larger) Hofstadter gap will beless (more) localized at the momentum boundary, thus requires a larger (smaller) edge width w and a smaller (larger)projection threshold P c . In calculating main text Fig. 2(b), we have chosen w = min { (cid:96) − , . √ Ω BZ } , and the edgeprojection threshold P c = 0 .
5. Here the factor 1 . √ Ω BZ is numerically tested to be a good choice foreliminating most edge states in the main Hofstadter gaps of the TBG model calculated here. Generically the optimalorder 1 factor in front of √ Ω BZ depends on the models calculated (and can even be chosen to be dependent on themagnetic field B and the energy of the eigenstate). Generically, we suggest to choose w to be around the order of thesmaller one of (cid:96) − and √ Ω BZ .Besides, we note that the Hofstadter butterfly with edge states deleted may have remaining edge state “hairs”near the edges of the Hofstadter gaps (i.e., not a clean Hofstadter butterfly), as shown in Fig. 2(b). This is becausethe momentum space edge states approaching the Hofstadter gap edges are more and more delocalized from themomentum space boundary, thus some of such edge states cannot satisfy the criteria (cid:104) P κ b ,w (cid:105) > P c and thus cannot bedeleted. By making the momentum boundary radius κ b larger (which is computationally more expensive) and choosea larger edge width w , one can reduce such edge state “hairs” and improve the clearness of the Hofstadter butterfly. F. Number of states per Brillouin zone in a fixed λ sector Here we discuss the number of occupied states per Brillouin zone (i.e., per reciprocal lattice “unit cell”) ρ K of theHamiltonian H λ,αβ Q (cid:48) α Q β ,mn in a fixed λ sector (coset) given a Fermi energy, which we used in the main text Eq. (8).Consider a magnetic field B corresponding to flux per unit cell ϕ = B Ω = (cid:96) − Ω, where Ω = | d × d | is thezero-magnetic-field unit cell area. Assume the Fermi energy is (cid:15) F , and the number of occupied states (below theFermi energy (cid:15) F ) per zero-magnetic-field unit cell area Ω in the real space is ρ .To find out the number of occupied states in each λ sector, we first count how many λ sectors there are. In aninfinite real space and with an infinite reciprocal lattice Q (i.e., without reciprocal lattice cutoff), the set of λ in Eq.(43) is an infinite set, the size of which cannot be perceived easily. To make the set of λ finite, we assume the systemhas a large but finite real space area Ω tot = N Ω Ω (where N Ω is the number of zero-magnetic-field unit cells), andtake a finite reciprocal lattice number cutoff N Q (as we did in the numerical calculation described in Sec. II E), butwe keep the Landau level number cutoff N L → ∞ . By Eq. (40), for each orbital α , each λ sector consists of thesub-Hilbert spaces with R ˆ τ eigenvalues λ − (cid:96) ˆ τ · (ˆ z × Q α ) with Q α = Q + p α (defined in Eq. (29)) running over all Q within the reciprocal lattice cutoff N Q . This means each λ sector consists of the sub-Hilbert spaces of N Q different(discrete) eigenvalues of R ˆ τ . On the other hand, for each orbital α , the total number of (discrete) R ˆ τ eigenvalues inthe system is equal to the degeneracy of a single Landau level in the free space, which is N R = Ω tot π(cid:96) . (62)This is because in the free space, the Landau level states of a definite LL band can be completely labeled by theeigenvalue of the guiding center along a certain direction (here chosen to be R ˆ τ ). Since each λ sector allows N Q different eigenvalues of R ˆ τ , and different λ sectors are orthogonal to each other in the Hilbert space, we conclude thenumber of λ sectors (cosets) in our method are related by N λ = N R N Q = Ω tot π(cid:96) N Q . (63)Meanwhile, since there are ρ occupied states per zero-magnetic-field unit cell area Ω = | d × d | , the total numberof occupied states is given by N occ = ρ Ω tot Ω . (64)Since the number of R ˆ τ eigenvalues is N R , we can define the number of occupied states in each R ˆ τ eigenvalue sectoras ρ K = N occ N R = ρ π(cid:96) Ω = 2 πϕ ρ . (65)5Note that ρ K is nothing but the LL filling fraction. On the other hand, the occupied states should evenly belong toeach λ sector, since different λ sectors have the same spectrum. Therefore, we conclude the number of states occupiedin each λ sector (in the N L → ∞ and finite N Q case considered here, i.e., N Q /N L →
0) is N ( λ ) K = N occ N λ = N occ N R N Q = 2 πϕ ρN Q = ρ K N Q , (66)which is independent of λ . Since N Q is the reciprocal lattice cutoff, or the number of Brillouin zones (BZs) we keepin the reciprocal lattice (note that we assumed N L → ∞ which does not give a momentum space cutoff), we see that ρ K = (2 π/ϕ ) ρ can be understood as the number of occupied states per BZ in a fixed λ sector.At rational flux ϕ = 2 πp/q , according to the Diophantine equation (109), we have ρ = ν/q where ν is the numberof occupied magnetic bands in the magnetic unit cell, and thus we find the number of occupied states per BZ (with λ fixed) is ρ K = (2 π/ϕ ) ρ = ν/p . Accordingly, the Diophantine equation (109) can be rewritten as t ν + s ν πϕ = ρ K , (67)i.e., main text Eq. (8).Since the total number of occupied states N ( λ ) K = N Q ρ K in a λ sector is independent of λ , hereafter we simplydenote it as N K , as appears in the main text Eqs. (9) and (10). G. Numerical determination of number of occupied states
This subsection discusses how to determine the number of occupied states N K in a gap in numerical calculations.We note that N K here is counted below the mid-gap energy of a gap, which may disperse as a function of magneticfield B (e.g., the inclined dashed lines in the main text Fig. 2(c)-(d)). Generically, N K counted in this way willslightly depend on where the mid-gap energy is chosen, which determines how many in-gap edge states are includedin N K in addition to the bulk states. However, the relative error in counting N K will tend to zero as N L and N Q increase, since the ratio between the number of in-gap edge states and the number of bulk band states will tend tozero when the momentum space area increases.In the main text, we have shown that in the calculation with LL cutoff N L and reciprocal cutoff N Q (in a fixed λ sector), there are two regimes:i) The regime ϕ/ π < N Q /N L , for which the momentum space boundary is at radius κ b = √ N L /(cid:96) , and themomentum space bulk area is A K = 2 πN L /(cid:96) , so the number of occupied states N K = ρ K A K / Ω BZ in the fixed λ sector (note that this is different from Eq. (66) where ϕ/ π > N Q /N L ) satisfies Eq. (9) in the main text, namely, N K = N L ( t ν ϕ/ π + s ν ) , (68)where t ν is the Chern number of the gap, and s ν is another integer which we show in Sec. VI can be understood asa dual Chern number for the momentum space. Accordingly, the in-gap spectral flow rate in a gap is given byd N K d( ϕ/ π ) = N L t ν , (69)which allows us to determine t ν of a gap in this regime by counting the number of states flowing across the midgapenergy per flux number ϕ/ π , as described in the example given below main text Eq. (9) (shown in main text Fig.2c). Further, if the number of occupied states N K in the gap at some flux ϕ is known (counted relative to somereference point, which will be discussed below), we could also derive s ν of the gap from Eq. (68).ii) The regime ϕ/ π > N Q /N L , for which the momentum space boundary is at radius κ b = (cid:112) N Q Ω BZ /π , and themomentum space bulk area is A K = N Q Ω BZ , so the number of occupied states N K = ρ K A K / Ω BZ in the fixed λ sector satisfies main text Eq. (10), namely, N K = N Q ( t ν + 2 πs ν /ϕ ) . (70)Accordingly, the in-gap spectral flow rate in a gap is given byd N K d(2 π/ϕ ) = N Q s ν , (71)6which allows us to determine s ν of a gap in this regime by counting the number of states (without deleting momentumspace edge states by the projector method in Sec. II E) flowing across the midgap energy per inverse flux number2 π/ϕ , as described in the example given below main text Eq. (10) (shown in main text Fig. 2d). Further, if thenumber of occupied states N K in the gap at some flux ϕ is known, we could also derive t ν of the gap from Eq. (70).In either regime, finding out N K at some ϕ in the numerical results allows us to derive both t ν and s ν of a gap. Innumerical calculations, however, the number of occupied states N K below a gap should be counted relative to somereference energy level where N K = 0 is defined. Therefore, we need to find out where the N K = 0 energy level isdefined, which depends on models. This is discussed below.
1. Determination of N K = 0 A generic method to find out N K = 0 is the following: first, find a gap which extends over both the small ϕ regime ϕ/ π < N Q /N L and the large ϕ regime ϕ/ π > N Q /N L (this is the condition for this generic method to work), whichwe call the reference gap. By extracting out the spectral flow rate in the small and large ϕ regimes, one can derive t ν and s ν of the reference gap from Eqs. (69) and (71), respectively. Then, N K of this reference gap at any ϕ can bedetermined from Eqs. (68) and (70). Then, the number of occupied states N K in any other gap j at a given ϕ canbe determined in the numerical calculation by counting the number of states from the midgap energy of the referencegap to the midgap energy of gap j , plus the number of occupied states in the reference gap which is known. Fromthe reference gap, one could determine which energy level corresponds to N K = 0.As an example, in the Hofstadter butterfly of the TBG model in the main text Fig. 2(a) and 2(b) (see Sec. III fordetails), we can take the (1 ,
0) gap (labeled in main text Fig. 2(b)) as such a reference gap. This Hofstadter spectrum iscalculated by setting N Q = 37 and N L = 60. In the regime ϕ/ π < N Q /N L , the spectral flow rate across the mid-gapenergy of gap (1 ,
0) can be counted along the dashed line in Fig. 2(c) to be d N K /d ( ϕ/ π ) ≈ / .
25 = 64 = N L t ν (16levels when ϕ/ π increases from 0.25 to 0.5), so we find t ν = 1 being the integer closest to 64 /N L . Then in the regime ϕ/ π > N Q /N L , there are no levels flowing in the (1 ,
0) gap except for some horizontal lines, which are spurious Diraczero modes of the TBG model and should not be counted in N K (see Sec. III B for detailed explanation). Therefore,the spectral flow rate in this regime is d N K /d (2 π/ϕ ) = 0 = N Q s ν , which leads to s ν = 0. Thus the gap is labeled byquantum numbers ( t ν , s ν ) = (1 , N K = N L (1 × . ϕ/ π = 0 . < N Q /N L . Therefore, the N K = 0 level can be found by counting 30 levels downwards from the midgapenergy of the (1 ,
0) gap at flux ϕ/ π = 0 .
5, which turns out to be approximately the level at zero energy (see themain text Fig. 2(c)).We further discuss the following two special cases, where the reference point N K = 0 can be determined moreeasily:i) Models with a kinetic energy bounded from below, e.g., a single-orbital model with a Hamiltonian (cid:101) H ( r ) = (cid:15) ( − i ∇ ) + (cid:80) j V j e i q j · r with a quadratic kinetic energy (cid:15) ( k ) = k / m ≥
0, where m is the electron effective mass. Inthis case, one must have Chern number t ν = 0 below the lowest energy band of the entire spectrum, i.e., when nostates are occupied at all. By Eq. (70), one then finds N K = 0 below the lowest energy band in the limit ϕ → ∞ .Since there are no states at lower energies, there should be no spectra flows below the lowest band with respect to ϕ ,so one should have N K = 0 below the lowest energy band at any flux ϕ .ii) Models with a Dirac kinetic energy (which has no lower bound), such as the TBG model in Eq. (73) which has akinetic term (cid:15) ( k ) = v F σ ∗ · k (where σ ∗ = ( σ x , − σ y )). In this case, if there is no LL cutoff and reciprocal lattice cutoff,the energy spectrum of the system does not have a lower bound. With a LL cutoff N L and a reciprocal lattice cutoff N Q , the Hamiltonian size is M N L N Q for M intrinsic orbitals (each Dirac kinetic term (cid:15) ( k ) = v F σ ∗ · k has two intrinsicorbitals). In the ϕ → ∞ limit, the energies of the eigenstates are dominated by the kinetic term (cid:15) ( k ) = v F σ ∗ · k (since k is replaced by ˆ κ + k and ˆ κ ∝ (cid:96) − = √ B which goes to infinity as ϕ → ∞ ), which should give a (nearly)particle-hole symmetric spectrum because of the particle-hole symmetry of the Dirac kinetic term. Accordingly, allthe spectral flows should be (nearly) particle-hole symmetric about the half filling, which fixes N K = 0 at the halffilling point. Therefore, in the limit ϕ → ∞ , one has N K = 0 at the half filling of the Hamiltonian with cutoffs N L and N Q , e.g., the filling between the M N L N Q / M N L N Q / M N L N Q does not change with respect to ϕ (for fixed N L and N Q ), and N K = 0 is a referencefilling independent of ϕ , we conclude that N K = 0 is between the M N L N Q / M N L N Q /
2) + 1-thlevel for any ϕ . In particular, for the TBG model in Eq. (73), the energy spectrum is particle-hole symmetric at all ϕ , so the half filling point N K = 0 is at zero energy at any ϕ . We note that when further counting the number ofoccupied states N K of certain gaps relative to this half filling N K = 0 point, one needs to exclude the unphysicalspurious modes due to LL cutoff N L as discussed in Sec. III B.7 III. THE EXAMPLE OF THE TBG CONTINUUM MODEL
In this section, we explain the Hofstadter spectrum calculation of the one-valley TBG continuum model in Ref. [2],the results of which are given in the main text Fig. 2.
A. Description of the model
The model consists of the Dirac electrons of the same valley of two graphene layers, which are relatively twisted byangle θ . Besides, we only consider one spin, namely, the model is spinless. The model can be written in real space asa 4 × H TBG = (cid:32) − iv F σ ∗ · ∇ (cid:80) j =1 V j e i q j · r (cid:80) j =1 V † j e − i q j · r − iv F σ ∗ · ∇ (cid:33) , (72)where the upper (lower) two basis are the A and B sublattices of the upper (lower) monolayer graphene, σ ∗ = ( σ x , − σ y )are the Pauli matrices acting on A and B sublattices of the monolayer graphene lattice, the momenta q j are given by q = k θ (0 , − T , q = k θ (cid:32) √ , (cid:33) T , q = k θ (cid:32) − √ , (cid:33) T , and the interlayer hopping matrices V = w (1 + σ x ) , V = w (cid:32) − σ x − √ σ y (cid:33) , V = w (cid:32) − σ x + √ σ y (cid:33) , where 1 stands for the 2 × v F ≈ · nm, w = 110meV, and k θ = | q j | = (8 π/ a ) sin( θ/ a = 0 . θ = 2 . ◦ .The reciprocal vectors of the TBG continuum model are given by g = q − q and g = q − q . By Fouriertransforming the zero field Hamiltonian (72) into the momentum space, the Hamiltonian becomes a model in ahoneycomb reciprocal lattice, where the orbitals of layer 1 and layer 2 are located at the two different sublattices Q ∈ q + g Z + g Z and Q ∈ − q + g Z + g Z of the honeycomb reciprocal lattice, respectively. Namely, theorigin of the momentum in layer ζ = 1 , − ζ − q , which is an example of choosing orbital-dependentmomentum origins in Eq. (29). Such a shift has the advantage of making the symmetries of the momentum spaceHamiltonian more explicit, and thus is adopted in most literatures of TBG. Since the kinetic energy for sublattices { Q } and { Q } in Eq. (72) are identical, we can use a single notation Q to denote both reciprocal sublattice sites { Q , Q } , i.e., the full honeycomb reciprocal lattice sites, and rewrite the Hamiltonian as H Q (cid:48) Q ( k ) = v F σ ∗ · ( k + Q ) δ Q (cid:48) Q + (cid:88) j =1 ( V j δ Q (cid:48) , Q + q j + h.c. ) . (73)One only needs to remember that the two different sublattices of the reciprocal lattice correspond to layers 1 and 2,respectively. The Hamiltonian under magnetic field in our basis is then given by the substitution k → ˆ κ + k , withˆ κ = √ (cid:96) ( a + a † , − ia + ia † ). In the calculation of main text Fig. 2, we set a twist angle θ = 2 . ◦ , take cutoffs N Q = 37and N L = 60, and choose the central momentum k at the Γ point of the first TBG BZ (A different choice of k only affects the spectrum at extremely small magnetic fluxes | ϕ/ π | < /N L , in which regime our calculation reducesto the LL calculation for the k · p model expanded at momentum k ). The same spectrum with edge states presentand with them deleted by the edge projection criteria of Eq. (61) in a larger energy range is shown in Fig. 4(a)-(b).The edge states in Fig. 4(b) are deleted following the method described in Sec. II E, where we used edge width w = min { (cid:96) − , . √ Ω BZ } for the edge projector P κ b ,w , and the edge projection threshold P c = 0 .
5. More examples ofthe TBG Hofstadter spectra calculated with our method can be found in Ref. [12].As another example, we also calculated the spectrum of the TBG Hamiltonian with a Dirac mass term added (whichcan arise from hBN substrate alignment in TBG [24, 25]), i.e., a model Hamiltonian H (cid:48) Q (cid:48) Q ( m D , k ) = [ v F σ ∗ · ( k + Q ) + m D σ z ] δ Q (cid:48) Q + (cid:88) j =1 ( V j δ Q (cid:48) , Q + q j + h.c. ) . (74)The calculated Hofstadter butterfly with edge states present and deleted for m D = 200meV are shown in Fig. 4(c)-(d).We will comment more on this case of nonzero Dirac mass in Sec. III B.8 (a) (b) (c) (d) ϕ/2π ϕ/2π ϕ/2π ϕ/2π FIG. 4. Hofstadter butterfly and spectral flow of θ = 2 . ◦ TBG with N Q = 37 and N L = 60 in a larger energy interval thanthat of main text Fig. 2, where (a)-(b) are calculated for a massless Dirac kinetic term, while (c)-(d) are calculated for a Diracmass m D = 200meV. The horizontal axis linear coordinate is equal to ϕ/ π in the range [0 , − π/ϕ in therange [1 , ∞ ]. For convenience, we still label the horizontal axis by the values of ϕ . The edge states in (b) and (d) are deletedfollowing the method described in Sec. II E, where we used w = min { (cid:96) − , . √ Ω BZ } for the edge projector P κ b ,w , and the edgeprojection threshold P c = 0 . B. Spurious zero modes
In the main text, Fig. 2(a) (see also Fig. 4(a)) shows many nondispersive horizontal levels at large fluxes ϕ , whichare roughly distributed in energies from − . . m D = 200meV,these nondispersive horizontal levels are still present at large ϕ , but are distributed in an energy range − . . ϕ is small, these levels become dispersive with respect to ϕ and tend to higher energies, merging eitherwith the Hofstadter bulk bands or with the dispersive momentum edge states. These nondispersive horizontal levelsat large fluxes are understood as spurious zero modes of Dirac fermions due to the LL cutoff N L , which can be derivedas follows.In the large ϕ limit, the magnetic length (cid:96) →
0, and the massless Dirac kinetic term of the TBG Hamiltonian (73)at momentum site Q will tend to v F σ ∗ · (ˆ κ + k + Q ) = √ v F (cid:96) (cid:18) a − p + (cid:96)/ √ a † − p − (cid:96)/ √ (cid:19) → √ v F (cid:96) (cid:18) aa † (cid:19) , (75)where p ± = k x + Q x ± i ( k y + Q y ). Therefore, up to error O ( p ± (cid:96) ), the Dirac kinetic term on each site Q under LLcutoff N L has two zero energy modes: ψ = (cid:18) | (cid:105) (cid:19) , ψ s = (cid:18) | N L (cid:105) (cid:19) , (76)where we have used the fact that a | (cid:105) = 0 and a † | N L (cid:105) = 0. Note that a † | N L (cid:105) = 0 is only true because of the LLcutoff. The first mode ψ is the physical Dirac zero mode, while the second mode ψ s is an unphysical spurious zeromode due to the LL cutoff N L . With a reciprocal lattice cutoff of N Q BZs, we have 2 N Q spurious zero modes on the2 N Q honeycomb reciprocal lattice (since there are two Dirac cones in each BZ). These spurious zero modes are locatedat radius | ˆ κ | = (cid:112) (cid:104) ˆ κ (cid:105) ≈ √ N L (cid:96) in the momentum space (main text Fig. 1(b)). For large fluxes ϕ > πN Q /N L , themomentum boundary is at radius κ b ≈ min { √ N L (cid:96) , (cid:113) N Q Ω BZ π } = (cid:113) N Q Ω BZ π < √ N L (cid:96) , thus these spurious modes atradius | ˆ κ | ≈ √ N L (cid:96) are outside the momentum space boundary and are not physical states of the Hofstadter butterfly.9Moreover, these spurious modes on different reciprocal sites Q have zero Dirac kinetic energy, and only hop amongnearest reciprocal sites with an amplitude ≈ ψ † s V j ψ s = w = 110 meV. Therefore, at large flux ϕ , they behave asa honeycomb tight-binding model in the reciprocal lattice with hopping amplitude w , which has 2 N Q energy levelsindependent of ϕ distributed between − w and 3 w . These levels give the horizontal lines at large ϕ in the maintext Fig. 2(a) (as well as Fig. 4(a)).When a Dirac mass term m D σ z is added to the TBG Hamiltonian as given in Eq. (74), one finds the spuriouszero mode ψ s in Eq. (76) at each momentum site Q no longer has a Dirac zero kinetic energy; instead, it is shiftedto energy m D . Therefore, with the hopping w among nearest momentum sites, one expects the 2 N Q spurious zeromodes to be distributed in the energy range between m D − w and m D + 3 w . This is exactly the case in Fig. 4(c).At small ϕ , one can no longer ignore the p ± terms in Eq. (75), thus the spurious modes ψ s in Eq. 76 is no longer azero energy mode of the Dirac Hamiltonian in Eq. (75). Therefore, one expect these spurious modes to disperse withrespect to ϕ at small ϕ , in agreement with the numerical result.When counting the number of occupied states in a fixed λ sector (in Eqs. (9) and (10) of the main text), theseunphysical spurious modes at large magnetic fields should be excluded. Further, since these spurious modes areoutside the momentum space boundary (when ϕ/ π > N Q /N L ), they will be identified as edge states by the boundaryprojector P κ b ,w in Eq. (61) (their expectation values (cid:104) P κ b ,w (cid:105) are close to 1), and thus will be removed in the edge-stateremoved spectrum (main text Fig. 2(b), and Fig. 4(b) and (d)). IV. BASIS COMPLETENESS AND MATRIX ELEMENTS FOR TIGHT-BINDING MODELS
The open momentum space method can also be applied to the numerical calculation of the Hofstadter butterflyof tight-binding models, as we demonstrated in the main text Fig. 3. To understand why the method is also validfor tight-binding models, in this section, we construct the complete and orthonormal basis for tight binding modelsemployed by our method. We prove that under the basis we construct, the tight-binding Hamiltonian under magneticfield (Peierls substitution) is block diagonalized into identical blocks, and each block has matrix elements given bythe zero-magnetic-field momentum space Hamiltonian H αβ ( k ) with the simple substitution k → ˆ κ + k , where k isthe quasi-momentum, k is an arbitrary momentum vector, ˆ κ = √ (cid:96) ( a + a † , − ia + ia † ), and a and a † are the Landaulevel raising and lowering operators. A. Standard Peierls substitution
We denote the Wannier orbital α in the unit cell labeled by lattice vector D ∈ d Z + d Z in real space as | D , α (cid:105) ,and use u α to denote the position of orbital α in a unit cell. In the continuum space, we have (cid:104) r , α | D , α (cid:105) = W α ( r − D − u α ) , (77)where W α ( r ) is the Wannier function of orbital α , and | r , α (cid:105) is the underlying continuum space basis at position r ofintrinsic orbital α as defined by Eq. (25). For the discussion of tight-binding models here, without loss of generality,we shall choose the gauge where all p α = 0 in the definition of | r , α (cid:105) in Eq. (25) (recall that p α is the momentumspace origin of orbital α , which can be chosen freely), namely, | r , α (cid:105) = c † α ( r ) | (cid:105) . (78)This ensures that the Hamiltonian matrix in the real space basis is invariant under lattice vector translation T d i (instead of changing by a unitary transformation as shown in Eq. (30) when p α (cid:54) = 0).We first comment that the derivation of the standard Peierls substitution [21] requires an approximation that theWannier orbitals | D , α (cid:105) are infinitely localized, namely, W α ( r ) = δ ( r ). This is because the Peierls substitution onlypicks up the gauge phase factor connecting two points of Wannier positions, and is independent of the details (shapes,sizes, etc) of Wannier functions, which can be true only if each Wannier function is infinitely small and thus does notfeel the magnetic field inside the orbital itself. This does not require the Wannier charge density to be localized onthe site position D , instead it can be a delta function localized at any position D + u α away from the site position D .However, in our paper here, we shall keep the Wannier function W α ( r ) in Eq. (77) generic , instead of assumingit is a delta function. This makes our proof of the method the most generic, which applies to nonstandard Peierlssubstitutions discussed in Sec. IV B, too.0The tight-binding model in real space then generically takes the form H = (cid:88) j,α,β t αβj T αβ D j + u α − u β , (79)where t αβj are the hopping amplitudes, and we have defined T αβ D j + u α − u β = (cid:88) D e i (cid:82) cj,αβ A ( D + r ) · d r | D + D j , α (cid:105)(cid:104) D , β | (80)as the translation operator from the Wannier orbital | D , β (cid:105) at position D + u β to the Wannier orbital | D + D j , α (cid:105) at position D + D j + u α under the Peierls substitution of the gauge potential A ( r ) in the continuum space, with D j ∈ d Z + d Z , and c j,αβ being the straight line segment from u β to D j + u α .—– At zero magnetic field , we can choose the gauge A ( r ) = 0, and the Hamiltonian can be written into themomentum space by Fourier transformation as H αβ ( k ) = (cid:88) j t αβj e − i k · ( D j + u α − u β ) , (81)where k is the quasi-momentum which takes values in the Brillouin zone, and the basis is the Bloch basis | k , α (cid:105) = N Ω (cid:80) D e i k · ( D + u α ) | D , α (cid:105) , with N Ω being the number of unit cells in real space.—– At nonzero magnetic field , the gauge potential A ( r ) satisfies ∂ x A y ( r ) − ∂ y A x ( r ) = B , where B is a uniformmagnetic field in the continuum space. We now proceed to define a complete orthonormal basis for Hamiltonian (79)based on the continuum space.Before starting, we first note that the continuum space has a Hilbert space spanned by the real space basis | r , α (cid:105) of all positions r , which is much larger than the Hilbert space of the tight binding model spanned by the Wannierbasis | D , α (cid:105) defined in Eq. (77). In the following, we shall first define a complete basis for the continuum space; thenwe project the basis into the sub-Hilbert space of the tight-binding model spanned by | D , α (cid:105) using a projector, andprove that the resulting projected basis forms a complete orthonormal basis for the tight-binding model.In the continuum space, we can define the kinematic momentum operator Π = − i ∇ − A ( r ), and the guiding centeroperator R = r − (cid:96) (cid:126) ˆ z × Π , as we did in Sec. I. Similar to Sec. I, we define a guiding center R ˆ τ along the ˆ τ direction,with ˆ τ · (ˆ z × g )ˆ τ · (ˆ z × g ) irrational. Moreover, we denote the reciprocal lattice of the tight-binding model as Q ∈ g Z + g Z . Wecan then define lowering and raising operators following the same procedure as we did in Sec. II a Q = (cid:96) √ x − Q x − k ,x + i (Π y − Q y − k ,y )] , a † Q = (cid:96) √ x − Q x − k ,x − i (Π y − Q y − k ,y )] , (82)where k is the center momentum which can be chosen freely. Afterwards, we can define a basis | λ, Q , n, α (cid:105) in thecontinuum space satisfying R ˆ τ | λ, Q , n, α (cid:105) = [ λ − (cid:96) ˆ τ · (ˆ z × Q )] | λ, Q , n, α (cid:105) , a † Q a Q | λ, Q , n, α (cid:105) = n | λ, Q , n, α (cid:105) , (83)where λ takes values in the quotient set defined in Eq. (43), namely, λ ∈ Λ ˆ τ = R / [ (cid:96) ˆ τ · (ˆ z × g ) Z + (cid:96) ˆ τ · (ˆ z × g ) Z ].The basis definition (83) follows exactly the same derivation of the basis | λ, Q α , n, α (cid:105) in Eq. (40) in Sec. II, exceptthat here we have chosen the gauge that the momentum origins of all orbitals α are at p α = (see Eq. (25) for thedefinition of p α , and see Eq. (78) for our gauge choice for tight binding models here), thus Q α = Q for all α (recallthe definition of Q α in Eq. (II)). As we have proved in Eq. (44), the basis | λ, Q , n, α (cid:105) in Eq. (83) satisfies (cid:104) λ (cid:48) , Q (cid:48) , n (cid:48) , β | λ, Q , n, α (cid:105) = δ βα δ λλ (cid:48) δ Q (cid:48) Q δ n (cid:48) n , (84)thus forms a complete orthonormal basis for the Hilbert space of the continuum space, where λ ∈ Λ ˆ τ = R / [ (cid:96) ˆ τ · (ˆ z × g ) Z + (cid:96) ˆ τ · (ˆ z × g ) Z ] as defined in Eq. (43), Q ∈ g Z + g Z are the reciprocal vectors, and n ≥ N Q (cid:88) λ, Q ,n | λ, Q , n, α (cid:105)(cid:104) λ, Q , n, α | = α = (cid:90) d r | r , α (cid:105)(cid:104) r , α | , (85)where N Q is the number of reciprocal sites Q (which tends to infinity), and α stands for the identity matrix in theorbital α subspace in the continuum space.1Next, based on the continuum space basis | λ, Q , n, α (cid:105) , we would like to define a complete basis for the Hilbert spaceof the tight-binding model spanned by Wannier orbitals | D , α (cid:105) . We define the basis for the tight-binding model as | λ, n, α (cid:105) = (cid:88) D | D , α (cid:105)(cid:104) D + u α , α | λ, , n, α (cid:105) , (86)where | D , α (cid:105) is the Wannier orbital, | D + u α , α (cid:105) denotes the position eigenbasis at position r = D + u α (the centerof the Wannier orbital), and λ ∈ Λ ˆ τ = R / [ (cid:96) ˆ τ · (ˆ z × g ) Z + (cid:96) ˆ τ · (ˆ z × g ) Z ] as defined in Eq. (43). One may wonderwhy only the continuum space states | λ, , n, α (cid:105) at Q = are used in defining the basis (86). In fact, by Eq. (48), wehave the following identity for any lattice vector D and reciprocal vector Q (in the equation below, r is understoodas a number denoting the continuous space coordinates) (cid:104) D + u α , α | λ, Q , n, α (cid:105) = (cid:90) d r δ ( r − D − u α ) (cid:104) r , α | e i Q · r | λ, , n, α (cid:105) = e i Q · ( D + u α ) (cid:104) D + u α , α | λ, , n, α (cid:105) = e i Q · u α (cid:104) D + u α , α | λ, , n, α (cid:105) . (87)Thus, one can equivalently rewrite the basis definition as | λ, n, α (cid:105) = 1 (cid:112) N Q (cid:88) Q , D e − i Q · u α | D , α (cid:105)(cid:104) D + u α , α | λ, Q , n, α (cid:105) , (88)which is expressed using the continuum space basis of all reciprocal sites Q .By the orthonormal relation (84) and the basis expression in Eq. (86), it is easy to see the subset of basis | λ, n, α (cid:105) is orthonormal: (cid:104) λ (cid:48) , n (cid:48) , β | λ, n, α (cid:105) = δ βα (cid:88) D (cid:104) λ (cid:48) , , n (cid:48) , α | D + u α , α (cid:105)(cid:104) D + u α , α | λ, , n, α (cid:105) = δ βα (cid:104) λ (cid:48) , , n (cid:48) , α | (cid:88) Q e i Q · ( r − u α ) | λ, , n, α (cid:105) = δ βα (cid:88) Q e − i Q · u α (cid:104) λ (cid:48) , , n (cid:48) , α | λ, Q , n, α (cid:105) = δ λλ (cid:48) δ n (cid:48) n δ βα , (89)where in the 2nd line r is understood as the position operator (not number), and we have used Eq. (48) which implies | λ, Q , n, α (cid:105) = e i Q · r | λ, , n, α (cid:105) . Besides, from the 1st line to the 2nd line we have used the following identity for eachorbital α : (cid:88) D | D + u α , α (cid:105)(cid:104) D + u α , α | = (cid:90) d r | r , α (cid:105)(cid:104) r , α | (cid:88) D δ ( r − D − u α )= (cid:90) d r | r , α (cid:105)(cid:104) r , α | (cid:88) Q e i Q · ( r − u α ) = (cid:88) Q e i Q · ( r − u α ) , (90)where r on the right-most hand side of the 2nd line is understood as the position operator (instead of a number).Furthermore, by Eq. (87), we can prove the completeness of the basis (86) for the tight-binding model as follows: (cid:88) λ,n,α | λ, n, α (cid:105)(cid:104) λ, n, α | = (cid:88) λ,n,α (cid:88) D , D (cid:48) | D , α (cid:105)(cid:104) D + u α , α | λ, , n, α (cid:105)(cid:104) λ, , n, α | D (cid:48) + u α , α (cid:105)(cid:104) D (cid:48) , α | = 1 N Q (cid:88) λ, Q ,n,α (cid:88) D , D (cid:48) | D , α (cid:105) e − i Q · u α (cid:104) D + u α , α | λ, Q , n, α (cid:105)(cid:104) λ, Q , n, α | D (cid:48) + u α , α (cid:105) e i Q · u α (cid:104) D (cid:48) , α | = (cid:88) α (cid:88) D , D (cid:48) | D , α (cid:105)(cid:104) D + u α , α | N Q (cid:88) λ, Q ,n | λ, Q , n, α (cid:105)(cid:104) λ, Q , n, α | | D (cid:48) + u α , α (cid:105)(cid:104) D (cid:48) , α | = (cid:88) D , D (cid:48) ,α | D , α (cid:105)(cid:104) D + u α , α | (cid:18)(cid:90) d r | r , α (cid:105)(cid:104) r , α | (cid:19) | D (cid:48) + u α , α (cid:105)(cid:104) D (cid:48) , α | = (cid:88) D ,α | D , α (cid:105)(cid:104) D , α | . (91)This proves that the basis | λ, n, α (cid:105) forms a complete orthonormal basis for the Hilbert space of the tight-bindingmodel spanned by Wannier orbitals | D , α (cid:105) .Now we discuss the translation operator in Eq. (80) under magnetic field. First, we prove the following identity ofthe position basis | r , α (cid:105) in the continuum space. Assume c f is a path (not necessarily straight) from position r to2position r f in the continuum space. By path partitioning c f into N small segments δ r j = r j − r j − (1 ≤ j ≤ N ),and note that e − i Π · δ r = e i A ( r ) · δ r e − δ r ·∇ in the δ r → P e − i (cid:82) cf Π · d r | r , α (cid:105) = lim N →∞ N (cid:89) j =1 e − i Π · δ r j | r , α (cid:105) = lim N →∞ (cid:90) d r N | r N , α (cid:105)(cid:104) r N , α | e − i Π · δ r N (cid:90) d r N − | r N − , α (cid:105)(cid:104) r N − , α | e − i Π · δ r N − · · · (cid:90) d r | r , α (cid:105)(cid:104) r , α | e − i Π · δ r | r , α (cid:105) = lim N →∞ N (cid:89) j =1 e i A ( r + (cid:80) i ≤ j δ r i ) · δ r j | r + N (cid:88) j =1 δ r j , α (cid:105) = e i (cid:82) cf A ( r ) · d r | r f , α (cid:105) , (92)where P stands for path ordering, and we have used the fact that e − δ r ·∇ | r , α (cid:105) = | r + δ r , α (cid:105) . Therefore, note that thedefinition of the basis | λ, n, α (cid:105) in Eq. (86) contains the position eigenstate | D + u α , α (cid:105) at position r = D + u α , wecan rewrite the action of the translation operator in Eq. (80) on the basis | λ, n, α (cid:105) as T αβ D j + u α − u β | λ, n, β (cid:48) (cid:105) = δ ββ (cid:48) (cid:88) D e i (cid:82) cαβ A ( D + r ) · d r | D + D j , α (cid:105)(cid:104) D + u β , β | λ, , n, β (cid:105) = δ ββ (cid:48) (cid:88) D | D + D j , α (cid:105)(cid:104) D + D j + u α , β |P e − i (cid:82) cj,αβ Π · d r | λ, , n, β (cid:105) , (93)where notations of the form | D + u α , β (cid:105) denotes the position basis of orbital β in the continuum space at position D + u α (note that β need not be equal to α ). In contrast, | D + D j , α (cid:105) denotes the Wannier orbital basis of the unitcell at D + D j , as we have defined. We note that in Eq. (93), from the 1st line to the 2nd line we have used thefact that the position basis in the continuum space satisfies (cid:104) D + u β , β | = (cid:104) D + D j + u α , β |P e − i (cid:82) cj,αβ Π · d r , which wehave proved in Eq. (92).We then define lowering and raising operators a = a = (cid:96) √ [Π x − k ,x + i (Π y − k ,y )] and a = a † , where a , a † arethe lowering and raising operators in Eq. (82) at reciprocal site Q = . By further defining ˆ κ = √ (cid:96) ( a + a † , − ia + ia † ),we can rewrite the kinematic momentum as Π = ˆ κ + k , where k is the center momentum in Eq. (82). We can thenfurther simplify Eq. (93) as T αβ D j + u α − u β | λ, n, β (cid:48) (cid:105) = δ ββ (cid:48) (cid:88) D | D + D j , α (cid:105)(cid:104) D + D j + u α , β |P e − i (cid:82) cj,αβ (ˆ κ + k ) · d r | λ, , n, β (cid:105) , = δ ββ (cid:48) (cid:88) D | D + D j , α (cid:105)(cid:104) D + D j + u α , α |P e − i (cid:82) cj,αβ (ˆ κ + k ) · d r | λ, , n, α (cid:105) , = δ ββ (cid:48) (cid:88) D | D , α (cid:105)(cid:104) D + u α , α |P e − i (cid:82) cj,αβ (ˆ κ + k ) · d r | λ, , n, α (cid:105) = δ ββ (cid:48) P e − i (cid:82) cj,αβ (ˆ κ + k ) · d r | λ, n, α (cid:105) , (94)where P stands for path ordering. From the 1st line to the 2nd line of Eq. (94), we have used the fact that the matrixelements of the operator P e − i (cid:82) cj,αβ (ˆ κ + k ) · d r (acting from the right on the state (cid:104) D + D j + u α , β | ) only depend onthe position r = D + D j + u α of the state (since ˆ κ is defined in terms of Π which displaces r ), and are independentof the orbital index β . So we can simply change the orbital index β in the 1st line into α in the 2nd line. In thederivation from the 3rd line to the 4-th line, we have also used the fact that the matrix elements of operator ˆ κ onlydepends on the LL quantum number n , so in the last line of Eq. (94) one should understand the operator ˆ κ as solelyacting on the quantum number n of basis | λ, n, α (cid:105) . More concretely, the LL lowering and raising operators in ˆ κ actas a | λ, n, α (cid:105) = √ n | λ, n − , α (cid:105) and a † | λ, n, α (cid:105) = √ n + 1 | λ, n + 1 , α (cid:105) .For standard Peierls substitution, the path c j,αβ in Eq. (94) is a straight line segment (see the paragraph belowEq. (80)), so the integral on the exponent of P e − i (cid:82) cj,αβ (ˆ κ + k ) · d r is along the straight line segment from D + u β to D + D j + u α . Since the matrix elements of ˆ κ only depend on n when acting on | λ, n, α (cid:105) , the matrix ˆ κ + k shouldbe a constant along the path of integration, and thus we can further simplify Eq. (94) into T αβ D j + u α − u β | λ, n, β (cid:48) (cid:105) = δ ββ (cid:48) e − i (ˆ κ + k ) · ( D j + u α − u β ) | λ, n, α (cid:105) . (95)3 E ϕ /2 π FIG. 5. The first 4 periods of Hofstadter butterfly of the tight-binding model H ( k ) = − cos k x − cos k y in the square latticecalculated using our method of substitution k → ˆ κ + k , where we take a LL cutoff N L = 100. The cutoff N L breaks theperiodicity of the Hofstadter spectrum, and the larger the magnetic field B is, the less clear the Hofstadter butterfly is. Forsufficiently large LL cutoffs (e.g., N L ∼ It is evident from Eq. (95) that the matrix elements of T D j + u α − u β is diagonal in quantum number λ and independentof λ . Therefore, we can rewrite the translation operator in a fixed λ subspace from basis | λ, n, β (cid:105) to basis | λ, n (cid:48) , α (cid:105) ( n, n (cid:48) ∈ Z ) are nonnegative integers) as T λ,αβ D j + u α − u β = e − i (ˆ κ + k ) · ( D j + u α − u β ) , (96)where ˆ κ acting on the LL quantum number n (recall that Π = ˆ κ + k ). Accordingly, the tight-binding Hamiltonianin magnetic field in a fixed λ subspace from basis | λ, n, β (cid:105) to basis | λ, n (cid:48) , α (cid:105) ( n, n (cid:48) ∈ Z ) takes the form H λ,αβ = (cid:88) j t αβj T λ,αβ D j + u α − u β = (cid:88) j t αβj e − i (ˆ κ + k ) · ( D j + u α − u β ) , (97)which is exactly the zero field momentum space Hamiltonian (81) with the substitution k → ˆ κ + k .As an example of the above method, we numerically calculate the Hofstadter butterfly of the simplest squarelattice model with one orbital per site and nearest hopping amplitude − /
2, which has a momentum space Hamil-tonian H ( k ) = − cos k x − cos k y . The Hamiltonian matrix with LL cutoff N L is constructed by replacing k x bythe Hermitian matrix k ,x N L + ( a + a † ) / ( √ (cid:96) ) and similarly for k y , where a is a matrix with matrix elements[ a ] mn = (cid:104) m | a | n (cid:105) = √ nδ m,n − ( m, n are integers from 0 to N L ). The cosine of a Hermitian matrix W can be calcu-lated by first diagonalizing the matrix into W = U † D W U where D W is diagonal and U is unitary. Then cos W canbe calculated efficiently using the identity cos W = U † (cos D W ) U and the fact that cos D W is simply the cosine ofeach element of the diagonal matrix D W .The spectra with cutoffs N L = 100 and N L = 500 are shown in the main text Fig. 3. Fig. 3 shows that thehigher cutoff N L is, the clearer the Hofstadter butterfly is. Furthermore, the cutoff N L breaks the periodicity of theHofstadter spectrum with respect to magnetic flux per unit cell ϕ . Fig. 5 shows the first 4 periods of the numericalHofstadter butterfly with N L = 100. In particular, the larger ϕ is, the less clear the Hofstadter butterfly is. This isbecause in Eq. (97) the operator ˆ κ ∝ (cid:96) − ( a, a † ); for a fixed cutoff N L , the error in the LL operators ( a, a † ) is fixed,so the error in operator ˆ κ is larger for larger ϕ (which corresponds to smaller (cid:96) ).Lastly, we note that throughout our derivation, we have not used the details of the Wannier function W α ( r ) in Eq.(77). Therefore, our derivation holds generically as long as the Peierls substitution is valid, independent of the shapeof the Wannier orbitals. B. Nonstandard Peierls substitution
In some models the Wannier orbitals cannot be approximated as localized at one position (but dominantly localizesat several positions respecting the symmetries of the lattice), so the standard Peierls approximations are no longervalid. However, in certain models, nonstandard Peierls substitution can be derived under certain approximations4[12]. Accordingly, the Peierls substitution might not be along the straight line segment path c αβ , but may be alonga nonstraight path or the sum of multiple nonstraight paths from D + u β to D + D j + u α . For instance, in the4-band tight-binding model for TBG in Ref. [26], the Wannier orbital is centered at AB and BA stackings of TBG butextends to the 3 closest AA stacking centers, and it is shown to have an approximate nonstandard Peierls substitutiongiven by the summation of contributions of multiple broken-line paths [12].For such a tight-binding model with nonstandard Peierls substitution, the tight-binding Hamiltonian is still of theform of Eq. (79), except that the translation operator in Eq. (79) now reads T αβ D j + u α − u β = (cid:88) D ,µ e i (cid:82) cµj,αβ A ( D + r ) · d r | D + D j , α (cid:105)(cid:104) D , β | , (98)where the phase factor involves the summation of multiple different paths (e.g., broken-lines) c µj,αβ (labeled by index µ = 1 , , · · · ) from position u β to D j + u α . For instance, in the 4-band tight-binding model for TBG studied in[12], the Peierls substitution between the nearest neighbors is given by the summation of the gauge phase factor of 2different paths (i.e., µ = 1 ,
2) from one AB site to another AB site via the nearest 2 AA sites.Our method can still apply to such tight-binding models with nonstandard Peierls substitutions. To see this, wefirst note that throughout our derivations in Sec. IV A for the standard Peierls substitution case, we do not requireat all the Wannier orbitals in Eq. (77) to be localized. Therefore, in the nonstandard Peierls substitution case here,we can still define the eigenbasis by Eq. (86), which still satisfies the orthogonality in Eq. (89) and the completenessin Eq. (91).The action of the translation operator is still given by Eq. (94), except that one need to sum over all the paths c µj,αβ (i.e., replace c j,αβ in Eq. (94) by c µj,αβ and sum over the path index µ ). However, since c µj,αβ are no longer straightpaths, Eq. (94) cannot be further reduced to Eq. (95). Therefore, the translation operator in a fixed λ sector frombasis | λ, n, β (cid:105) to | λ, n (cid:48) , α (cid:105) can be expressed by Eq. (94) as T λ,αβ D j + u α − u β = (cid:88) µ P e − i (cid:82) cµαβ (ˆ κ + k ) · d r , (99)where ˆ κ = √ (cid:96) ( a + a † , − ia + ia † ) acts on the LL quantum number n , and P stands for path ordering. The Hamiltonianin a fixed λ subspace then reads H λ,αβ = (cid:80) j t αβj T λ,αβ D j + u α − u β , the matrix elements of which are independent of λ . V. REVIEW OF THE DIOPHANTINE EQUATION
In this section, we briefly review the proof of the Diophantine equation.For a lattice model with magnetic flux ϕ = 2 π pq per unit cell, where p and q are two coprime numbers, the energyspectrum forms a set of Hofstadter bands. In particular, each Hofstadter gap is characterized by two integers t ν and s ν , which satisfy the Diophantine equation t ν p + s ν q = ν , (100)where ν is an integer. In particular, t ν is the Chern number of the Hofstadter gap.Here we briefly review how the Diophantine equation is proved following Ref. [15] (See also Ref. [27] Chapter 5.3 fora different proof). In fact, one can prove an equivalent statement, that each Hofstadter band satisfies a Diophantineequation σp + mq = 1 , (101)where σ is the Chern number of the Hofstadter band, and m is another integer characterizing the band.The proof is as follows. First, we note that a lattice Hamiltonian H (either a continuum model in Eq. (27)or a tight-binding model as in Eq. (79) we considered) in a uniform magnetic field B = B ˆ z still has translationsymmetries along the Bravais lattice vectors d and d . However, the translation symmetry operators are not thesimple translation operators T d j = e − i Π · d j ( j = 1 ,
2) (102)in magnetic field in the continuum space, since the operator T d j do not commute with Hamiltonian H as one caneasily verify (since H contains operator Π , and [Π x , Π y ] (cid:54) = 0). Instead, the (magnetic) translation symmetry operatorswhich commute with the Hamiltonian H are given by (for simplicity here we choose all the p α = in Eq. (29)): (cid:101) T d j = T d j e − i ( (cid:96) − ˆ z × r ) · d j = e − i ( Π + (cid:96) − ˆ z × r ) · d j = e − i(cid:96) − (ˆ z × R ) · d j , [ (cid:101) T d j , H ] = 0 , ( j = 1 ,
2) (103)5 g /q g | ψ > n, k | ψ > n, k + g =e | ψ > n, k i σ q d k | ψ > n, k + g /q =| ψ > n, k FIG. 6. Illustration of the gauge choices in Eq. (106) for proving the Diophantine equation, where the magnetic flux is ϕ/ π = p/q . where R is the guiding center operator. Namely, the translation operator T d j and the translation symmetry operator (cid:101) T d j differ by a unitary transformation e − i ( (cid:96) − ˆ z × r ) · d j . At zero magnetic field, T d j and (cid:101) T d j become the same.At rational flux ϕ = 2 π pq per unit cell, it is straightforward to show that the translation symmetry operators satisfythe commutation relation (cid:101) T d (cid:101) T d = e − i πp/q (cid:101) T d (cid:101) T d . (104)Note that in contrast we have T d T d = e i πp/q T d T d . One can therefore define two magnetic translation symmetryoperators (cid:101) T q d = ( (cid:101) T d ) q and (cid:101) T d which commute with each other, i.e., [ (cid:101) T q d , (cid:101) T d ] = 0. Since they also commutewith the Hamiltonian H , we can define the Bloch wave function eigenstates | ψ n, k (cid:105) of a Hofstadter band n of H , withquasimomentum k defined by (cid:101) T q d | ψ n, k (cid:105) = e iq d · k | ψ n, k (cid:105) , (cid:101) T d | ψ n, k (cid:105) = e i d · k | ψ n, k (cid:105) . (105)The magnetic BZ is then a parallelogram spanned by momentum vectors g /q and g , where g i satisfies g i · d j = 2 πδ ij ( i, j = 1 , k in the magnetic BZ, the energy eigenvalues are nondegenerate), which is generically true whenthere is no other symmetries (which may protect degeneracies) except for the translation symmetries. When theHofstadter band n has a Chern number σ , one can choose the Bloch wave function | ψ n, k (cid:105) as a continuous function of k satisfying | ψ n, k + g /q (cid:105) = | ψ n, k (cid:105) , | ψ n, k + g (cid:105) = e iσq d · k | ψ n, k (cid:105) , (106)as illustrated in Fig. 6. Here we do not restrict k within the magnetic BZ. One can easily verify that the above choicegives a Berry phase σq d · ( g /q ) = 2 πσ circulating the boundary of the magnetic BZ (starting from momentum ( , )to ( , g ) to ( g /q, g ) to ( g /q, ) and then back to ( , )), which is required by the Chern number σ .By Eq. (104) and the definition of quasimomentum k in Eq. (105), the state (cid:101) T d | ψ n, k (cid:105) should be equal to theBloch state | ψ n, k + p g /q (cid:105) up to a phase factor, and should be degenerate with the state | ψ n, k (cid:105) . In consistency withEq. (106), we can in general choose the phase such that (cid:101) T d | ψ n, k (cid:105) = e imq d · k | ψ n, k + p g /q (cid:105) , (107)where m is some coefficient to be determined. First, one can prove that m ∈ Z is an integer: this is because the firstequation in (106) requires e imq d · k | ψ n, k + p g /q (cid:105) = (cid:101) T d | ψ n, k (cid:105) = (cid:101) T d | ψ n, k + g /q (cid:105) = e imq d · ( k + g /q ) | ψ n, k + g /q + p g /q (cid:105) = e imq d · ( k + g /q ) | ψ n, k + p g /q (cid:105) , namely, e imq d · ( g /q ) = e i πm = 1. We can then apply (cid:101) T d by q times to obtain (cid:101) T q d | ψ n, k (cid:105) = e imq d · k | ψ n, k + p g (cid:105) = e i ( σp + mq ) q d · k | ψ n, k (cid:105) = e iq d · k | ψ n, k (cid:105) . (108)Since k can take any value, we conclude that the Diophantine equation 101 for the band has to hold.It is then straightforward to see the Diophantine equation 100 holds in a Hofstadter gap: ν is simply the numberof Hofstadter bands below the gap counted from some reference point of filling (see Sec. II G for discussion of thereference point of filling). The Hofstadter band between the ν -th gap and the ( ν − σ = t ν − t ν − , and the other quantum number m = s ν − s ν − .6Since ν has the physical meaning of number of occupied magnetic bands, and the magnetic unit cell has an area q times of the original unit cell, we can define the number of electrons per original unit cell as ρ = ν/q . One can thendivide the Diophantine equation by q to obtain our main text Eq. (7): t ν ϕ π + s ν = ρ , (109)which holds even when the flux quanta ϕ/ π is irrational. The derivative of this equation with respect to ϕ gives theStreda formula 2 π dn r dϕ = t r [20]. VI. QUANTIZED LORENTZ SUSCEPTIBILITY FROM s ν In this section, we show that the quantum number s ν , which can be viewed as a dual Chern number in themomentum space, yields a quantized Lorentz susceptibility γ xy = eBs ν . This is analogous to the Hall conductancegiven by the Chern number t ν . A. Review of the quantized Hall conductance
We first briefly review how the Chern number t ν leads to a quantized Hall conductance σ xy = t ν e h [10]. We assumethe flux per unit cell is ϕ = 2 πp/q , so that we can define the magnetic BZ quasimomentum ( k x , k y ) by Eq. (105).Using the Kubo formula, the Hall conductance under a uniform electric field E y can be written as σ xy = i ∂∂ω (cid:90) ∞−∞ d ω (cid:48) (cid:104) G ω + ω (cid:48) ˆ j x G ω (cid:48) ˆ j y (cid:105) (cid:12)(cid:12)(cid:12) ω → (110)where ˆ j is the spatially uniform current operator which is conjugate to a spatially uniform gauge field A ., while G ω = ( ω − H ) − is the Green’s function in the magnetic field B at energy ω , and H is the (first quantized) fullsingle-particle Hamiltonian H under magnetic field in the real space (under the position basis | r , α (cid:105) ). For example,for continuum models H is as defined in Eq. (34), while for tight-binding models H is as given in Eq. (79) embeddedinto the continuum space by the infinitely localized Wannier orbital assumption (77). The coefficient σ xy correspondsto a Hall response j x = − σ xy dA y /dt = σ xy E y , where E y = − dA y /dt is the uniform electric field in the y direction.We denote the Bloch eigenstates of H in the n -th band as | ψ n, k (cid:105) satisfying H | ψ n, k (cid:105) = (cid:15) n, k | ψ n, k (cid:105) , (111)where the magnetic BZ quasimomentum k is defined by Eq. (105) using the magnetic translation symmetry operators (cid:101) T q d = e − iq(cid:96) − (ˆ z × R ) · d and (cid:101) T d = e − i(cid:96) − (ˆ z × R ) · d , and (cid:15) n, k is the energy of Bloch state | ψ n, k (cid:105) . Since H is diagonal inthe magnetic BZ quasimomentum k , the Green’s function (in the magnetic field B ) is also diagonal in k , and can bewritten as G ω = (cid:88) n, k | ψ n, k (cid:105)(cid:104) ψ n, k | ω − (cid:15) n, k . (112)We comment that under the continuum space basis | λ, Q α , n, α (cid:105) in Eq. (40) employed by our Hofstadter method, theBloch state | ψ n, k (cid:105) does not live in a definite λ sector. This is can be seen by noting that | ψ n, k (cid:105) is the eigenstate of (cid:101) T q d and (cid:101) T d , while we have (cid:101) T q d R ˆ τ (cid:101) T † q d = R ˆ τ + q ˆ τ · d = R ˆ τ − p(cid:96) ˆ τ · (ˆ z × g ) and (cid:101) T d R ˆ τ (cid:101) T † d = R ˆ τ + ˆ τ · d = R ˆ τ + pq (cid:96) ˆ τ · (ˆ z × g ).Since the quantum number λ ∈ Λ ˆ τ = R / [ (cid:96) ˆ τ · (ˆ z × g ) Z + (cid:96) ˆ τ · (ˆ z × g ) Z ] as defined in Eq. (43), we find that theaction of (cid:101) T q d preserves the coset of λ , while (cid:101) T d maps the sector of coset λ to the sector of coset λ + pq (cid:96) ˆ τ · (ˆ z × g ).Therefore, each Bloch state | ψ n, k (cid:105) must be the superposition of the degenerate eigenstates in multiple of q sectorslabeled by quantum numbers λ + m pq (cid:96) ˆ τ · (ˆ z × g ) ( λ belongs to a certain set, and m = 0 , · · · , q − | λ, Q α , n, α (cid:105) will be complicated, since the basis state | λ, Q α , n, α (cid:105) withdefinite quantum numbers λ, Q α , n is not translationally invariant.For any state preserving the translation symmetry, the expectation value of the current operator ˆ j satisfies (recallwe have set charge e = 1) (cid:104) ˆ j (cid:105) = (cid:104) d r d t (cid:105) = − i (cid:104) [ r , H ] (cid:105) . (113)7The uniform current operator ˆ j preserves the translation symmetry and is thus diagonal in magnetic BZ quasimo-mentum k . In contrast, the position operator r does not commute with translation symmetry operators (cid:101) T q d and (cid:101) T d ,thus r is not fully diagonal in magnetic BZ quasimomentum k , and accordingly [ r , H ] is not fully diagonal in k (seeRef. [28]). Therefore, Eq. (113) shows that the current operator ˆ j is given by the part of − i [ r , H ] that is diagonal inmagnetic BZ quasimomentum k , namely,ˆ j = (cid:88) m,n, k ˆ j m k ,n k | ψ m, k (cid:105)(cid:104) ψ n, k | , ˆ j m k ,n k = − i (cid:104) ψ m, k | [ r , H ] | ψ n, k (cid:105) . (114)For a given magnetic BZ quasimomentum k (within the magnetic BZ spanned by vector g /q and g ), we can definea reduced magnetic BZ momentum space Hamiltonian ˆ H ( k ) (we label it by a hat, since it is distinguished from thefull Hamiltonian H in field in Eq. (34) for continuum models or in Eq. (79) for tight-binding models):ˆ H ( k ) = P ( e − i k · r He i k · r ) P , P = qN Ω (cid:88) j ,j ∈ Z (cid:101) T qj d + j d , (115)where k is understood as a parameter, r is the position operator (understood as (cid:80) α (cid:82) r | r , α (cid:105)(cid:104) r , α | d r when expandedin position basis, with α being the orbital), P is the projection operator into the magnetic translationally invariantsubspace, (cid:101) T d = e − i(cid:96) − (ˆ z × R ) · d is the translation operator of displacement d by the guiding center operator R asdefined in Eq. (103), and N Ω is the number of zero-magnetic-field unit cells (there are N Ω /q magnetic unit cellsat magnetic flux per unit cell ϕ = 2 πp/q ). Note that any state with a nonzero momentum k (cid:54) = 0 in the magneticBZ, for which the translation operator (cid:101) T qj d + j d has an eigenvalue e i ( qj d + j d ) · k , will be annihilated by P (since (cid:80) j ,j e i ( qj d + j d ) · k = N Ω q δ k , ). Thus, unlike the full Hamiltonian H in the real space which has a large Hilbert spacedimension equal to the number of bands times the number of k in the magnetic BZ, the reduced Hamiltonian ˆ H ( k )is defined at a fixed k and has a smaller Hilbert dimension equal to the number of bands, which can be understoodas the full Hamiltonian H projected into the sub-Hilbert space at magnetic BZ momentum k . An explicit example ofˆ H ( k ) for a tight-binding model is given in Eqs. (158) and (159), where one can see ˆ H ( k ) is simply the magnetic BZmomentum space Hamiltonian under the Fourier transformed basis of the Wannier orbitals. The Hamiltonian ˆ H ( k )has eigenstates | u n, k (cid:105) = (cid:88) α (cid:90) d r e − i k · r | r , α (cid:105)(cid:104) r , α | ψ n, k (cid:105) , ˆ H ( k ) | u n, k (cid:105) = (cid:15) n, k | u n, k (cid:105) , (116)where n runs over all the bands in the magnetic BZ. Note that the state | u n, k (cid:105) has is magnetically translationallyinvariant, namely, (cid:101) T q d | u n, k (cid:105) = (cid:101) T d | u n, k (cid:105) = | u n, k (cid:105) . Therefore, P | u n, k (cid:105) = | u n, k (cid:105) , and | u n, k (cid:105) can be understood asthe periodic part of the Bloch wave function | ψ n, k (cid:105) . In particular, from the definition of ˆ H ( k ) in Eq. (115) and thedefinition of current operator ˆ j in Eq. (114), we have (calculated by inserting the position basis) (cid:104) ψ m, k | ˆ j | ψ n, k (cid:105) = − i (cid:88) α,α (cid:48) (cid:90) d r d r (cid:48) (cid:104) ψ m, k | r , α (cid:105)(cid:104) r , α | ( r H − H r (cid:48) ) | r (cid:48) , α (cid:48) (cid:105)(cid:104) r (cid:48) , α (cid:48) | ψ n, k (cid:105) = − i (cid:88) α,α (cid:48) (cid:90) d r d r (cid:48) (cid:104) u m, k | r , α (cid:105)(cid:104) r , α | e − i k · r ( r H − H r (cid:48) ) e i k · r (cid:48) | r (cid:48) , α (cid:48) (cid:105)(cid:104) r (cid:48) , α (cid:48) | u n, k (cid:105) = (cid:104) u m, k | ∂∂ k (cid:18) P (cid:90) d r d r (cid:48) | r , α (cid:105)(cid:104) r , α | e − i k · r He i k · r (cid:48) | r (cid:48) , α (cid:48) (cid:105)(cid:104) r (cid:48) , α (cid:48) | P (cid:19) | u n, k (cid:105) = (cid:104) u m, k | ∂ ˆ H ( k ) ∂ k | u n, k (cid:105) , (117)where we have used the first equation in (116) and the fact that P | u n, k (cid:105) = | u n, k (cid:105) . Therefore, by carrying out the ω (cid:48) integral in Eq. (110), we arrive at σ xy = i e (cid:126) (cid:88) n,m (cid:90) MBZ d k π n F ( (cid:15) n, k ) − n F ( (cid:15) m, k )( (cid:15) m, k − (cid:15) n, k ) (cid:104) ψ n, k | ˆ j x | ψ m, k (cid:105)(cid:104) ψ m, k | ˆ j y | ψ n, k (cid:105) = i e (cid:126) (cid:88) n,m (cid:90) MBZ d k π n F ( (cid:15) n, k ) − n F ( (cid:15) m, k )( (cid:15) m, k − (cid:15) n, k ) (cid:104) u n, k | ∂ ˆ H ( k ) ∂k x | u m, k (cid:105)(cid:104) u m, k | ∂ ˆ H ( k ) ∂k y | u n, k (cid:105) = i e (cid:126) (cid:88) n ∈ occ ,m (cid:54)∈ occ (cid:90) MBZ d k π (cid:104) u n, k | ∂ ˆ H ( k ) ∂k x | u m, k (cid:105)(cid:104) u m, k | ∂ ˆ H ( k ) ∂k y | u n, k (cid:105) − (cid:104) u n, k | ∂ ˆ H ( k ) ∂k y | u m, k (cid:105)(cid:104) u m, k | ∂ ˆ H ( k ) ∂k x | u n, k (cid:105) ( (cid:15) m, k − (cid:15) n, k ) , (118)8where the integral of k runs over one magnetic BZ (MBZ), n F ( (cid:15) ) = − sgn( (cid:15) )2 is the zero temperature Fermi-Diracdistribution function, n ∈ occ means n runs over all the occupied bands, and m (cid:54)∈ occ means m runs over all theempty bands. Since we are in a Hofstadter gap, n F ( (cid:15) n, k ) only depends on the band index n . Then, using the factthat (cid:104) u m, k | ∂ ˆ H ( k ) ∂k i | u n, k (cid:105) = ∂ k i ( (cid:104) u m, k | ˆ H ( k ) | u n, k (cid:105) ) − (cid:104) ∂ k i u m, k | ˆ H ( k ) | u n, k (cid:105) − (cid:104) u m, k | ˆ H ( k ) | ∂ k i u n, k (cid:105) = 0 − (cid:15) n, k (cid:104) ∂ k i u m, k | u n, k (cid:105) − (cid:15) m, k (cid:104) u m, k | ∂ k i u n, k (cid:105) = ( (cid:15) n, k − (cid:15) m, k ) (cid:104) u m, k | ∂ k i u n, k (cid:105) = ( (cid:15) m, k − (cid:15) n, k ) (cid:104) ∂ k i u m, k | u n, k (cid:105) , (119)we can rewrite the Hall conductance as the TKNN formula [10] σ xy = i e (cid:126) (cid:88) n ∈ occ (cid:88) m (cid:54)∈ occ (cid:90) MBZ d k π (cid:0) (cid:104) ∂ k x u n, k | u m, k (cid:105)(cid:104) u m, k | ∂ k y u n, k (cid:105) − (cid:104) ∂ k y u n, k | u m, k (cid:105)(cid:104) u m, k | ∂ k x u n, k (cid:105) (cid:1) = i e (cid:126) (cid:88) n ∈ occ (cid:88) m (cid:90) MBZ d k π (cid:0) (cid:104) ∂ k x u n, k | u m, k (cid:105)(cid:104) u m, k | ∂ k y u n, k (cid:105) − (cid:104) ∂ k y u n, k | u m, k (cid:105)(cid:104) u m, k | ∂ k x u n, k (cid:105) (cid:1) = i e (cid:126) (cid:88) n ∈ occ (cid:90) MBZ d k π (cid:0) (cid:104) ∂ k x u n, k | ∂ k y u n, k (cid:105) − (cid:104) ∂ k y u n, k | ∂ k x u n, k (cid:105) (cid:1) = e (cid:126) (cid:88) n ∈ occ (cid:90) MBZ d k π F nxy ( k ) = t ν e h , (120)where we have used the fact that t ν in the Streda formula (Eq. (109)) is the total Chern number of occupied bands, (cid:80) m | u m, k (cid:105)(cid:104) u m, k | is the identity matrix for the hatted Hamiltonian ˆ H ( k ) at magnetic BZ quasimomentum k , and F nxy ( k ) is the Berry curvature of the n -th band. The corresponding U(1) Berry gauge field is A n ( k ) = i (cid:104) u n, k | ∂ k u n, k (cid:105) , F nxy ( k ) = ∂ x A ny ( k ) − ∂ y A nx ( k ) . (121) B. Quantized Lorentz susceptibility from the Kubo formula
We now try to find a quantity that can be viewed as the momentum space dual of the Hall conductance. The Hallconductance gives a current density j x = ∂H/∂ k x in response to a uniform electric field E y = d k y / d t where k is thecanonical momentum. Therefore, a natural dual response can be defined by exchanging the roles of the real spaceposition r and the momentum k . Such a position-momentum dual response (dual to the Hall conductance response)corresponds to a force density (force per unit cell) F x = − ∂H/∂x acting on the lattice in response to a positionpumping, i.e., velocity v y = d y/ d t of the system, which is nothing but the Lorentz force. Therefore, we name thecoefficient γ xy in the response F x = γ xy v y (122)as the Lorentz susceptibility.We first identify the spatially uniform force operator ˆ F representing the Lorentz force felt by the electrons as thelattice moves rigidly. Recall that an electron in a magnetic field is centered at the guiding center R . In a translationallyinvariant state, if the guiding centers R of the electrons drift at velocity (cid:104) d R d t (cid:105) , the electrons will feel a Lorentz force (cid:104) ˆ F (cid:105) = B ˆ z × (cid:104) d R d t (cid:105) = − iB ˆ z × (cid:104) [ R , H ] (cid:105) = − i (cid:104) [ Π , H ] (cid:105) − iB ˆ z × (cid:104) [ r , H ] (cid:105) , (123)where in the last equality we have used the definition of guiding center R = r − (cid:96) ˆ z × Π . For adiabatic processes,the electrons are in equilibrium, so the Lorentz force (cid:104) ˆ F (cid:105) the electrons felt has to be balanced by a force −(cid:104) ˆ F (cid:105) thelattice exerts on the electrons. According to Newton’s third law, the electrons will then exert a force (cid:104) ˆ F (cid:105) on thelattice. Therefore, experimentally one could measure the force felt by the lattice to find the Lorentz force (cid:104) ˆ F (cid:105) felt bythe electrons.Since the force operator ˆ F is spatially uniform and preserves the translation symmetry, it is diagonal in magnetic BZquasimomentum k . On the right hand side of Eq. (123), the operator [ Π , H ] is readily diagonal in magnetic BZ quasi-momentum k , since Π and hence H both commute with the translation symmetry operators (cid:101) T q d = e − iq(cid:96) − (ˆ z × R ) · d (cid:101) T d = e − i(cid:96) − (ˆ z × R ) · d . The other operator [ r , H ] is not diagonal in k (Ref. [28]), but its diagonal part is simplythe current operator ˆ j as we discussed in Sec. VIA. Therefore, we conclude the force operator is given byˆ F = − i [ Π , H ] + B ˆ z × ˆ j , (124)which is diagonal in magnetic BZ quasimomentum k . Note that here H is the full Hamiltonian defined in Eq. (34)for continuum models or in Eq. (79) for tight-binding models (embedded in the continuum space).In particular, if there is no periodic lattice potential for the electrons at all, H will be solely a function of Π (i.e.,the kinetic term of electrons in the free space), and one will have (cid:104) ˆ F (cid:105) = − iB ˆ z × (cid:104) [ R , H ] (cid:105) ≡ since [ R , Π ] = 0.Physically, this is because without a lattice potential, the electrons cannot feel a force in respond to the movement ofthe lattice and thus cannot exert a force on the lattice.Similar to σ xy which is given by the Kubo formula of the spatially uniform current operator ˆ j , we can define theLorentz susceptibility in Eq. (122) by the Kubo formula of the spatially uniform force operator ˆ F as γ xy = − i ∂∂ω (cid:90) ∞−∞ d ω (cid:48) (cid:104) G ω + ω (cid:48) ˆ F x G ω (cid:48) ˆ F y (cid:105) , (125)where G ω is the Green’s function defined in Eq. (112).In the below, we first calculate the quantized value of the Lorentz susceptibility in Sec. VI B 1, and then derive adual formula for the Lorentz susceptibility analogous to the TKNN formula for Hall conductance in Sec. VI B 2.
1. Calculation of the Lorentz susceptibility from the Diophantine equation
By the expression of ˆ F in Eq. (124), we can rewrite the Lorentz susceptibility defined in Eq. (125) as γ xy = − ie Ω (cid:88) n ∈ occ ,m (cid:54)∈ occ (cid:90) MBZ d k π (cid:15) m, k − (cid:15) n, k ) (cid:104) (cid:104) ψ n, k | ˆ F x | ψ m, k (cid:105)(cid:104) ψ m, k | ˆ F y | ψ n, k (cid:105) − ( x ↔ y ) (cid:105) = − ie Ω (cid:88) n ∈ occ ,m (cid:54)∈ occ (cid:15) m, k − (cid:15) n, k ) (cid:90) MBZ d k π (cid:104) − B (cid:16) (cid:104) ψ n, k | ˆ j y | ψ m, k (cid:105)(cid:104) ψ m, k | ˆ j x | ψ n, k (cid:105) − ( x ↔ y ) (cid:17) − iB (cid:16) (cid:104) ψ n, k | [Π x , H ] | ψ m, k (cid:105)(cid:104) ψ m, k | ˆ j x | ψ n, k (cid:105) − (cid:104) ψ n, k | ˆ j x | ψ m, k (cid:105)(cid:104) ψ m, k | [Π x , H ] | ψ n, k (cid:105)− (cid:104) ψ n, k | [Π y , H ] | ψ m, k (cid:105)(cid:104) ψ m, k | ˆ j y | ψ n, k (cid:105) + (cid:104) ψ n, k | ˆ j y | ψ m, k (cid:105)(cid:104) ψ m, k | [Π y , H ] | ψ n, k (cid:105) (cid:17) − (cid:104) ψ n, k | [Π x , H ] | ψ m, k (cid:105)(cid:104) ψ m, k | [Π y , H ] | ψ n, k (cid:105) + (cid:104) ψ n, k | [Π x , H ] | ψ m, k (cid:105)(cid:104) ψ m, k | [Π y , H ] | ψ n, k (cid:105) (cid:105) , (126)where the integral of k runs over one magnetic BZ (MBZ). Note that so far everything is expressed in the eigenbasis | ψ m, k (cid:105) of the full Hamiltonian H . The first two terms of correlations of ˆ j x and ˆ j y in Eq. (126) simply yield the Berrycurvature, which can be calculated by transforming from the eigenbasis | ψ m, k (cid:105) of the full Hamiltonian H into thereduced basis | u m, k (cid:105) of the reduced magnetic BZ momentum space Hamiltonian ˆ H ( k ) as we have shown in Sec. VIA.We now examine the remaining terms of Eq. (126) (which can be calculated simply in the eigenbasis | ψ m, k (cid:105) of thefull Hamiltonian H ). Since [ Π , H ] is diagonal in k , and ˆ j is the part of operator − i [ r , H ] that is diagonal in k , wecan rewrite the terms of correlations of [Π x , H ] and ˆ j x as (cid:88) n ∈ occ ,m (cid:54)∈ occ (cid:104) ψ n, k | [Π x , H ] | ψ m, k (cid:105)(cid:104) ψ m, k | ˆ j x | ψ n, k (cid:105) ( (cid:15) m, k − (cid:15) n, k ) − (cid:16) [Π x , H ] ↔ ˆ j x (cid:17) = (cid:88) n ∈ occ (cid:88) m (cid:54)∈ occ , k (cid:48) (cid:104) ψ n, k | [Π x , H ] | ψ m, k (cid:48) (cid:105)(cid:104) ψ m, k (cid:48) | ( − i [ x, H ]) | ψ n, k (cid:105) ( (cid:15) m, k (cid:48) − (cid:15) n, k ) − (cid:16) [Π x , H ] ↔ − i [ x, H ] (cid:17) = (cid:88) n ∈ occ (cid:88) m (cid:54)∈ occ , k (cid:48) i ( (cid:15) m, k (cid:48) − (cid:15) n, k ) (cid:104) ψ n, k | Π x | ψ m, k (cid:48) (cid:105)(cid:104) ψ m, k (cid:48) | x | ψ n, k (cid:105) ( (cid:15) m, k (cid:48) − (cid:15) n, k ) − (cid:16) Π x ↔ x (cid:17) = (cid:88) n ∈ occ (cid:88) m, k (cid:48) i (cid:104) ψ n, k | Π x | ψ m, k (cid:48) (cid:105)(cid:104) ψ m, k (cid:48) | x | ψ n, k (cid:105) − (cid:16) Π x ↔ x (cid:17) = (cid:88) n ∈ occ i (cid:104) ψ n, k | [Π x , x ] | ψ n, k (cid:105) , (127)0where we have used the fact that (cid:80) m, k (cid:48) | ψ m, k (cid:48) (cid:105)(cid:104) ψ m, k (cid:48) | is identity for the entire Hilbert space. A similar equality holdsfor the terms of correlations of [Π y , H ] and ˆ j y . Lastly, by the same derivation technique, we have (cid:88) n ∈ occ ,m (cid:54)∈ occ (cid:104) ψ n, k | [Π x , H ] | ψ m, k (cid:105)(cid:104) ψ m, k | [Π y , H ] | ψ n, k (cid:105) ( (cid:15) m, k − (cid:15) n, k ) = − (cid:88) n ∈ occ (cid:104) ψ n, k | [Π x , Π y ] | ψ n, k (cid:105) . (128)Therefore, we can rewrite Eq. (126) as γ xy = e Ω (cid:88) n ∈ occ (cid:90) MBZ d k π (cid:16) − B F nxy ( k ) − iB (cid:104) ψ n, k | [Π x , x ] | ψ n, k (cid:105) + iB (cid:104) ψ n, k | [Π y , y ] | ψ n, k (cid:105) − i (cid:104) ψ n, k | [Π x , Π y ] | ψ n, k (cid:105) (cid:17) = e Ω (cid:88) n ∈ occ (cid:90) MBZ d k π (cid:16) − B F nxy ( k ) + B (cid:17) = eB (cid:16) ρ − t ν ϕ π (cid:17) = eBs ν , (129)where we have used commutation relations [Π x , x ] = [Π y , y ] = − i , [Π x , Π y ] = i(cid:96) − = iB , and the Diophantine equationrewritten in the form of Eq. (109). As we discussed in both the main text and the supplementary Sec. V, ρ = ν/q isthe number of electrons per zero-magnetic-field unit cell. This proves that the Lorentz susceptibility is quantized interms of s ν , and corresponds to a Lorentz force per original (zero field) unit cell F x = γ xy v y = eBs ν v y . (130)
2. A formula for Lorentz susceptibility similar to the TKNN formula for Hall conductance
In Eq. (117), we know the matrix elements of the spatially uniform current operator ˆ j can be calculated by ∂ ˆ H ( k ) ∂ k ,where ˆ H ( k ) is the reduced Hamiltonian defined in Eq. (115). Here for the force operator ˆ F , we can find a similarexpression as follows. First, we define another reduced Hamiltonianˆ (cid:101) H ( d ) = P ( (cid:101) T − d H (cid:101) T d ) P = P ( e i(cid:96) − (ˆ z × R ) · d He − i(cid:96) − (ˆ z × R ) · d ) P , P = qN Ω (cid:88) j ,j ∈ Z (cid:101) T qj d + j d , (131)where (cid:101) T d = (cid:101) T †− d = e − i(cid:96) − (ˆ z × R ) · d (similar to e i k · r ) is the translation operator generated by guiding center R , and P is still the projection operator into the magnetic translationally invariant subspace as defined in Eq. (115). Using theBloch eigenstates | ψ n, k (cid:105) of the full Hamiltonian H , we can define a set of states | w n, d (cid:105) = (cid:101) T − d | ψ n,(cid:96) − ˆ z × d (cid:105) = e i(cid:96) − (ˆ z × R ) · d | ψ n,(cid:96) − ˆ z × d (cid:105) , (132)where | ψ n,(cid:96) − ˆ z × d (cid:105) is the Bloch wave function at momentum k = (cid:96) − ˆ z × d . More discussions and an explicit exampleof the Hamiltonian ˆ (cid:101) H ( d ) and wavefunction | w n, d (cid:105) are given in Sec. VI C. Making use of the fact that (cid:101) T d (cid:101) T d (cid:48) = e − i(cid:96) − ˆ z · ( d × d (cid:48) ) (cid:101) T d (cid:48) (cid:101) T d = e − i(cid:96) − ˆ z · ( d × d (cid:48) ) / (cid:101) T d + d (cid:48) (since [ R x , R y ] = − i(cid:96) ), we have (cid:101) T q d | w n, d (cid:105) = (cid:101) T q d (cid:101) T − d | ψ n,(cid:96) − ˆ z × d (cid:105) = e i(cid:96) − ˆ z · ( q d × d ) (cid:101) T − d (cid:101) T q d | ψ n,(cid:96) − ˆ z × d (cid:105) = e − iq d · ( (cid:96) − ˆ z × d ) (cid:101) T − d e iq d · ( (cid:96) − ˆ z × d ) | ψ n,(cid:96) − ˆ z × d (cid:105) = | w n, d (cid:105) , (cid:101) T d | w n, d (cid:105) = (cid:101) T d (cid:101) T − d | ψ n,(cid:96) − ˆ z × d (cid:105) = e i(cid:96) − ˆ z · ( d × d ) (cid:101) T − d (cid:101) T d | ψ n,(cid:96) − ˆ z × d (cid:105) = e − i d · ( (cid:96) − ˆ z × d ) (cid:101) T − d e i d · ( (cid:96) − ˆ z × d ) | ψ n,(cid:96) − ˆ z × d (cid:105) = | w n, d (cid:105) , (133)where we have used the definition of Bloch momentum in Eq. (105). Therefore, we find the states | w n, d (cid:105) aretranslationally invariant, namely, P | w n, d (cid:105) = | w n, d (cid:105) . (134)It is then easy to see that the eigenstates of the reduced Hamiltonian ˆ (cid:101) H ( d ) are given by | w n, d (cid:105) , namely,ˆ (cid:101) H ( d ) | w n, d (cid:105) = P (cid:101) T − d H (cid:101) T d | w n, d (cid:105) = P (cid:101) T − d H | ψ n,(cid:96) − ˆ z × d (cid:105) = P (cid:101) T − d (cid:15) n,(cid:96) − ˆ z × d | ψ n,(cid:96) − ˆ z × d (cid:105) = (cid:15) n,(cid:96) − ˆ z × d | w n,(cid:96) − ˆ z × d (cid:105) , (135)1where (cid:15) n,(cid:96) − ˆ z × d are the eigenenergies at momentum k = (cid:96) − ˆ z × d . In particular, we can further prove that the forceoperator satisfies (recall Eq. (123)) (cid:104) ψ n,(cid:96) − ˆ z × d | ˆ F | ψ n,(cid:96) − ˆ z × d (cid:105) = − i (cid:104) ψ n,(cid:96) − ˆ z × d | B [ˆ z × R , H ] | ψ n,(cid:96) − ˆ z × d (cid:105) = − i (cid:104) w n, d | (cid:96) − e i(cid:96) − (ˆ z × R ) · d [ˆ z × R , H ] e − i(cid:96) − (ˆ z × R ) · d | w n, d (cid:105) = −(cid:104) w n, d | ∂∂ d (cid:16) P e i(cid:96) − (ˆ z × R ) · d He − i(cid:96) − (ˆ z × R ) · d P (cid:17) | w n, d (cid:105) = −(cid:104) w n, d | ∂ ˆ (cid:101) H ( d ) ∂ d | w n, d (cid:105) . (136)Eq. (136) is then the dual analog to Eq. (117) (see Sec. VI C for an explicit example).We then consider the periodicity of the wave function | w n, d (cid:105) in the space of parameter d . We first note that by thedefinition of the projection operator P in Eq. (131), one has P (cid:101) T d /p = qN Ω (cid:101) T d /p (cid:88) j ,j ∈ Z e − i(cid:96) − ˆ z · [( qj d + j d ) × ( d /p )] (cid:101) T qj d + j d = (cid:101) T d /p P , (137)where we have used the fact that (cid:96) − ˆ z · ( d × d ) = 2 πp/q . Therefore, we findˆ (cid:101) H ( d + d /p ) = P ( (cid:101) T − d − d /p H (cid:101) T d + d /p ) P = P (cid:101) T − d /p (cid:101) T − d H (cid:101) T d (cid:101) T d /p P = (cid:101) T − d /p ( P (cid:101) T − d H (cid:101) T d P ) (cid:101) T d /p = (cid:101) T − d /p ˆ (cid:101) H ( d ) (cid:101) T d /p , (138)namely, the reduced Hamiltonian ˆ (cid:101) H ( d ) and ˆ (cid:101) H ( d + d /p ) differ by a d -independent unitary transformation (cid:101) T d /p .Moreover, the lattice translation symmetry tells us thatˆ (cid:101) H ( d + d ) = P ( (cid:101) T − d − d H (cid:101) T d + d ) P = P (cid:101) T − d (cid:101) T − d H (cid:101) T d (cid:101) T d P = P (cid:101) T − d H (cid:101) T d P = ˆ (cid:101) H ( d ) , (139)where we have used the fact that the full Hamiltonian satisfies [ (cid:101) T d , H ] = 0 (translational symmetry), as given in Eq.(103). Therefore, we conclude that the eigenstates of ˆ (cid:101) H ( d ) satisfies | w n, d + d /p (cid:105) = e iϕ ( d ) (cid:101) T − d /p | w n, d (cid:105) , | w n, d + d (cid:105) = e iϕ ( d ) | w n, d (cid:105) , (140)where ϕ i ( d ) are some d -dependent phases. This effectively defines a dual ”Brillouin zone” Ω M in the d space withperiods d and d /p (which is a torus). The unitary transformation (cid:101) T − d /p serves as an embedding matrix for mappingstate | w n, d (cid:105) to | w n, d + d /p (cid:105) . Note that Ω M is simple 1 /p fraction of the zero magnetic field unit cell. Furthermore,this allows us to define a U(1) Berry gauge field and its field strength for band n on the dual ”Brillouin zone” Ω m (which is a closed manifold): (cid:101) a n ( d ) = − i (cid:104) w n, d | ∂ d w n, d (cid:105) , (cid:101) f nxy ( d ) = ∂ x (cid:101) a ny ( d ) − ∂ y (cid:101) a nx ( d ) . (141)In the following, we will show that the Lorentz susceptibility is given by a formula similar to the TKNN formula [10],but with the usual Berry gauge field in Eq. (121) replaced by the new gauge field we defined in Eq. (141).We first note that the Lorentz susceptibility in Eq. (126) can be re-expressed in the d space (by variable substitution k = (cid:96) − ˆ z × d ) as γ xy = − ie Ω (cid:88) n ∈ occ ,m (cid:54)∈ occ (cid:90) MBZ d k π (cid:15) m, k − (cid:15) n, k ) (cid:104) (cid:104) ψ n, k | ˆ F x | ψ m, k (cid:105)(cid:104) ψ m, k | ˆ F y | ψ n, k (cid:105) − ( x ↔ y ) (cid:105) = − i eBp πq (cid:88) n ∈ occ ,m (cid:54)∈ occ (cid:90) d ∈M { q d /p, d /p } d d (cid:104) (cid:104) ψ n,(cid:96) − ˆ z × d | ˆ F x | ψ m,(cid:96) − ˆ z × d (cid:105)(cid:104) ψ m,(cid:96) − ˆ z × d | ˆ F y | ψ n,(cid:96) − ˆ z × d (cid:105) − ( x ↔ y ) (cid:105) ( (cid:15) m,(cid:96) − ˆ z × d − (cid:15) n,(cid:96) − ˆ z × d ) , (142)where Ω is the zero magnetic field unit cell area, MBZ stands for the magnetic BZ spanned by g /q and g , and M { q d /p, d /p } denotes a torus with periods q d /p and d /p , which is nothing but the MBZ rotated by π/ (cid:96) (this can be seen by noting that d = Ω2 π ˆ z × g , and d = − Ω2 π ˆ z × g , and B Ω = (cid:96) − Ω = 2 πp/q ).2Then, by Eq. (107) we know that (cid:101) T d | ψ n, k (cid:105) is the same as | ψ n, k − p g /q (cid:105) up to a phase factor. Since p and q arecoprime, there exists an integer j so that jp ≡ − q ), and thus (cid:101) T j d | ψ n, k (cid:105) = e iα k | ψ n, k + jp g /q (cid:105) = e iα (cid:48) k | ψ n, k − g /q (cid:105) = e iα (cid:48) k | ψ n, k + (cid:96) − ˆ z × d /p (cid:105) , where α k and α (cid:48) k are some phase factors. Making use of the fact that [ ˆ F , (cid:101) T d ] = 0 (translationallyinvariant), we then have (cid:104) ψ n,(cid:96) − ˆ z × d | ˆ F | ψ m,(cid:96) − ˆ z × d (cid:105) = (cid:104) ψ n,(cid:96) − ˆ z × d | (cid:101) T − j d ˆ F (cid:101) T j d | ψ m,(cid:96) − ˆ z × d (cid:105) = (cid:104) ψ n,(cid:96) − ˆ z × ( d + d /p ) | ˆ F | ψ m,(cid:96) − ˆ z × ( d + d /p ) (cid:105) , (143)namely, the quantity (cid:104) ψ n,(cid:96) − ˆ z × d | ˆ F | ψ m,(cid:96) − ˆ z × d (cid:105) invariant under displacement d → d + d /p . Therefore, the integralin region M { q d /p, d /p } in Eq. (142) is equal to the same integral in a region M { d , d /p } times a global factor q/p .Note that region M { d , d /p } is nothing but the torus dual ”Brillouin zone” Ω M on which we defined the U(1) Berrygauge field in Eq. (141). We can then rewrite the Lorentz susceptibility in Eq. (142) as γ xy = − i eB π (cid:88) n ∈ occ ,m (cid:54)∈ occ (cid:90) d ∈ Ω M d d (cid:104) (cid:104) ψ n,(cid:96) − ˆ z × d | ˆ F x | ψ m,(cid:96) − ˆ z × d (cid:105)(cid:104) ψ m,(cid:96) − ˆ z × d | ˆ F y | ψ n,(cid:96) − ˆ z × d (cid:105) − ( x ↔ y ) (cid:105) ( (cid:15) m,(cid:96) − ˆ z × d − (cid:15) n,(cid:96) − ˆ z × d ) = − i eB π (cid:88) n ∈ occ ,m (cid:54)∈ occ (cid:90) d ∈ Ω M d d (cid:20) (cid:104) ψ n,(cid:96) − ˆ z × d | ∂ ˆ (cid:101) H ( d ) ∂d x | ψ m,(cid:96) − ˆ z × d (cid:105)(cid:104) ψ m,(cid:96) − ˆ z × d | ∂ ˆ (cid:101) H ( d ) ∂d y | ψ n,(cid:96) − ˆ z × d (cid:105) − ( x ↔ y ) (cid:21) ( (cid:15) m,(cid:96) − ˆ z × d − (cid:15) n,(cid:96) − ˆ z × d ) , (144)where we have used Eq. (136). Then, by noting that (cid:104) w m, d | ∂ ˆ (cid:101) H ( d ) ∂d i | w n, d (cid:105) = ∂ d i ( (cid:104) w m, d | ˆ (cid:101) H ( d ) | w n, d (cid:105) ) − (cid:104) ∂ d i w m, d | ˆ (cid:101) H ( d ) | w n, d (cid:105) − (cid:104) w m, d | ˆ (cid:101) H ( d ) | ∂ d i w n, d (cid:105) = 0 − (cid:15) n,(cid:96) − ˆ z × d (cid:104) ∂ d i w m,d | w n, d (cid:105) − (cid:15) m,(cid:96) − ˆ z × d (cid:104) w m, d | ∂ d i w n, d (cid:105) = ( (cid:15) n,(cid:96) − ˆ z × d − (cid:15) m,(cid:96) − ˆ z × d ) (cid:104) w m, d | ∂ d i w n, d (cid:105) = ( (cid:15) m,(cid:96) − ˆ z × d − (cid:15) n,(cid:96) − ˆ z × d ) (cid:104) ∂ k i w m, d | w n, d (cid:105) , (145)we can then derive the following: γ xy = − i eB π (cid:88) n ∈ occ ,m (cid:54)∈ occ (cid:90) d ∈ Ω M d d (cid:0) (cid:104) ∂ d x w n, d | w m, d (cid:105)(cid:104) w m, d | ∂ d y w n, d (cid:105) − (cid:104) ∂ d y w n, d | w m, d (cid:105)(cid:104) w m, d | ∂ k x w n, d (cid:105) (cid:1) = − i eB π (cid:88) n ∈ occ ,m (cid:90) d ∈ Ω M d d (cid:0) (cid:104) ∂ d x w n, d | w m, d (cid:105)(cid:104) w m, d | ∂ d y w n, d (cid:105) − (cid:104) ∂ d y w n, d | w m, d (cid:105)(cid:104) w m, d | ∂ d x w n, d (cid:105) (cid:1) = − i eB π (cid:88) n ∈ occ (cid:90) d ∈ Ω M d d (cid:0) (cid:104) ∂ d x w n, d | ∂ d y w n, d (cid:105) − (cid:104) ∂ d y w n, d | ∂ d x w n, d (cid:105) (cid:1) = eB (cid:88) n ∈ occ (cid:90) d ∈ Ω M d d π (cid:101) f nxy ( d ) , (146)where (cid:101) f nxy ( d ) is the new Berry curvature in the d space we defined in Eq. (141). This is then a dual analogy to theTKNN formula. In particular, our results for γ xy in Sec. VI B 1 implies (cid:88) n ∈ occ (cid:90) d ∈ Ω M d d π (cid:101) f nxy ( d ) = s ν . (147)Therefore, s ν can be viewed as a dual Chern number, where the integral is within a torus “dual magnetic Brillouinzone” Ω M with periods d and d /p .— Eq. (147) can also be proved by the gauge choices of the Bloch wave functions in Eqs. (106) and (107) (which weemployed to prove the Diophantine equation), namely, one can fix | ψ n, k + g /q (cid:105) = | ψ n, k (cid:105) , | ψ n, k + g (cid:105) = e − iσq d · k | ψ n, k (cid:105) and (cid:101) T d | ψ n, k (cid:105) = e imq d · k | ψ n, k + p g /q (cid:105) for a Hofstadter band n between the ( ν − ν -th gap, where σ = t ν − t ν − and m = s ν − s ν − . Then by Eq. (132), we find | w n, d + d /p (cid:105) = (cid:101) T − d − d /p | ψ n,(cid:96) − ˆ z × ( d + d /p ) (cid:105) = e i(cid:96) − ˆ z · ( d × d ) / p (cid:101) T − d /p (cid:101) T − d | ψ n,(cid:96) − ˆ z × d − g /q (cid:105) = e i(cid:96) − ˆ z · ( d × d ) / p (cid:101) T − d /p (cid:101) T − d | ψ n,(cid:96) − ˆ z × d (cid:105) = e i(cid:96) − ˆ z · ( d × d ) / p (cid:101) T − d /p | w n, d (cid:105) , | w n, d + d (cid:105) = (cid:101) T − d − d | ψ n,(cid:96) − ˆ z × ( d + d ) (cid:105) = e − i(cid:96) − ˆ z · ( d × d ) / (cid:101) T − d (cid:101) T − d | ψ n,(cid:96) − ˆ z × d + p g /q (cid:105) = e − i(cid:96) − ˆ z · ( d × d ) / − imq d · ( (cid:96) − ˆ z × d ) (cid:101) T − d | ψ n,(cid:96) − ˆ z × d (cid:105) = e − i(cid:96) − ˆ z · ( d × d ) / imq(cid:96) − ˆ z · ( d × d ) | w n, d (cid:105) . (148)3Therefore, for band n between the ( ν − ν -th gap, we find the dual Berry gauge field in Eq. (141)satisfies (cid:90) d ∈ Ω M d d π (cid:101) f nxy ( d ) = 12 π (cid:73) ∂ Ω M (cid:101) a n ( d ) · d d = 12 π [ mq(cid:96) − ˆ z · ( d × d /p ) − (cid:96) − ˆ z · ( d × d /p ) / (cid:96) − ˆ z · ( d × d ) / p ] = mq πp (cid:96) − Ω= m = s ν − s ν − . (149)Summing over all occupied bands then proves Eq. (147). C. Relation between the dual wave functions, and a model example
This subsection is devoted for a better understanding of the relation between the two wave functions | u n, k (cid:105) definedin Eq. (116) and | w n, d (cid:105) defined in Eq. (132). To do this, we examine their form in the real space basis | r (cid:105) . If wedenote the Bloch wave function of momentum k as | ψ n, k (cid:105) , by the definition of the periodic wavefunction | u n, k (cid:105) in Eq.(116), we have ψ n, k ( r ) = (cid:104) r | ψ n, k (cid:105) , u n, k ( r ) = (cid:104) r | u n, k (cid:105) = e − i k · r ψ n, k ( r ) . (150)To find the real space wave function of | w n, d (cid:105) = (cid:101) T − d | ψ n,(cid:96) − ˆ z × d (cid:105) as defined in Eq. (132), we first note that (cid:101) T d | r (cid:105) = e − i(cid:96) − (ˆ z × R ) · d | r (cid:105) = e − i(cid:96) − (ˆ z × r ) · d − i Π · d | r (cid:105) = e i (cid:82) r + d r [ A ( r (cid:48) ) − (cid:96) − ˆ z × r (cid:48) ] · d r (cid:48) | r + d (cid:105) = e i(cid:96) − (ˆ z × d ) · ( r + d ) e i (cid:82) r + d r A ( r (cid:48) ) · d r (cid:48) | r + d (cid:105) , (151)which can be seen by following a derivation similar to Eq. (92), where the integral on the exponent is along thestraight line segment from r to r + d . Therefore, we find w n, d ( r ) = (cid:104) r | w n, d (cid:105) = (cid:104) r | (cid:101) T − d | ψ n,(cid:96) − ˆ z × d (cid:105) = e − i(cid:96) − (ˆ z × d ) · ( r + d ) e − i (cid:82) r + d r A ( r (cid:48) ) · d r (cid:48) (cid:104) r + d | ψ n,(cid:96) − ˆ z × d (cid:105) = e − i (cid:82) r + d r A ( r (cid:48) ) · d r (cid:48) u n,(cid:96) − ˆ z × d ( r + d ) . (152)This is the explicit form of the wave function w n, d ( r ).Without going into the real space position basis, we can also derive that (here r below denotes the position operator) | w n, d (cid:105) = (cid:101) T − d | ψ n,(cid:96) − ˆ z × d (cid:105) = e i(cid:96) − (ˆ z × R ) · d e i(cid:96) − (ˆ z × d ) · r | u n,(cid:96) − ˆ z × d (cid:105) = e i Π · d | u n,(cid:96) − ˆ z × d (cid:105) , (153)where we have used the fact that [(ˆ z × R ) · d , (ˆ z × r ) · d ] = 0, and R − r = − (cid:96) ˆ z × Π . By this relation, the dual Berrycurvature (cid:101) f nxy ( d ) defined in Eq. (141) can be explicitly related to the Berry curvature F nxy ( k ) in Eq. (121). We have | ∂ d i w n, d (cid:105) = ( ∂ d i e i Π · d ) | u n,(cid:96) − ˆ z × d (cid:105) + e i Π · d | ∂ d i u n,(cid:96) − ˆ z × d (cid:105) = e i Π · d [ i (Π i + (cid:15) ij d j / | u n,(cid:96) − ˆ z × d (cid:105) + | ∂ d i u n,(cid:96) − ˆ z × d (cid:105) ] , (154)where (cid:15) ji ( i, j = x, y ) is the Levi-Civita tensor ( (cid:15) xy = − (cid:15) yx = 1). Therefore, we find the dual Berry curvature (cid:101) f nxy ( d ) = − i (cid:104) ∂ d x w n, d | ∂ d y w n, d (cid:105) + i (cid:104) ∂ d y w n, d | ∂ d x w n, d (cid:105) = − i (cid:104) u n,(cid:96) − ˆ z × d | [Π x , Π y ] | u n,(cid:96) − ˆ z × d (cid:105) − i (cid:104) ∂ d x u n,(cid:96) − ˆ z × d | ∂ d y u n,(cid:96) − ˆ z × d (cid:105) + i (cid:104) ∂ d y u n,(cid:96) − ˆ z × d | ∂ d x u n,(cid:96) − ˆ z × d (cid:105)− (cid:104) ∂ d y u n,(cid:96) − ˆ z × d | Π x | u n,(cid:96) − ˆ z × d (cid:105) − (cid:104) u n,(cid:96) − ˆ z × d | Π x | ∂ d y u n,(cid:96) − ˆ z × d (cid:105) + (cid:104) ∂ d x u n,(cid:96) − ˆ z × d | Π y | u n,(cid:96) − ˆ z × d (cid:105) + (cid:104) u n,(cid:96) − ˆ z × d | Π y | ∂ d x u n,(cid:96) − ˆ z × d (cid:105) = (cid:96) − − (cid:96) − F nxy ( (cid:96) − ˆ z × d ) + ∂ d x ( (cid:104) u n,(cid:96) − ˆ z × d | Π y | u n,(cid:96) − ˆ z × d (cid:105) ) − ∂ d y ( (cid:104) u n,(cid:96) − ˆ z × d | Π x | u n,(cid:96) − ˆ z × d (cid:105) )= B − B F nxy ( (cid:96) − ˆ z × d ) + ∂ d x ( (cid:104) w n, d | Π y | w n, d (cid:105) ) − ∂ d y ( (cid:104) w n, d | Π x | w n, d (cid:105) ) . (155)This gives the relationship between the dual Berry curvature (cid:101) f nxy ( d ) and the Berry curvature F nxy ( k ). In particular, wenote that the two terms ∂ d x ( (cid:104) w n, d | Π y | w n, d (cid:105) ) − ∂ d y ( (cid:104) w n, d | Π x | w n, d (cid:105) ) in the last line of Eq. (155) are total derivativeswhich do not contribute to the dual Chern number s ν , since (cid:104) w n, d | Π i | w n, d (cid:105) are well-defined periodic functions (physicalquantities) in the dual Brillouin zone Ω M spanned by d and d /p .4
1. An example
For simplicity, we consider a single-orbital tight-binding model on a square lattice with a magnetic flux per unitcell ϕ = 2 π/
3. Assume the lattice vectors are d = (1 ,
0) and d = (0 , | D n ,n (cid:105) are deltafunctions at lattice sites D n ,n = ( n , n ) ∈ d Z + d Z in the continuous space. At zero magnetic field, we assumeall the nearest bonds have hopping − t (with t being a real number), while all the longer range hoppings are zero.We then add a uniform magnetic field B = 2 πp/q = 4 π/ ϕ = B | d × d | = 4 π/ A ( r ) = (0 , Bx ). This correspond to p = 2 and q = 3. Accordingly, the magnetic length (cid:96) satisfies (cid:96) − = 4 π/ B can then be written as H = − t (cid:88) n ,n ∈ Z (cid:16) | D n +1 ,n (cid:105)(cid:104) D n ,n | + e i πn / | D n ,n +1 (cid:105)(cid:104) D n ,n | + h.c. (cid:17) . (156)Therefore, we have e − i k · r He i k · r = − t (cid:88) n ,n ∈ Z (cid:16) e − ik x | D n +1 ,n (cid:105)(cid:104) D n ,n | + e − ik y + i πn / | D n ,n +1 (cid:105)(cid:104) D n ,n | + h.c. (cid:17) , (157)where r is the position operator (instead of a parameter), while k = ( k x , k y ) is the momentum parameter. By Eq.(151), the projector P in Eq. (131) in the Hilbert space of Wannier basis | D n ,n (cid:105) can be rewritten as P = 3 N Ω (cid:88) j ,j ∈ Z (cid:101) T j d + j d = (cid:88) j =1 | j (cid:105)(cid:104) j | , | j (cid:105) = (cid:114) N Ω (cid:88) n ,n ∈ Z | D n + j − ,n (cid:105) , ( j = 1 , , . (158)The three basis | j (cid:105) ( j = 1 , ,
3) correspond exactly to the three orbitals in a magnetic unit cell. Therefore, under thethree basis | j (cid:105) , we find the reduced momentum space Hamiltonian ˆ H ( k ) is given byˆ H ( k ) = P e − i k · r He i k · r P = − t k y ) e ik x e − ik x e − ik x k y + 2 π/ e ik x e ik x e − ik x k y − π/ . (159)This is exactly the momentum space Hamiltonian we are familiar with.We now consider the other reduced Hamiltonian ˆ (cid:101) H ( d ). By Eq. (151), for d = ( d x , d y ) we have (cid:101) T − d H (cid:101) T d = − t (cid:88) n ,n ∈ Z (cid:16) | D n +1 ,n − d (cid:105)(cid:104) D n ,n − d | + e − i πd x / i πn / | D n ,n +1 − d (cid:105)(cid:104) D n ,n − d | + h.c. (cid:17) , (160)where we have used A ( r ) = (0 , Bx ) in the continuous space. We note that the Hamiltonian is transformed into ashifted Wannier basis | D n +1 ,n − d (cid:105) . Note that by Eq. (151), the lattice translation operators act as (cid:101) T d | D n ,n − d (cid:105) = e − i πd y | D n +3 ,n − d (cid:105) , (cid:101) T d | D n ,n − d (cid:105) = | D n ,n +1 − d (cid:105) . (161)Therefore, the projector P in the Hilbert space of Wannier basis | D n +1 ,n − d (cid:105) is P = 3 N Ω (cid:88) j ,j ∈ Z (cid:101) T j d + j d = (cid:88) j =1 | j, d (cid:105)(cid:104) j, d | , | j, d (cid:105) = (cid:114) N Ω (cid:88) n ,n ∈ Z e − i π (3 n + j ) d y / | D n + j − ,n − d (cid:105) , (162)where j = 1 , ,
3. Under the basis | j, d (cid:105) , we findˆ (cid:101) H ( d ) = P (cid:101) T − d H (cid:101) T d P = − t πd x / e − i πd y / e i πd y / e i πd y / πd x / π/ e − i πd y / e − i πd y / e i πd y / πd x / − π/ . (163)Therefore, under the Landau gauge, we find ˆ (cid:101) H ( d ) takes a form analogous to ˆ H ( k ) with k = (4 π/ z × d , except thatthe basis | j, d (cid:105) is d dependent. Therefore, if the eigenstates of ˆ H ( k ) are | u n, k (cid:105) = (cid:80) j =1 u n,j ( k ) | j (cid:105) , the eigenstates ofˆ (cid:101) H ( d ) will be | w n, d (cid:105) = (cid:80) j =1 u n,j ( π ˆ z × d ) | j, d (cid:105) , in agreement with Eq. (152).The Chern numbers of the 3 bands of ˆ H ( k ) can be calculated to be σ n = {− , , − } (from the lowest band tothe highest band). According to Eq. (155), we find the dual Chern numbers of the 3 bands of ˆ (cid:101) H ( d ) are given by m n = − σ n = { , − , } (from the lowest band to the highest band).5 - d - dd d A ( r ) FIG. 7. An example of tight binding model on a square lattice in a magnetic field, where Landau gauge is chosen.
D. Physical understanding of the quantized Lorentz susceptibility
The quantized Lorentz susceptibility γ xy = eBs ν can be easily understood from the following physical picture.Consider a lattice model in the x - y plane moving with a velocity v = v y ˆ y in the ˆ z direction magnetic field B , andthe Fermi level is in a Hofstadter gap ν . In the rest frame of the lattice system, the magnetic field B will move witha velocity − v = − v y ˆ y , which produces an electric field E = v y B ˆx . This electric field then produces a Hall currentdensity (cid:101) j y = − σ xy v y B in the rest frame of the lattice. If we go back to the laboratory frame where the system moveswith velocity v = v y ˆ y , we would find a total current density j y = eρv y Ω + (cid:101) j y = (cid:16) eρ Ω − σ xy B (cid:17) v y , where ρ is the number of occupied electron per unit cell, and Ω is the unit cell area. Therefore, the Lorentz force perunit cell is given by F x = j y Ω B = ( eρ − σ xy Ω B ) Bv y = (cid:18) eρ − t ν e h Ω B (cid:19) Bv y = eB (cid:16) ρ − t ν ϕ π (cid:17) v y = eBs ν v y , (164)where we have used the Diophantine equation. Therefore, we find the Lorentz susceptibility is γ xy = eBs νν