Science, Faculty of
Physics and Astronomy, Department of
DSpace
UBCV
Egan, Shannon
Master of Science - MSc
University of British Columbia
Twisted moiré materials have opened up a new platform for realizing the quantum anomalous Hall (QAH) effect, with experiments reporting its observation in both twisted graphene and transition metal dichalcogenide heterostructures. A fascinating feature of these moiré QAH devices is the ability to switch the direction of Hall resistivity with very small magnetic fields or currents, making them a candidate for stable and programmable magnetic memory. Describing this phenomenon theoretically requires knowledge of both the orbital magnetization and orbital susceptibility, of which the latter is less well understood. This motivated us to derive a numerical formula for the zero-field orbital susceptibility of a general multi-band moiré continuum model, expressed in terms of quantities which reflect the geometry of the moiré bands. We use this formula to study the conditions for which magnetic switching can occur in moiré bands of twisted homobilayer MoS₂. We predict that the conduction band states of twisted MoS₂ give rise to two pairs of flat Chern bands in each valley: |C| = ±1 bands generated by moiré potential alone; and |C| = ±2 bands generated by a combination of moiré potential and spin-orbit coupling. All of the flat bands carry a large orbital magnetization ~1 Bohr magneton per moiré unit cell. Our Hartree-Fock analysis shows that at 1/2-filling, interactions generically drive the system into a spontaneous time-reversal broken valley-polarized state, yielding an "orbital Chern insulator" with quantized Hall resistivity. By calculating the Fermi sea contribution to orbital susceptibility, we show that the paramagnetism necessary to induce magnetic switching in the orbital Chern insulator state is confined to the gap between the |C| = ±2 Chern bands. The chemical potential and temperature dependence of the orbital susceptibility may help explain outstanding puzzles in the behaviour of moiré quantum Hall devices.
Orbital susceptibility and magneticswitching in twisted moire´heterostructuresbyShannon EganB.Sc., McGill University, 2020A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinThe Faculty of Graduate and Postdoctoral Studies(Physics)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)August 2022© Shannon Egan 2022The following individuals certify that they have read, and recommend to the Faculty ofGraduate and Postdoctoral Studies for acceptance, the thesis entitled:Orbital susceptibility and magnetic switching in twisted moire´ heterostructuressubmitted by Shannon Egan in partial fulfillment of the requirements forthe degree of Master of Sciencein PhysicsExamining Committee:Marcel Franz, Professor, Physics, UBCSupervisorMona Berciu, Professor, Physics, UBCSupervisory Committee MemberiiAbstractTwisted moire´ materials have opened up a new platform for realizing the quantum anoma-lous Hall (QAH) effect, with experiments reporting its observation in both twisted grapheneand transition metal dichalcogenide heterostructures. A fascinating feature of these moire´QAH devices is the ability to switch the direction of Hall resistivity with very small mag-netic fields or currents, making them a candidate for stable and programmable magneticmemory. Describing this phenomenon theoretically requires knowledge of both the orbitalmagnetization and orbital susceptibility, of which the latter is less well understood. Thismotivated us to derive a numerical formula for the zero-field orbital susceptibility of ageneral multi-band moire´ continuum model, expressed in terms of quantities which reflectthe geometry of the moire´ bands. We use this formula to study the conditions for whichmagnetic switching can occur in moire´ bands of twisted homobilayer MoS2. We predictthat the conduction band states of twisted MoS2 give rise to two pairs of flat Chern bandsin each valley: C = ±1 bands generated by moire´ potential alone; and C = ±2 bandsgenerated by a combination of moire´ potential and spin-orbit coupling. All of the flatbands carry a large orbital magnetization ∼ 1µB per moire´ unit cell. Our Hartree-Fockanalysis shows that at 1/2-filling, interactions generically drive the system into a spon-taneous time-reversal broken valley-polarized state, yielding an “orbital Chern insulator”with quantized Hall resistivity. By calculating the Fermi sea contribution to orbital sus-ceptibility, we show that the paramagnetism necessary to induce magnetic switching inthe orbital Chern insulator state is confined to the gap between the C = ±2 Chern bands.The chemical potential and temperature dependence of the orbital susceptibility may helpexplain outstanding puzzles in the behaviour of moire´ quantum Hall devices.iiiLay SummaryThe ability to isolate a single layer of crystal in the laboratory, creating an effectivelytwo-dimensional (2D) material, has proven an important advancement in both materialsscience and fundamental physics. Experimentalists now have the ability to stack and twistisolated monolayers at extremely precise angles (∼ 0.1 degree), leading to unprecedentedcontrol over properties of layered materials. Here we propose a new model for a twistedbilayer of the 2D hexagonal material MoS2, where the spin of conduction electrons playsan important role, and characterize its magnetic properties. Unusually strong interactionsbetween electrons change the ground state of the system, leading to an exotic Quantum“Anomalous” Hall state and magnetism which remains stable without any external mag-netic field. Understanding this phenomenon could pave the way to engineering magneticmemory devices in twisted moire´ materials.ivPrefaceThe analysis presented in Chapter 2, as well as Appendix A-B referenced therein, pre-viously appeared in a published paper [B. T. Zhou, S. Egan, and M. Franz, “Moire´ flatChern bands and correlated quatum anomalous Hall states generated by spin-orbit cou-plings in twisted homobilayer MoS2,” Phys. Rev. Research, vol. 4, p. L012032, Mar2022]. My collaborator Dr. Benjamin T. Zhou designed the research program for this pa-per, conceived of the model (Section 2.1), calculated the band structure and Chern number(Section 2.2), performed the Hartree-Fock calculation (Section 2.3), and wrote the initialmanuscript. All figures appearing in this Chapter also appear in the paper. I reproducedand verified the results in Figure 2.2, and made significant edits to the manuscript textalongside Professor Marcel Franz.Section 3.1, and subsection 3.2.1 of Section 3.2 reproduce theory and derivations fromthe literature on orbital magnetism, with all works appropriately cited within. The gauge-invariant and numerically tractable formula for orbital susceptibility of a separable bandmodel derived in Section 3.2.2 is an original result, although it is derived using a moregeneral formula and methods from the literature. In contrast to other works, I took specialcare to discuss the numerical implementation of the susceptibility formula, which we foundlacking in existing literature.All results presented in Chapter 4.3 are the original work of the author, Shannon Egan,based on the model (presented in Section 2.1) by Benjamin T. Zhou, who also providedguidance alongside my supervisor Marcel Franz. The work in this Chapter is unpublished.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 Correlated quantum anomalous Hall states in twisted homobilayer MoS2 42.1 Continuum model for twisted homobilayer MoS2 . . . . . . . . . . . . . . . 42.2 Moire´ flat Chern bands and skyrmionic spin textures . . . . . . . . . . . . 62.3 Correlated QAH state at 1/2 filling . . . . . . . . . . . . . . . . . . . . . . 83 Theory of orbital magnetism . . . . . . . . . . . . . . . . . . . . . . . . . . . 123.1 Orbital magnetization of crystalline systems . . . . . . . . . . . . . . . . . 123.2 Orbital susceptibility . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143.2.1 Gauge-invariant Green’s function formalism . . . . . . . . . . . . . 153.2.2 Multi-band formula for orbital susceptibility . . . . . . . . . . . . . 184 Orbital magnetism in topological moire´ bands . . . . . . . . . . . . . . . 234.1 Large out-of-plane magnetization in twisted MoS2 . . . . . . . . . . . . . . 234.2 Ginzburg-Landau theory for magnetic switching . . . . . . . . . . . . . . . 284.3 Orbital susceptibility of the Fermi sea . . . . . . . . . . . . . . . . . . . . . 305 Discussion and outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35viTable of ContentsAppendicesA Calculation of effective inter-layer tunneling parameter w0 . . . . . . . 42B Self-consistent Hartree-Fock calculations . . . . . . . . . . . . . . . . . . . 46B.1 Hartree-Fock energy functional and Lagrange equations . . . . . . . . . . . 46B.2 Mean-field approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48B.3 Hartree-Fock ground states . . . . . . . . . . . . . . . . . . . . . . . . . . . 49B.4 Physically accessible regimes of (J0, J1) . . . . . . . . . . . . . . . . . . . . 51C Gauge-invariant orbital magnetic moment and Berry curvature . . . . 53D Derivation of numerical formula for orbital susceptibility . . . . . . . . 54D.1 Tr′ Xˆ in the band basis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54D.2 Energy integral . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56D.3 Numerical evaluation of χorb . . . . . . . . . . . . . . . . . . . . . . . . . . 60viiList of Tables2.1 Parameters used inH0 for twisted homobilayer MoS2 with monolayer latticeconstant a = 3.16A˚. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6viiiList of Figures2.1 Schematic crystal structure of twisted bilayer MoS2 with (a) top view and(b) side view. Electrons in different layers experience opposite effectiveelectric fields due to local lack of mirror symmetry. (c) Moire´ Brillouinzone (MBZ) formed by twisting hexagonal Brillouin zones of each layer.The in-plane spinor field exhibits opposite vorticities at Km,+ and Km,− ofMBZ due to layer-dependent Rashba effects. . . . . . . . . . . . . . . . . . . 52.2 Moire´ bands of valley ξ = + formed by CB states of twisted homobilayerMoS2 at θ0 = 1.4◦ with (a) αso = 0 and (b) αso ̸= 0. The colorbars indicatethe out-of-plane spin expectation value ⟨Sz⟩ in units of ℏ/2. (c) Momentum-space profile of Berry curvature Ωz in units of A˚2 in the 2nd moire´ bandin (b). Large Ωz of order 104A˚2 emerges as a result of the nontrivial gapinduced by SOC. (d) Momentum-space spin texture of the 2nd moire´ bandin (b). White arrows indicate the in-plane spinor field, which has oppositevorticities at Km,+ and Km,− and leads to a nontrivial skyrmion numberS = +2. Parameters used are presented in Table 2.1. . . . . . . . . . . . . . 72.3 (a) J0 − J1 phase diagram of the system with the 2nd moire´ band at 1/2-filling. Under repulsive density-density interactions, only the J1 ≤ J0regime can be accessed (see Appendix B.4), where the valley-polarized(VP)state is generically favored and the system becomes a correlated QAH insu-lator. (b) Evolutions of the order parameters ∆+− and ∆M ≡ ∆+−∆− asthe ratio J0/J1 is tuned along the line-cut (green dashed line in (a)) acrossthe phase boundary. A first-order transition happens at the critical pointJ0 ≃ J1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92.4 (a)-(b) Non-interacting bands of twisted homobilayer MoS2 for (a) valleyξ = +, and (b) valley ξ = − with Rashba strength αso = 30 meV· A˚.Dashed lines in (a)-(b) indicate location of the Fermi energy EF when the2nd moire´ bands are 1/2-filled and interaction effects are neglected. In theabsence of correlations, both bands from valley ξ = + and valley ξ = − arehalf-filled. (c)-(d) Correlated band structures when the 2nd moire´ bandsare half-filled, where the system enters the valley-polarized state. Electroncorrelations induce an insulating gap ∆CQAH ∼ J0, where we set J0 = 10meVs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11ixList of Figures4.1 Distribution of out-of-plane mtotnk = mnk +mΩnk in the MBZ for the ξ = +1moire´ flat bands labelled in Fig.2.2(b). The chemical potential µ is set tothe top of each band. While the exact colour scale varies for each band, zerois always represented by a pale yellow, and the range is set to [−cmaxn , cmaxn ],where cmaxn = maxk∈MBZ∣∣mtotnk∣∣. . . . . . . . . . . . . . . . . . . . . . . . . 244.2 Comparison of bulk mnk and edge mΩnk magnetic moment contributions in(a) C = +2 and (b) C = −2 bands of valley ξ = +1. The chemical potentialµ is set to the top of each band, and the colour scale is defined in the sameway as Fig.4.1. M bulkz (Medgez ) is defined as the average of mznk (mz,Ωnk ) overthe MBZ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 254.3 Mz vs. µ for each moire´ Chern band (ξ = +1) identified in Fig2.2(b) atT = 0.2K. M bulkz and Medgez are defined as in Fig.4.2. The magnetizationplotted here represents the contribution from the states of a single band, notthe total magnetization of all filled states. The magnetization contributedby the bulk term plateaus when µ is above the band edge (grey shadedregion), while the edge contribution continues to scale linearly in proportionto Chern number. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 274.4 Fermi sea orbital susceptibility χsea vs. µ for non-interacting moire´ bandsof valley ξ = +1 at different temperatures. Each point incorporates thecontribution of all states which are filled at a given µ, by summing overeach k ∈ MBZ and band n. Grey shaded regions indicate the width ofeach band, with the slight overlap between the 2nd and 3rd band (C = ±2)appearing in darker grey in panel (a) and (c). (a)-(b) and (c)-(d) containthe same data, but we artificially increase the energy of the 3rd and 4thband by 1 meV in to highlight the paramagnetic susceptibility at the top ofthe 2nd band. The width and strength of the paramagnetic region is highlysensitive to temperature and gap size. . . . . . . . . . . . . . . . . . . . . . 32A.1 Interlayer tunneling strength t(q) as a function of q calculated by numericalFourier transform. Dashed vertical line in black indicates the value w ≈ 10meV extrapolated at q = K where Ka/π = 4/3. . . . . . . . . . . . . . . . . 44xList of Abbreviations2D Two-dimensionalBM Bistritzer-MacDonaldBZ Brillouin zoneCB Conduction bandCS-DFT Current- and spin-density functional theoryDoS Density of statesGL Ginzburg-LandauhBN Hexagonal boron nitrideIVC Inter-valley coherentMATBG Magic-angle twisted bilayer grapheneMBZ Moire´ Brillouin zoneQAH Quantum anomalous HallQAHE Quantum anomalous Hall effectSOC Spin-orbit couplingTBG Twisted bilayer grapheneTMD Transition metal dichalcogenideTMG Twisted multilayer grapheneVB Valence bandVP Valley polarizedxiAcknowledgementsI am grateful to my supervisor, Professor Marcel Franz, for his mentorship and support,especially through difficult periods of the COVID-19 pandemic. Many thanks for keepingme busy with interesting and diverse problems, and for listening to and encouraging myideas. I fondly remember our ski day discussions, and the many laps on Spanky’s Ladder.A special thanks to Dr. Benjamin T. Zhou, a great physicist and mentor whose ideashave formed the foundation of this work. I also thank all members of the Franz group,past and present, for welcoming me so warmly to UBC.xiiI dedicate this work first to my parents and sister, who still don’t understand what I dobut support me regardless.To my brother, whose company and cooking made Vancouver feel like home.Last but not least, to my partner Thomas, whose constant support helped me stay focusedand gain confidence throughout this process.xiiiChapter 1IntroductionWhen 2-dimensional (2D) crystals are stacked at a small twist angle, a moire´ patternforms possessing a different set of symmetries from the aligned bilayer. At certain com-mensurate twist angles, an emergent translational invariance can arise so that the twistedheterostructure forms a moire´ lattice, with a much larger unit cell of ∼ 10-100 atomsand correspondingly smaller moire´ Brillouin zone (MBZ). The discovery of flat bands andinteraction-driven states like superconductivity and correlated insulators in twisted bi-layer graphene (TBG) [1, 2] brought these materials to the forefront of condensed matterphysics. One of the exotic phenomena these systems can realize is the quantum anomalousHall effect (QAHE), in which a quantized Hall conductivity σxy is observed in the absenceof magnetic field1.In order for σxy to be quantized at zero field, the system must have bands characterizedby an integer-valued topological invariant called the Chern number, which is defined asthe integral over the Brillouin zone (BZ) of the Berry curvature Ωnk:Cn =∫BZd2k2πΩnk , (1.1)where n is the band index. If any of the bands are partially filled, i.e. the system ismetallic, then the Hall conductivity is given by the sum of Ωnk over occupied states:σxy =e2h∑n∫d2k2πnF (εnk)Ωnk , (1.2)where nF is the Fermi-Dirac distribution for a given chemical potential µ; and εnk is theenergy of the state at wave-vector k in band n. In the metallic state σxy can take any realvalue. In an insulating system, where µ is tuned within a band gap, Equation 1.2 reducesto the sum over Chern numbersσxy =e2h∑nCn ≡ e2hC (1.3)and σxy is guaranteed to be an integer multiple ofe2h , thus realizing the QAHE wheneverthe total Chern number C is non-zero. We often characterize this “Chern insulator” stateby its Hall resistivity Rxy =he21C instead of conductivity, since the former can be measureddirectly in experiment. Note that C vanishes in the presence of time-reversal symmetry(T ) [4], thus the presence of QAHE necessarily implies that T is broken intrinsically bythe system.1See David Tong’s Lectures on the Quantum Hall Effect for an excellent introduction to the basic theoryand phenomenology of Quantum Hall Effect [3]1Chapter 1. IntroductionThe discovery of a correlated quantum anomalous Hall (QAH) state in TBG alignedto hexagonal boron nitride (hBN) [5] was an exciting proof of concept; demonstrating thata moire´ heterostructure with flat bands could enter a T -broken state driven by Coulombinteractions, yielding a quantized Hall resistivity Rxy at zero magnetic field. Prior to this,QAH effect had only been observed in topological insulators decorated with magneticdopants whose ordered magnetic moments break T symmetry [6–8]. More recently, QAHhas been reported in untwisted AB stacked bilayer graphene [9], AB-stacked MoTe2 /WSe2 heterobilayer [10], and magic-angle TBG (MATBG) without hBN alignment [11].The experimental signature of the QAH state is a hysteretic response of the Hallresistivity as a function of out-of-plane magnetic field. The hysteresis effect is characterizedby an initial Rxy = 0 measured at zero field, which then increases as the field is turnedon up to a plateau at he2C , where C is the total Chern number, and persists as the fieldis tuned back to zero. If the field is then reversed, the Hall resistivity changes sign andapproaches − he2C , which again holds as the field returns to zero. The quantization of Rxyis a unique feature of the Chern insulating state.The properties of the hysteresis curve suggest that:1. the T -breaking which leads to non-zero Hall conductance is a spontaneous symmetrybreaking, since Rxy = 0 until a magnetic field is applied and Rxy can be polarizedin either direction;2. The electronic states which give rise to QAHE also exhibit some form of ferromag-netism, enabling magnetic switching of the Hall resistivity.The mechanism for ferromagnetism is presumably related to large orbital magnetizationof the topological bands, estimated at ∼ 1 − 10µB per moire´ unit cell in magic-angleTBG and twisted multilayer graphene (TMG) [12, 13]. Direct experimental evidence fororbital magnetization in TBG is mounting [14–16], making TBG the first example of aferromagnet driven mostly by the angular momentum of electrons, and not by their spin.When finite orbital magnetization arises concurrently with non-zero Chern number andQAHE, the state can be dubbed an “orbital Chern insulator.”In this work, we predict that in twisted homobilayers of the transition metal dichalco-genide (TMD) MoS2, the moire´ bands arising from monolayer conduction band (CB) statescan realize an orbital Chern insulator with C = ±2 and a large orbital magnetization of∼1µB per moire´ unit cell. Motivated by experimental observation of magnetic/electricalswitching of the Hall resistivity in twisted graphene heterostructures [5,17], we seek a the-oretical framework to describe the orbital magnetization’s response to an applied magneticfield. This leads us to derive a new formula for multi-band orbital susceptibility appli-cable to continuum models of moire´ materials, which depends heavily on the geometricproperties of individual Bloch bands and on interband coupling.In Chapter 2 we present our proposal for topological flat bands generated by spin-orbitcoupling (SOC) in twisted homobilayer MoS2 CB states, and demonstrate via Hartree-Fock analysis that the interacting ground state at 1/2-filling is a spontaneously T -brokenstate with non-zero Chern number C = ±2. In Chapter 3, we introduce the existing theoryof orbital magnetism and derive a numerically tractable susceptibility formula for multi-band continuum model Hamiltonians, expressed in terms of quantities reflecting the band2Chapter 1. Introductiongeometry such as Berry curvature and orbital magnetic moment. We will then use thesederivations in Chapter 4.3 to calculate the orbital magnetization of the moire´ bands, anddemonstrate the properties of the Fermi sea contribution to orbital susceptibility. Thisanalysis leads to interesting conclusions about the possible mechanisms for small-fieldmagnetic switching in moire´ quantum anomalous Hall devices.3Chapter 2Correlated quantum anomalousHall states in twisted homobilayerMoS2Hexagonal materials whose low energy states lie close to the Brillouin zone (BZ) corners(denoted K), such as graphene and transition metal dichalcogenides (TMDs), have pro-vided a convenient paradigm for theoretical study of moire´ physics. The monolayer bandstructure can be approximated by a continuum model close to each of the ±K points,which effectively act as a new degree of freedom called “valley.” When two monolayers areassembled at a small angular twist, the continuum model in each layer must be modified bya periodic moire´ potential, due to displacement of the atoms in the opposite layer, and in-terlayer tunneling. These are the basic ingredients of the celebrated Bistritzer-Macdonald(BM) model [18] used to describe twisted bilayer graphene (TBG). Many works haveshown that the moire´ patterns in TBG lead to flat bands at certain “magic” angles, whereinteractions can play a significant role in determining the system’s ground state [19–22].In this Chapter we establish a novel topological phase generated by SOC in CB moire´bands in twisted TMDs, which is essentially different from those in twisted graphene whereSOC is negligible, and those in the valence band (VB) moire´ bands of twisted TMDs whereSOC plays no essential role in creating nonzero Chern numbers [23]. Focusing on thespecific case of a twisted homobilayer MoS2 we predict that an interplay among angulartwist, local symmetry breaking and spin-orbit coupling in the CB states creates moire´ flatbands with nonzero Chern numbers C = ±2, and density-density interactions genericallydrive the system at 1/2-filling into a valley-polarized correlated QAH state.2.1 Continuum model for twisted homobilayer MoS2Stacking two monolayers of MoS2 with a small relative twist angle θ0 results in a moire´superlattice with lattice constant aM = a/θ0, where a is the lattice constant of eachmonolayer. The system has a symmetry group of D3⊗{Uv(1), T }, where Uv(1) stands forvalley conservation, T for time-reversal symmetry, and the point group D3 is generatedby the three-fold rotation about the z-axis (C3z) and a two-fold rotation about the y-axis(C2y) (Fig.2.1(a)). Uv(1) is an emergent symmetry due to the moire´ potential, whicheffectively decouples states near opposite K points of the untwisted BZ. This will allowus to describe states from opposite valleys with two independent Hamiltonians.For an aligned bilayer at θ0 = 0, the mirror symmetry Mz about the horizontal2D mirror plane, which is respected in a monolayer, is broken locally on each layer. This42.1. Continuum model for twisted homobilayer MoS2Figure 2.1: Schematic crystal structure of twisted bilayer MoS2 with (a) top view and (b)side view. Electrons in different layers experience opposite effective electric fields due tolocal lack of mirror symmetry. (c) Moire´ Brillouin zone (MBZ) formed by twisting hexag-onal Brillouin zones of each layer. The in-plane spinor field exhibits opposite vorticitiesat Km,+ and Km,− of MBZ due to layer-dependent Rashba effects.generates uniform electric fields of opposite signs in different layers, which remain effectiveupon a small angular twist with θ0 ∼ 1◦ (Fig.2.1(b)). Based on models of monolayer MoS2with broken Mz [24–26], the effective Hamiltonian for CB minima at ±K valleys of twotwisted layers, without moire´ potential or interlayer coupling, can be written as:Hξ0,l =∑k,αβc†k,l,αhξl,αβ(k)ck,l,β, (2.1)hξl (k) =ℏ2(k − ξKl)22m∗− µ+ ξβsoσz+ (−1)lαso[(k − ξKl)× σ] · zˆ,where ξ = ± denotes the valley index, l = 1(2) labels the top (bottom) layer, α, βlabel the spin indices with the σ-matrices representing the usual Pauli matrices for spins.K1 ≡Km,+ andK2 ≡Km,− are the shifted +K-points on layer 1 and layer 2 (Fig.2.1(c))due to the angular twist of θ0 equivalent to rotating the top (bottom) layer by an angleθ0/2 (−θ0/2) about the z-axis. In Eq.(2.1), m∗ ≈ 0.5me is the CB effective mass at K, µis the Fermi energy measured from band minima. The αso-term and βso-term account forthe Rashba SOC and the Ising SOC, respectively. The Ising SOC has a valley-dependentsign due to time-reversal symmetry and causes a splitting of ∆so = 2|βso| in the spin-upand spin-down CB states. The Rashba SOC has a layer-dependent sign imposed by C2ywhich swaps the two layers.Spatial modulations due to the formation of moire´ pattern can be captured by in-troducing coupling terms between states at k and k + GM,j [22, 23], where GM,j =− 4π√3aM(cos (j−1)π3 , sin(j−1)π3 ) are the reciprocal vectors of the moire´ superlattice. Themomentum-space moire´ potential is given byHM =∑k6∑j=1∑l=1,2∑α=↑,↓c†k,l,αVl,GM,jck+GM,j ,l,α, (2.2)52.2. Moire´ flat Chern bands and skyrmionic spin texturesTable 2.1: Parameters used in H0 for twisted homobilayer MoS2 with monolayer latticeconstant a = 3.16A˚.µ αso βso VM ψ w00 meV 80 meV·A˚ -1.5 meV 10 meV -89.6◦ 10 meVwhere Vl,GM,j = VM (cosψ−i(−1)jl sinψ) are complex coupling parameters with amplitudeVM and phase ψ.The effective inter-layer tunneling for states at K-valleys is modeled following thegeneral recipe of two-center approximation [18,23], which can be written asHξT = −w0∑kc†k,1(ck,2 + ck+ξGM,2,2 + ck+ξGM,3,2) + h.c., (2.3)where w0 denotes the effective inter-layer tunneling amplitude near K, which is extrap-olated to be w0 ≈ 10 meV using a combined approach of empirical models and Slater-Koster methods (see Appendix A for details). HM and HξT are modelled by constructinga 4N × 4N matrix, where the degrees of freedom are 2 layers, 2 spins, and N coupledk-points in the original BZ. We use a hexagonal mesh of k points, where the hexagon isdefined by its corners k+RGM,j , with R an integer parameter. We find good convergenceof the band structure and Berry curvature for R ∼ 5-6. Note that the matrix elementsof HM and HξT are constant in k, so that all of the explicit k-dependence is contained inhξl (k).The total non-interacting continuum model thus reads: H0 =∑ξ=±∑l=1,2Hξ0,l +HξT +HM with parameters listed in Table 2.1.2.2 Moire´ flat Chern bands and skyrmionic spin texturesFrom general considerations, level spacings among the lowest moire´ bands correspondroughly to the quantization energy ∆E ∼ ℏ2π2/(2m∗a2M ) of an electron confined in asuperlattice cell with size a2M = a2/θ20. For θ0 ∼ 1◦ − 2◦, ∆E ∼ 1 − 10 meV in twistedMoS2, which is comparable to the spin-orbit splitting ∆so ≡ 2|βso| = 3 meV caused bythe Ising SOC [27]. If spin is conserved (e.g., in the absence of Rashba SOC), spin-up andspin-down moire´ bands are expected to cross in general.By setting αso = 0 in H0, we find that for ξ = +, the 1st spin-down and the 2nd spin-upmoire´ bands cross each other for θ0 ∈ (1.25◦, 1.5◦). As a specific example, spin-resolvedmoire´ bands at θ0 = 1.4◦ with αso = 0 are shown in Fig.2.2(a). With spin conservationin this case, the spin Chern number is well-defined, and the spin-up and spin-down bandsinvolved in the level crossing have Chern numbers (C1,↓ = +1, C2,↑ = −1), which originatefrom the layer-pseudospin magnetic fields as in moire´ bands from VB states [23]. Uponturning on αso, crossing points are gapped out by the layer-dependent Rashba SOC whichmix the spin-up and spin-down species (Fig.2.2(b)). As a result of this mixing, the originalC1,↓ and C2,↑ from layer pseudospins annihilate each other; meanwhile, a new pair offlat bands (2nd and 3rd bands in Fig.2.2(b)) with Chern numbers C = ±2 emerge. The62.2. Moire´ flat Chern bands and skyrmionic spin texturesFigure 2.2: Moire´ bands of valley ξ = + formed by CB states of twisted homobilayer MoS2at θ0 = 1.4◦ with (a) αso = 0 and (b) αso ̸= 0. The colorbars indicate the out-of-plane spinexpectation value ⟨Sz⟩ in units of ℏ/2. (c) Momentum-space profile of Berry curvature Ωzin units of A˚2 in the 2nd moire´ band in (b). Large Ωz of order 104A˚2 emerges as a resultof the nontrivial gap induced by SOC. (d) Momentum-space spin texture of the 2nd moire´band in (b). White arrows indicate the in-plane spinor field, which has opposite vorticitiesat Km,+ and Km,− and leads to a nontrivial skyrmion number S = +2. Parameters usedare presented in Table 2.1.72.3. Correlated QAH state at 1/2 fillingnontrivial gap induced by Rashba SOCs is further signified by giant Berry curvatures oforder 104A˚2 in the MBZ (Fig.2.2(c)).The unusual Chern numbers C = ±2 arise from a unique combination of angular twistand local mirror-symmetry breaking: as moire´ bands at opposite Km,+ and Km,− pointsoriginate from states in different layers (Fig.2.1(c)), the layer-dependent Rashba SOCs(Eq.(2.1)) create an in-plane spinor field with opposite vorticities at Km,+ and Km,− inthe MBZ. This unsual pattern causes the in-plane spin to wind twice as one goes around aloop enclosing the Γm-point in the MBZ (Fig.2.2(d)). As we confirm numerically, the spintextures of the 2nd and 3rd moire´ bands are characterized by nonzero skyrmion numbersS = ∫ d2knˆk · (∂kx nˆk × ∂ky nˆk)/4π = ±2 (nˆk: spin orientation of state at k), and thesetwo bands describe a two-level system with one-to-one correspondence between Chernnumbers C and skyrmion number S [4].2.3 Correlated QAH state at 1/2 fillingGiven the narrow bandwidth W ≈ 1 meV of moire´ Chern bands (Fig.2.2(b)), the char-acteristic Coulomb energy scale gC ≃ e2/ϵaM ≈ 10 meV overwhelms the kinetic energyand correlation physics is expected to arise (assuming θ0 = 1.4◦ and dielectric constantϵ = 10 [28]). Motivated by the observation that twisted graphene at 3/4−filling can behaveas an SU(4)-quantum Hall ferromagnet [5,11,22,29], and the fact that all CB moire´ bandsin Fig.2.2(b) involves only a two-fold valley degeneracy (note that spin SU(2)-symmetryis already broken by SOC), we examine correlation effects when the 2nd moire´ band is1/2-filled. For simplicity, we drop the band index n in the following.Including density-density interactions, the low-energy effective Hamiltonian is: Heff =H0,eff + HI,eff , where H0,eff =∑k,ξ=±Eξ(k)c†ξ(k)cξ(k) is the non-interacting part withband energies Eξ(k), and the interacting Hamiltonian compatible with Uv(1) and T readsHI,eff =∑kk′q,ξξ′V (q)AΛξ(k + q,k)Λξ′(k′ − q,k′) (2.4)× c†ξ(k + q)c†ξ′(k′ − q)cξ′(k′)cξ(k),where V (q) is the Fourier transform of the two-body interaction, A is the total areaof the system, and Λξ(k + q,k) ≡ ⟨uξ,k+q|uξ,k⟩ (|uξ,k⟩: periodic part of the Blochstates) are the form factors arising from projecting the density-density interaction toactive bands. The general trial Hartree-Fock ground state at 1/2-filling has the form:|Φ⟩ =∏k∈MBZ[w+(k)c†+(k) +w−(k)c†−(k)] |Φ0⟩, where |Φ0⟩ is the vaccuum state with alllower-lying bands filled, and wξ(k) are the variational parameters.Minimizing EΦ ≡ ⟨Φ|Heff |Φ⟩ with mean-field approximations for Fock exchange terms(see Appendix B.1) leads to a set of self-consistent equations:∆+ +∆− = J0,∆+ −∆− = J02N∑k(∆+−∆−)−δ(k)D(k) ,∆+− = J12N∑k∆+−D(k) .(2.5)82.3. Correlated QAH state at 1/2 fillingFigure 2.3: (a) J0−J1 phase diagram of the system with the 2nd moire´ band at 1/2-filling.Under repulsive density-density interactions, only the J1 ≤ J0 regime can be accessed (seeAppendix B.4), where the valley-polarized(VP) state is generically favored and the systembecomes a correlated QAH insulator. (b) Evolutions of the order parameters ∆+− and∆M ≡ ∆+ −∆− as the ratio J0/J1 is tuned along the line-cut (green dashed line in (a))across the phase boundary. A first-order transition happens at the critical point J0 ≃ J1.Here, J0 and J1 denote the mean-field intra-valley and inter-valley exchange couplingconstants, ∆ξ = J0n¯ξ (n¯ξ: mean occupancy for valley ξ) and ∆+− are the intra-valleyand inter-valley order parameters, which characterize the average energy gains from intra-valley and inter-valley exchange interactions, respectively. δ(k) ≡ E+(k)−E−(k) denotesenergy difference between bands from two valleys ξ = ±,D(k) =√∆2+− + [δ(k)−∆M ]2/4,where we introduce the valley magnetic order ∆M ≡ ∆+ −∆− = J0(n¯+ − n¯−).Solutions of Eq.(2.5) are divided into two classes:(i) the valley-polarized (VP) state [22,30], in which one out of the two nearly degeneratemoire´ bands from two valleys is fully filled, with spontaneous T -symmetry breaking:∆M = ξJ0,∆+− = 0, wξ(k) = 1;(ii) the inter-valley coherent (IVC) state, in which the Slater determinant is formedby coherent superpositions of states from the two valleys, with spontaneous Uv(1)-symmetry breaking [22,30,31]: ∆M = 0,∆+− ≈ J1/2, w+,−(k) ̸= 0.By comparing the energies of the VP and IVC states, we obtain the mean-field J0 − J1phase diagram for the range 0.3gC ≤ J0, J1 ≤ gC (Fig.2.3(a)). The IVC state is favoredthroughout the entire J1 ≥ J0 regime, while the VP phase takes up most of the J1 < J0regime. By tuning J0/J1 across the phase boundary (green dashed line in Fig.2.3(a)), afirst-order transition occurs at J0 ≃ J1 (Fig.2.3(b)).It is worth noting that given a general repulsive interaction V (q) > 0, the inequality2|Λ+(k′,k)||Λ−(k,k′)| ≤ |Λ+(k′,k)|2 + |Λ−(k,k′)|2 together with T -symmetry leads toJ0 ≥ J1, regardless of details of the form factors and the microscopic interaction (Appendix92.3. Correlated QAH state at 1/2 fillingB.3). Thus, electron correlations most likely drive the system at 1/2-filling into a T -brokenVP state, which realizes a correlated QAH state with C = ±2.Furthermore, the topological insulating gap to be manifested experimentally at 1/2-filling is not solely determined by the gap induced by Rashba SOCs in the non-interactingbands. Instead, it is largely governed by the correlation-induced insulating gap for thevalley-polarized state, which is approximately given by the intra-valley exchange couplingstrength J0 ∼ gC ≈ 10 meVs.Physically, for a VP state, the intra-valley exchange interactions among electrons inthe fully filled moire´ band lower the energy of each individual electron in this valley by anamount of J0. On the other hand, the energy of the empty band in the other valley remainsunaffected, as electron correlation cannot occur within an empty band. As a result, oncethe non-interacting bands acquire nonzero Chern numbers from any finite Rashba SOC,the correlated topological gap at 1/2-filling would be on the order of J0 = 10 meV. Toexplicitly demonstrate this, we use our Hartree-Fock formalism (see Eq.B.19) to calculatethe correlated band structures for both valley ξ = + and valley ξ = − when the pair of2nd moire´ bands are half-filled (Fig.2.4). As the spontaneous T -broken VP phase selectsone out of the two valleys randomly, we assume without loss of generality that the 2ndmoire´ band from ξ = + is fully filled, while its counterpart from ξ = − is empty. Clearly,the topological insulating gap (indicated by double arrow in Fig.2.4(c)-(d)) correspondsto the scale of J0, instead of the gap induced by Rashba SOC in the non-interacting bandsin Fig.2.4(a)-(b).102.3. Correlated QAH state at 1/2 fillingFigure 2.4: (a)-(b) Non-interacting bands of twisted homobilayer MoS2 for (a) valleyξ = +, and (b) valley ξ = − with Rashba strength αso = 30 meV· A˚. Dashed lines in (a)-(b) indicate location of the Fermi energy EF when the 2nd moire´ bands are 1/2-filled andinteraction effects are neglected. In the absence of correlations, both bands from valleyξ = + and valley ξ = − are half-filled. (c)-(d) Correlated band structures when the 2ndmoire´ bands are half-filled, where the system enters the valley-polarized state. Electroncorrelations induce an insulating gap ∆CQAH ∼ J0, where we set J0 = 10 meVs.11Chapter 3Theory of orbital magnetismIn the previous Chapter, we demonstrated that electron-electon interactions induce atransition to a valley polarized (VP) state with broken time reversal symmetry (T ) inour model of twisted homobilayer MoS2. This situation is analogous to magnetic ex-change interactions driving long range ordering in a ferromagnet, and indeed spontaneousT -breaking is usually associated with ferromagnetism. In both cases, interactions areessential for inducing a phase transition and stabilizing the symmetry-broken state.There is now strong experimental evidence for orbital ferromagnetism in TBG, suchas hysteretic response of the Hall resistivity Rxy under applied magnetic field (both quan-tized [5] and unquantized [29]), electrical switching of the magnetization [5,17], and evendirect measurements of the magnetization with scanning tunneling microscopy [14] andsuperconducting quantum interference devices [16]. These measurements suggest that anon-zero orbital magnetization manifests in tandem with non-zero Chern number, enablingmagnetic control of the quantum Hall state.While orbital magnetization in TBG and other moire´ systems has been extensivelystudied theoretically [12, 13, 32], the magnetization response to external fields (i.e. sus-ceptibility) is not well understood. The goal of this chapter is to derive a formula fororbital susceptibility in coupled multi-band models, and discuss its potential applicationsin describing the magnetic properties of moire´ systems. In Section 3.1 we review the so-called “modern theory of orbital magnetization,” while Section 3.2 introduces the theoryof orbital susceptibility. In Section 3.2.2 we use methods from the literature to derive anew formula for multi-band orbital susceptibility applicable to continuum models of moire´bands.3.1 Orbital magnetization of crystalline systemsElectrons in crystalline solids can possess a form of magnetism distinct from their intrinsicspin. The so-called “orbital magnetic moment” is analogous to a magnetic dipole momentgenerated by the motion of a classical charged particle, and interacts with external fieldsin similar ways. However, a problem arises when we try to compute matrix elements ofthe magnetic dipole operator er×p in Bloch bands: the operator r is ill-defined for Blochstates, which are delocalized in real space. We can circumvent this issue by repackagingthe Bloch states |ψnk⟩ into wave packets which are sharply peaked in both k and r, sothat they have a well-defined centre-of-mass, denoted kc and rc, in each space. Unlikeclassical point particles, these wave packets have a finite spread in real-space, and thuscan exhibit a self-rotation about rc which generates the orbital magnetic moment [33]mnk = −i e2ℏ⟨∂kunk| × (hk − εnk) |∂kunk⟩ , (3.1)123.1. Orbital magnetization of crystalline systemswhere ∂k is a gradient with respect to crystal momentum k, |unk⟩ is the cell-periodicpart of the Bloch state, hk is the k-dependent Hamiltonian, and εnk the correspondingeigenvalues at zero field. However, the self-rotation does not tell the whole story. Whilebulk wave packets carry no net current, and thus contribute angular momentum onlythrough self-rotation, wave packets near the boundary of a finite sample may have a non-zero centre-of-mass velocity parallel to the edge [33,34]. This edge current contributes tothe total magnetic moment like a classical current loop ∼ I ·A, with A the area of the loop,yet its effect can be expressed in terms of the bulk energy spectrum only. Thus orbitalmagnetization represents yet another case, similar to the bulk-boundary correspondenceof topological insulators [4], where surface properties are inextricably tied to properties ofthe bulk band structure.The total orbital magnetization at finite temperature and zero field can be expressedas an integral over the BZ [33]M(µ, T ) =∑n∫BZdk(2π)d[nF (εnk)mnk +eβℏln(1 + e−β(εnk−µ))Ωnk], (3.2)where nF (εnk) =(1 + eβ(εnk−µ))−1is the Fermi-Dirac distribution at β = 1/kBT , mnk isthe orbital magnetic moment defined in Eq.(3.1), and Ωnk is the Berry curvatureΩnk = i ⟨∂kunk| × |∂kunk⟩ . (3.3)At T = 0, we can show that Eq.(3.2) reduces toM(µ, T = 0) =∑n∫ µBZdk4π2[mnk +eℏ(µ− εnk)Ωnk], (3.4)where the upper limit indicates integration over energy states below the chemical potentialεnk < µ. We will denote the latter contribution to magnetic momentmΩnk =eℏ(µ−εnk)Ωnkand the total mtotnk =mnk +mΩnk.Eq.(3.4) was first derived by Xiao et al. based on the semiclassical theory of electrondynamics [35], in which electromagnetic fields are treated as a classical perturbation tothe band energy and eigenstates of a quantum Hamiltonian. They show that the semi-classical equations of motion in the presence of magnetic field and Berry curvature leadto a modified density of states (DoS), such that all k-space integrals acquire an additionalfactor ∫dk(2π)d→∫dk(2π)d(1 +eℏB ·Ωnk). (3.5)Including this DoS correction and the first-order Zeeman correction to the band energyε˜nk = εnk −mnk ·B, the thermodynamic grand potential density is given byF (µ, T,B) = − 1β∑n∫dk4π2(1 +eℏB ·Ωnk)ln(1 + e−β(εnk−mnk·B−µ)), (3.6)from which we can easily obtain Eq.(3.2) by applying the definitionM(µ, T ) = − ∂F (µ, T,B)∂B∣∣∣∣B=0. (3.7)133.2. Orbital susceptibilityThe above formulae reveal the close connection between orbital magnetization, a mea-surable observable, and the electronic band topology characterized by the Berry curvature.The doping-dependence of the orbital magnetization can also serve as a signature of non-zero Chern number. In an insulating state, we might expect a constant magnetizationwith respect to changes in µ, yet Eq.(3.4) tells us that:M(µ, T = 0)|µ∈gap ∼ C · µ+ const. , (3.8)where C =∑n Cn is the sum over the Chern number of each band n, given by the integralof the Berry curvature over the full Brillouin zone Cn ≡∫BZdk2πΩnk. Thus when totalChern number is non-zero, increasing the chemical potential within a gap leads to a linearscaling of the magnetization with slope ∝ C.The modern theory of orbital magnetization presented here has been well-establishedfor over a decade, with various approaches all converging to the zero-temperature mag-netization formula of Eq.(3.4). This includes an approach based on maximally-localizedWannier functions [34, 36]; a fully quantum derivation [37]; and a Green’s function for-malism which unifies approaches to both polarization and magnetization [38, 39]. Thelatter generalizes the DoS correction described above, and will prove useful for studyingelectric/magnetic response functions including orbital susceptibility.3.2 Orbital susceptibilityCalculating the zero-field magnetization from Eq.(3.2) informs us on how the orbital mag-netic moments of Bloch states add up as a band is filled, but does not tell us how thosemoments will respond when a magnetic field is turned on. This response is known as themagnetic susceptibility, denoted χ. We define the zero-field orbital susceptibility in termsof the orbital magnetization and/or grand potential as follows:χorb,ij = µ0∂Mi(µ, T,B)∂Bj∣∣∣∣B=0= −µ0 ∂2F (µ, T,B)∂Bi∂Bj∣∣∣∣B=0, (3.9)where µ0 is the vacuum permeability. In this work we will only consider 2D systems in atransverse magnetic field, and thus define χorb ≡ χorb,zz for simplicity.The main difficulty of deriving an orbital susceptibility formula consists in calculatingcorrections to the grand potential up to second order in the magnetic field. This leads toa proliferation of terms which makes the final results unwieldy and difficult to interpretphysically (see e.g. Refs. [40, 41]).Recently, approaches which complement the modern theory of orbital magnetizationand emphasize the geometric properties of orbital susceptibility have been developed [39,42, 43]. Gao et al. [43] employ the semi-classical approach developed in e.g. Ref. [44] tocalculate second-order corrections to the wave packet energy, and thereby the free energy.For our own exploration, we will derive the general formula for orbital susceptibility of tightbinding models introduced by Raoux et al. [39], using the gauge-invariant Green’s functionformalism for electrons in uniform magnetic fields described by Chen and Lee [38]. Fromthis starting point, we arrive at a numerical formula for multi-band orbital susceptibility,143.2. Orbital susceptibilityapplicable to BM-type continuum models of moire´ systems. The resulting formula is anew result with the following key features:1. Many contributions can be recast in terms of familiar quantities which reflect thegeometric properties of Bloch bands; including Berry curvature, quantum metrictensor, and orbital magnetic moment.2. It is transparent how to implement the formula numerically for a given band model,once the eigenstates and band energies are known.3. All terms can be expressed in terms of matrix multiplications for efficient computa-tion.3.2.1 Gauge-invariant Green’s function formalismIn this section, we will summarize Raoux et al.’s [39] derivation for the first and secondorder corrections to the Green’s function of a tight binding model in a uniform magneticfield. To begin, we re-write the grand potential asF (µ, T,B) = − 1βA∫ ∞−∞dω ln(1 + e−β(ω−µ))ρ(ω,B) , (3.10)where A is the area of the 2D system, and ρ(ω,B) is the DoS, which is modified inthe presence of the magnetic field. We can obtain Mz and χorb from the definitions inEqs.(3.7) and (3.9), respectively. The DoS can be expressed in terms of the retardedGreen’s function G(ω,B) of the real-space Hamiltonian in the magnetic field Hρ(ω,B) = − 1πImTrG(ω,B) , (3.11)whereG(ω,B) = (ω −H)−1 (3.12)and we will obtain a perturbative expansion for G in the field. Note that G is a matrix inthe space of lattice sites.To incorporate the effect of the magnetic field, we perform a Peierls substitution onthe zero-field Hamiltonian in real space:h =∑i,jtij |i⟩ ⟨j| (3.13)where i and j are site indices. The Peierls substitution adds a phase to each matrixelement tij → tijeiφij , withφij =eℏ∫ xjxiA(r) · dℓ . (3.14)A(r) is the vector potential for the static and uniform magnetic field, so that φij clearlyis gauge dependent. However, we can add up the phases around a triangle of lattice sites153.2. Orbital susceptibilityi,j,k to define a gauge-invariant magnetic flux, using Stokes theorem and B =∇×A(r):Φijk = φij + φjk + φki =eℏ∮ijkA(r) · dℓ = eℏ∫∫ijkB · dS , (3.15)where dℓ is the line element of the path around the triangle, and dS is the triangle’ssurface element. The resulting flux is given in Cartesian coordinates by:Φijk =eB2ℏ[(xi − xk)(yk − yj)− (yi − yk)(xk − xj)] . (3.16)Now, let g(ω) be the the Green’s function of h, and G(ω,B) that of the Hamiltonian inthe magnetic field:H =∑i,jtijeiφij |i⟩ ⟨j| . (3.17)We will define a “twisted” Green’s functiong˜ij ≡ eiφijgij , (3.18)and expand G in terms of g˜, before re-introducing g to obtain a gauge-invariant expressionfor the diagonal elements Gii. We can show that(ω −H)g˜ = 1−X , (3.19)whereXij = eiφij∑k(eiΦikj − 1)tikgkj . (3.20)Combining Eq.(3.19) with the definition of G, we can writeG = (ω −H)−1 = g˜(1−X )−1 = g˜∑n≥0X n ≡∑n≥0G(n) . (3.21)The DoS only depends on diagonal terms of G, which turn out to be gauge-invariant. Thezeroth term is G(0)ii = g˜ii = gii, as expected. By expanding X in powers of Φikj , we canshow thatG(1)ii = i∑jkΦjkigijtjkgki (3.22)andG(2)ii = −12∑jkΦ2jkigijtjkgki −∑jklmΦjkiΦlmigijtjkgkltlmgmi . (3.23)Defining position operators x =∑i xi |i⟩ ⟨i| and y =∑i yi |i⟩ ⟨i|, we can express G(1)ii andG(2)ii in terms of commutators of the position with h and g. We will use a shorthandnotation for commutators of position with an operator O, Ox ≡ − iℏ [x,O]. Using thesenotations and the expression for the flux in Eq.(3.16), the corrections to the Green’sfunction becomeG(1)ii = −ieB2ℏ(ghxghyg − ghyghxg)ii (3.24)163.2. Orbital susceptibilityandG(2)ii = −e2B28ℏ2[g(hxxgyy + hyygyy − 2hxygxy + 2(hxghyg − hyghxg)2)]ii, (3.25)where we have also used the relation for the Green’s function commutators gx = − iℏ [x, g] =ghxg (likewise for y). The hx terms have a physical interpretation as the velocity operatorat zero field vx ≡ x˙ = − iℏ [x, h].Eqs.(3.24,3.25) can be reinterpreted for a k-dependent Hamiltonian in momentumspace. The eigenfunctions of the Hamiltonian h are the Bloch states |ψnk(r)⟩ = eik·r |unk(r)⟩with eigenvalue εnk, where |unk(r)⟩ is a periodic function with the same periodicity asthe crystal. |unk(r)⟩ is then an eigenfunction of hk = e−ik·rheik·r, where r is the positionoperator, with the same eigenvalue εnk. Meanwhile, the single particle Green’s functioncan be expressed in the basis of Bloch states:g(ω) =A4π2∑m∫BZd2k|ψmk⟩ ⟨ψmk|ω − εmk + iη , (3.26)where η > 0 is a small parameter required to ensure convergence of the Fourier transform.Thus G(1)ii will be composed of products of expectation values ⟨ψnk|hi |ψmk⟩ (i = x, y).Now, consider the following expectation values of derivatives with respect to hk:⟨unk| ∂xhk |umk⟩ = ⟨unk| (−ixhk + ihkx) |umk⟩ = ⟨unk| e−ik·rhxeik·r |umk⟩ , (3.27)where we have used the fact that x and r commute. Therefore if we let g be the singleparticle Green’s function in momentum spaceg(k, ω) =∑m|umk⟩ ⟨umk|ω − εmk + iη (3.28)and include the integral A∫BZd2k4π2in the trace over G, then the hi appearing in Eqs.(3.24-3.25) can be interpreted as the ki-derivatives of hk. This hk could be, for example, aBM-type continuum model such as the one described in Section 2.1.The orbital magnetization can then be obtained from the integralMz(µ, T ) =ImπA∫ ∞−∞dωln(1 + e−β(ω−µ))βTr{ie2ℏ(ghxghyg − ghyghxg)}, (3.29)by substituting in Eq.(3.28), and performing the energy integration (see Appendix A ofRef. [39]). A direct comparison to Eq.(3.2) requires us to use the gauge-invariant formulaefor orbital magnetic moment and Berry curvature, given bymznk =ie2ℏ∑m ̸=n⟨unk|hxk |umk⟩ ⟨umk|hyk |unk⟩ − ⟨unk|hyk |umk⟩ ⟨umk|hxk |unk⟩(εnk − εmk) (3.30)andΩznk = i∑m̸=n⟨unk|hxk |umk⟩ ⟨umk|hyk |unk⟩ − ⟨unk|hyk |umk⟩ ⟨umk|hxk |unk⟩(εnk − εmk)2 , (3.31)173.2. Orbital susceptibilityrespectively (see Appendix C for derivation).Finally, after expressing G(2)ii as a derivative with respect to ω and performing anintegration by parts, we can replace all commutators with g and arrive at a general formulafor orbital susceptibility:χorb(µ, T ) = −µ0e212ℏ2ImπA∫ ∞−∞dω nF (ω) Tr{ghxxghyy − ghxyghxy− 4 (ghxghxghyghy − ghxghyghxghy)} , (3.32)where we have now dropped the k subscripts so that hi represents the derivative with re-spect to ki of the momentum-space Hamiltonian, g represents g(k, ω) defined in Eq.(3.28),and the Tr operator includes the integral over the BZ.3.2.2 Multi-band formula for orbital susceptibilityThe formula for orbital susceptibility derived in Section 3.2.1 is valid for any tight bindingmodel with a finite number of bands, and in principle can be applied directly once thederivatives of hk and its band structure are known. However the physical interpretations ofeach term are not transparent in the current form, and the energy integral is complicatedto implement numerically. We thus seek to establish a numerically tractable form ofEq.(3.32) that we can apply to continuum models of moire´ systems.First, we can integrate by parts to get an alternative form of Eq.(3.32) [39]:χorb(µ, T ) = −µ0e22ℏ2ImπA∫ ∞−∞dω nF (ω) Tr{ghxghyghxghy+12(ghxghy + ghyghx)ghxy} . (3.33)This form has been used previously in Refs. [45, 46] to study the orbital susceptibility ofmonolayer graphene and MoS2. The continuum model we wish to study depends on k onlythrough the parabolic dispersion assumed in each monolayer and is therefore separable,i.e. hxy = 0. What remains is identical to the famous Fukuyama formula [47]:χorb(µ, T ) = −µ0e22ℏ2ImπA∫ ∞−∞dω nF (ω) Tr{ghxghyghxghy︸ ︷︷ ︸≡Xˆ} . (3.34)Our task now is to express Tr Xˆ in the band basis, and perform the integral over ωanalytically so that the orbital susceptibility can be expressed directly in terms of theband structure. Using Eq.(3.28) and defining the trace over bands at a single k pointTr′ Xˆ =∑n ⟨unk| Xˆ |unk⟩, we have:Tr′ Xˆ =∑n,m,l,p⟨n|hx |m⟩ ⟨m|hy |l⟩ ⟨l|hx |p⟩ ⟨p|hy |n⟩(ω − εn)(ω − εm)(ω − εl)(ω − εp) (3.35)where we have suppressed the η parameters, dropped the k subscript on the eigenvaluesεn, and simplified notation for the eigenstates |unk⟩ → |n⟩. Each index is summed overall bands.183.2. Orbital susceptibilityTo compute the energy integral, we will need to perform a partial fraction decompo-sition of Tr′ Xˆ in order to make use of the following identity (see Appendix D):Im∫ ∞−∞dωf(ω)(ω − εn + iη)s = −π(s− 1)!f(s−1)(εn) (3.36)where s is an integer, f(ω) is any real function, and f (s−1)(ω) is its (s − 1)th derivative.These integrals will produce various higher derivatives of nF (ω) in the final susceptibilityformula. We leave all details of the derivation to Appendix D, and simply state the finalresult here:χorb(µ, T ) =µ0e2ℏ2∑n∫BZdk(2π)d{n′′′F (εn)(εxnεyn)212+ n′′F (εn)ℏeεxnεynReMnk− n′F (εn)[ℏ2e2(m2nk − (ReMnk)2)+ 2εxnεyngnxyk−∑m,lεxn ⟨n|hy |m⟩ ⟨m|hx |l⟩ ⟨l|hy |n⟩+ (x↔ y)(εn − εm)(εn − εl)]+ nF (εn)[ℏe(mnkΩnk − 2gnxykReMnk) + εxnεyn∑m⟨n|hx |m⟩ ⟨m|hy |n⟩+ c.c.(εn − εm)3−∑m,lεxn ⟨n|hy |m⟩ ⟨m|hx |l⟩ ⟨l|hy |n⟩+ (x↔ y) + c.c.(εn − εm)2(εn − εl)+∑m,l,p⟨n|hx |m⟩ ⟨m|hy |l⟩ ⟨l|hx |p⟩ ⟨p|hy |n⟩+ c.c.(εn − εm)(εn − εl)(εn − εp)]}.(3.37)Eq.(3.37) contains several new quantities, which we will catalogue in order of appearance.In the first line, the εin are ki-derivatives of the band energyεin ≡ ⟨n|hi |n⟩ = ∂iεn , (3.38)and Mnk is a complex tensor defined asMnk =eℏ∑m ̸=n⟨n|hx |m⟩ ⟨m|hy |n⟩εn − εm (3.39)of which the imaginary part is related to the out-of-plane orbital magnetic momentmnk = − eℏ Im∑m ̸=n⟨n|hx |m⟩ ⟨m|hy |n⟩εn − εm= −ImMnk . (3.40)We defined Mnk in analogy with the quantum geometric tensor introduced by BerryTnij(k) =∑m ̸=n⟨n|hi |m⟩ ⟨m|hj |n⟩(εn − εm)2 . (3.41)193.2. Orbital susceptibilityIts real part is the quantum metric tensor that appears in the second and third line ofEq.(3.37)gnxyk = ReTnxy(k), (3.42)and its imaginary part is related to the out-of-plane Berry curvatureΩnk = −2Im∑m ̸=n⟨n|hx |m⟩ ⟨m|hy |n⟩(εn − εm)2= −2ImTnij . (3.43)The third term in the second line can be re-written in terms of quantities of a single bandεxnεyn∑m⟨n|hx |m⟩ ⟨m|hy |n⟩+ c.c.(εn − εm)3 = 2εxnεynRe ⟨∂xn| (εn − hk)−1 |∂yn⟩ , (3.44)where we have used the relation ⟨m|hi|n⟩εn−εm = ⟨m|∂in⟩ for m ̸= n and (εn − hk)−1 |m⟩ =(εn − εm)−1 |m⟩. This leaves only three inseparable “interband” terms. We isolate theircontribution to the susceptibility as:χinter =µ0e2ℏ2∑n∫BZdk(2π)d{n′F (εn)∑m,lεxn ⟨n|hy |m⟩ ⟨m|hx |l⟩ ⟨l|hy |n⟩+ (x↔ y)(εn − εm)(εn − εl)+ nF (εn)[ ∑m,l,p⟨n|hx |m⟩ ⟨m|hy |l⟩ ⟨l|hx |p⟩ ⟨p|hy |n⟩+ c.c.(εn − εm)(εn − εl)(εn − εp)−∑m,lεxn ⟨n|hy |m⟩ ⟨m|hx |l⟩ ⟨l|hy |n⟩+ (x↔ y) + c.c.(εn − εm)2(εn − εl)]}. (3.45)These terms arise from “vertical mixing” of Bloch states at the same k-point but in differ-ent bands, which appears as a correction to the wave function in semiclassical perturbationtheory [43].We also define a contribution that was previously identified in two-band models [42],and has been studied in the 2D kagome lattice [48]:χΩ = µ0∑n∫BZdk(2π)d{− n′F (εn)m2nk + nF (εn)eℏmnkΩnk}. (3.46)The remaining terms all involve derivatives of the energy, as we can show by taking aderivative of Eq.(3.38) that ∂2xyεn ≡ εxyn = ℏeReMnk when hxy = 0. We collect them inone final termχε =µ0e2ℏ2∑n∫BZdk(2π)d{n′′′F (εn)(εxnεyn)212+ n′′F (εn)εxnεynεxyn − n′F (εn)[(εxyn )2 + 2εxnεyngnxyk]+ nF (εn)[2εxnεynRe ⟨∂xn| (εn − hk)−1 |∂yn⟩ − 2εxyn gnxyk]}.(3.47)203.2. Orbital susceptibilityBy integrating the first two terms by parts, we should be able to recover the “band byband” extension to the Peierl’s formula for orbital susceptibility in a single-band model [49]χLP ≡ µ0e2ℏ2∑n∫BZdk(2π)dn′F (εn)(εxxn εyyn − (εxyn )2), (3.48)which was defined by Raoux et al. [39] to represent the intra-band contributions to orbitalsusceptibility. Other terms which contribute to the geometric orbital susceptibility termsχg and χ˜g (defined in Ref. [42]) should also appear from the integration by parts.Thus the total susceptibility is given in terms of Eqs.(3.45-3.47) as:χorb = χΩ + χε + χinter . (3.49)χorb consists of both Fermi surface terms proportional to derivatives of nF (εn); andFermi sea terms proportional to nF (εn), which produce plateaus in χorb when µ is insidea band gap. As illustrated in Fig.2.4), the spontaneous T -breaking VP state which wouldinduce non-zero net Chern number is expected to realize a large correlation-induced gap.The Fermi sea terms are important for understanding magnetic response in an insulatingstate, as all other contributions go to zero. Yet these contributions have been overlookedin previous studies of orbital susceptibility in TBG [12,50].The central question is whether the above theory of orbital susceptibility can predictthe hysteretic response of the Hall resistivity Rxy observed in moire´ quantum Hall devices[5,11,15,17,29,51]. Broadly speaking, we can classify a material’s susceptibility as eitherparamagnetic (χ > 0), where the magnetic moments prefer to align with the appliedmagnetic field; or diamagnetic (χ < 0), where applying a field induces a magnetic momentwhich opposes the change in the field. These devices’ tendency to enter a ferromagneticstate after an external field is applied is consistent with paramagnetic susceptibility (χzz >0) in the Chern insulator; the magnetized state is then stabilized by electron-electroninteractions as the field is turned off, leading to spontaneous T -breaking. However, wecan show that at least one of the Fermi sea terms must be diamagnetic. Pie´chon et al. [42]demonstrated that χΩ individually obeys the susceptibility sum rule for tight bindingmodels [46] ∫ ∞−∞dµχorb(µ, T ) =∫ ∞−∞dµχΩ(µ, T ) = 0 , (3.50)which enforces that the susceptibility goes to zero when the chemical potential is far abovethe energy scale of the (finite) bands. The sum rule for χΩ implies that the sum of itsFermi sea terms mnkΩnk is always diamagnetic, since −n′F (εn) ≥ 0∀ εn, µ implies thatthe Fermi surface term is paramagnetic. Therefore, if in-gap paramagnetism is to appear,it must come from either χε or χinter.Although the formalism leading to Eq.(3.37) was developed for tight binding models,we argue that it is also valid for BM-type continuum models of moire´ materials. Liketight binding models, BM models couple a finite number of degrees of freedom (in mo-mentum space, rather than real space), and exhibit a discrete translational invariance onthe moire´ lattice scale. This means that the periodic Bloch wave functions and energies arewell-defined within a reduced moire´ Brillouin zone. Furthermore, Gutie´rrez-Rubio et al.213.2. Orbital susceptibilityshowed qualitative agreement between orbital susceptibility calculations using both tightbinding and effective continuum models for monolayer graphene and MoS2 [45]. Otherworks in moire´ systems have also shown agreement for energy dispersion and topologicalproperties between continuum and tight binding models [23].In Appendix D.3, we show that each contribution to the full susceptibility formulain Eq.(3.37) can be computed efficiently using matrix operations, even the more complexinterband terms involving sums over multiple band indices. This will allow us to calculatethe essential contribution of χinter to the orbital susceptibility in an insulating state.22Chapter 4Orbital magnetism in topologicalmoire´ bandsOrbital magnetization and its relationship with the QAH state has been studied exten-sively in the four-fold degenerate flat bands of TBG, where two valleys ξ = ±1 and twospins s =↑, ↓ generate a four-flavour manifold with approximate SU(4) symmetry [52,53].Electron-electron interactions may lift the degeneracy in valley, spin, or both; and the or-der in which these degeneracies are lifted will determine what state the system manifestsat different fillings. Multiple theoretical studies have suggested that different combina-tions of spin/valley symmetry breaking at the same 3/4 filling factor may lead to either aquantum anomalous Hall effect, or unquantized anomalous Hall conductance [12,21]; bothof which have been observed experimentally [5, 29].In the moire´ bands arising from MoS2 CB states, spin SU(2) symmetry is alreadybroken by SOC, and only a two-fold valley degeneracy remains. This makes it a simplermodel system for studying spontaneous T -breaking, and we may be able to infer generalprinciples about the magnetism of these two-valley topological flat band systems moreeasily. As we showed in Section 2.3, an interaction-induced valley polarized (VP) state withlarge correlated gap is likely to form whenever the 2nd moire´ band is 1/2-filled. Thus weexpect that the chemical potential can generally be tuned within a topological correlatedgap in order to yield quantized Hall conductance in tandem with orbital magnetism.In this Chapter, we will use the formulae for orbital magnetization (Mz) to showthat the flat moire´ bands of our twisted homobilayer MoS2 model each carry a largeout-of-plane orbital magnetization. We then model the magnetization’s interaction withan external field using a Ginzburg-Landau theory. Finally, we calculate the Fermi seacontribution to the susceptibility, and find paramagnetic susceptibility in a finite rangeof chemical potential. However, the magnitude of the susceptibility is not large enoughto explain magnetic switching of the magnetization/Hall resistivity induced by small (∼0.1T) magnetic fields. We discuss the implications of this on the nature of the orbitalsusceptibility in moire´ quantum Hall devices.4.1 Large out-of-plane magnetization in twistedhomobilayer MoS2 moire´ bandsRecent studies suggest that the flat bands in twisted graphene heterostructures carry verylarge orbital magnetic moment, with theoretical estimates of the orbital magnetization onthe order of 1-10µB per moire´ unit cell [12, 14, 32]. Similarly large orbital magnetizationcan be expected in flat bands of the twisted homobilayer MoS2 model presented in Chapter234.1. Large out-of-plane magnetization in twisted MoS2Figure 4.1: Distribution of out-of-plane mtotnk = mnk +mΩnk in the MBZ for the ξ = +1moire´ flat bands labelled in Fig.2.2(b). The chemical potential µ is set to the top of eachband. While the exact colour scale varies for each band, zero is always represented by apale yellow, and the range is set to [−cmaxn , cmaxn ], where cmaxn = maxk∈MBZ∣∣mtotnk∣∣.244.1. Large out-of-plane magnetization in twisted MoS2Figure 4.2: Comparison of bulk mnk and edge mΩnk magnetic moment contributions in (a)C = +2 and (b) C = −2 bands of valley ξ = +1. The chemical potential µ is set to the topof each band, and the colour scale is defined in the same way as Fig.4.1. M bulkz (Medgez ) isdefined as the average of mznk (mz,Ωnk ) over the MBZ.254.1. Large out-of-plane magnetization in twisted MoS22, given their large Berry curvature as shown in Fig.2.2(c) and the correspondence betweenΩnk and mtotnk highlighted in Section 3.1. We re-copy the definition of zero-temperatureorbital magnetization here for the reader’s convenience:M(µ, T = 0) =∑n∫ µBZdk4π2[mnk +eℏ(µ− εnk)Ωnk︸ ︷︷ ︸≡mΩnk](4.1)Using the continuum model defined in Section 2.1, we numerically evaluate the total out-of-plane orbital magnetic moment mz,totnk = mznk + mz,Ωnk for the non-interacting moire´bands plotted in Fig.2.2(b). Orbital magnetization is calculated numerically as the k-space average of the total magnetic moment Mz =1N∑kmz,totnk , and reported in unitsof µB per moire´ unit cell. We use a hexagonal grid of k-points within the MBZ withN = 3169. Note that we will usually drop the z superscript, so that any scalar form ofm, Ω refers to the out-of-plane component.We plot the distribution of mtotnk in the MBZ for each of the four flat bands in valleyξ = +1 in Fig.4.1 . All four Chern bands carry a large orbital magnetization, ∼ 1µBper moire´ unit cell. Bands in the opposite valley ξ = −1 are related by time reversalsymmetry, and therefore have equal and opposite orbital magnetization. Note that bandswith equal and opposite Chern number within the same valley do not necessarily haveequal magnitude equal and opposite orbital magnetization; nor is the latter quantized.Rather, the magnitude and sign of Mz depends on details of the band dispersion. Fig.4.2shows a comparison of mnk and mΩnk for the C = ±2 bands. We see that the edgecomponent mΩnk contributes significantly to both the distribution of mtotnk in the MBZand to the overall magnitude of the orbital magnetization, demonstrating that the bulkcomponent mnk alone is insufficient to describe the overall orbital magnetization. Thediscrepancy between bulk and total magnetization was famously demonstrated by Ceresoliet al. using real space simulations [36], leading them to propose an empirical formula forthe “itinerant” (edge) contribution.The edge contribution’s effect is even more pronounced when the chemical potentialis tuned within a band gap, as shown in Fig.4.3. M bulkz ≡ 1N∑kmnk plateaus as allavailable states in the band become filled, while M edgez ≡ 1N∑kmΩnk continues to increaselinearly in µ beyond the band edge. Note that the slope of Mz above the band edge isproportional to the Chern number C of the filled band, as anticipated in Section 3.1. Thelinear scaling of Mz vs. µ is present even at finite temperature (provided the gap is largeenough to resolve it), and is considered a signature of currents carried by topological edgestates in the gap.The results presented in Figs.4.1-4.3 represent isolated contributions from non-interactingbands in a single valley (ξ = +1). We emphasize that the net magnetization of the non-interacting system is always zero due to the overall T -symmetry of the Hamiltonian H0;thus the magnetization of bands in opposite valleys exactly cancel at any filling. Interac-tions, or an external field, must lift the degeneracy between bands in opposite valleys inorder for a non-zero Mz to arise. On the other hand, the large orbital and edge magneticmoment carried by each band mean that any degree of valley polarization would lead toan observable orbital magnetization on the order of ∼ 0.1-1µB/u.c. Thus we can think of264.1. Large out-of-plane magnetization in twisted MoS2Mz as an order parameter for the phase transition to a VP state. In the next section wewill construct a Ginzburg-Landau free energy for Mz in order to predict the coercive fieldrequired to switch the magnetization of the VP state.26.2 26.0 25.8 25.6 25.4 25.2µ (meV)0.40.20.00.20.40.60.81.0Mz (µB / u.c.)Mbulkz +MedgezMbulkzM edgez(a) 1st band, C = −125 24 23 22 21µ (meV)0.00.51.01.52.02.53.0Mz (µB / u.c.)Mbulkz +MedgezMbulkzM edgez(b) 2nd band, C = 224.0 23.5 23.0 22.5 22.0 21.5 21.0 20.5 20.0µ (meV)3.02.52.01.51.00.50.00.5Mz (µB / u.c.)Mbulkz +MedgezMbulkzM edgez(c) 3rd band, C = −222.0 21.5 21.0 20.5 20.0 19.5 19.0 18.5µ (meV)0.00.20.40.60.81.01.21.41.6Mz (µB / u.c.)Mbulkz +MedgezMbulkzM edgez(d) 4th band, C = 1Figure 4.3: Mz vs. µ for each moire´ Chern band (ξ = +1) identified in Fig2.2(b) atT = 0.2K. M bulkz and Medgez are defined as in Fig.4.2. The magnetization plotted hererepresents the contribution from the states of a single band, not the total magnetizationof all filled states. The magnetization contributed by the bulk term plateaus when µ isabove the band edge (grey shaded region), while the edge contribution continues to scalelinearly in proportion to Chern number.274.2. Ginzburg-Landau theory for magnetic switching4.2 Ginzburg-Landau theory for magnetic switchingBased on the analysis in Section 2.3, we predict that interactions drive the twisted MoS2system into a VP state when the 2nd moire´ band is half-filled. As this is a spontaneoussymmetry breaking effect and neither valley is preferred a priori, we expect magneticdomains with opposite magnetization and Chern number to form, each with their owntopological edge currents. The Hall resistivity Rxy observed in transport measurementsmay not be quantized, or even finite, until the sample is polarized by a magnetic field.This picture is consistent with the hysteretic Rxy measured in twisted graphene devicesat 3/4 filling [5, 17, 29], and with recent imaging of the orbital magnetism in TBG [16].The degeneracy of the two valleys also implies that the magnetization can be switchedwhen an opposing magnetic field is applied. We study a Ginzburg-Landau (GL) theoryto describe this phenomenon below, inspired by Ref. [12].We model the competition between valleys with the following GL free energy, usingnet orbital magnetization Mz (in µB per unit cell) as an order parameter:f = −a0M2z + b0M4z −MzBz . (4.2)The units are [f ] = eV, [Bz] = T, [a0] = µ−2B eV, and [b0] = µ−4B eV. The free energy isminimized by the equation:∂f∂Mz= 0 = −2a0Mz + 4b0M3z −Bz . (4.3)When the system is below the critical temperature for magnetic ordering a0 and b0 arepositive, such that at Bz = 0 the free energy has two degenerate minima with oppositemagnetization:±Mz,0 = ±√a02b0. (4.4)Mz,0 represents the magnitude of the magnetization when the system is in a spontaneousT -broken VP state. If we take the derivative of Equation 4.3 with respect to Bz, we canderive an expression for zero-field susceptibility in terms of the parameters:0 = −2a0∂Mz∂Bz+ 12b0M2z∂Mz∂Bz− 1 (4.5)χzz,0µ0≡ ∂Mz∂Bz∣∣∣∣Bz=0=1−2a0 + 12b0Mz,0 =14a0(4.6)whereχzz,0µ0has units µB/T, and represents the zero-field susceptibility in the VP state.Rearranging Equation 4.4 and 4.6, we obtain expressions for the GL theory parameters interms of Mz,0 and χ˜zz,0 ≡ χzz,0µ0 :a0 =14χ˜zz,0(4.7)b0 =18M±2z χ˜zz,0. (4.8)284.2. Ginzburg-Landau theory for magnetic switchingOnce a0 and b0 are known the GL theory is fully determined, and the only free parameteris the external field Bz.Our goal now is to determine the external field required to switch the sign of themagnetization once the system is in a VP state. In other words, sitting at zero field andmaximum |Rxy|, what field must be applied to “tilt” the system into the opposite valley?Mathematically, this is equivalent to finding the value of Bz at which Equation 4.3 hasonly a single real root. We see that the polynomial can be rearranged into the form of adepressed cubic t3 + pt+ q, for which the discriminant is∆ = −(4p3 + 27q2) . (4.9)The depressed cubic equation has three real roots for ∆ > 0, and only one real root for∆ < 0. In our case p = − a02b0 , q = Bz4b0 and ∆ is:∆ =a302b30− 27 B2z16b20(4.10)Setting ∆ = 0 and solving for Bz, we obtain the coercive field:Bc =4a03√a06b0=Mz,03√3χ˜zz,0(4.11)To estimate a0, b0, and Bc for our twisted MoS2 continuum model, we can independentlycalculate Mz,0 and χzz,0 from the theory developed in Chapter 3. However, note that a0and therefore χzz,0 must be positive (paramagnetic) in order for well-defined free energyminima to exist. This is by no means guaranteed, especially when the chemical potentialsits inside the correlated gap of the VP state (illustrated in Fig.2.4) where only Fermi seaterms should contribute; and depends greatly on the band details. We already showed inSection 3.2 that the Fermi sea contribution to χΩ (Eq.(3.46)) is necessarily diamagnetic.This puzzle motivates us to numerically evaluate the orbital susceptibility arising fromFermi sea terms, which would determine the magnetic response of the system in the VPstate if the chemical potential is tuned inside the correlated gap.294.3. Orbital susceptibility of the Fermi sea4.3 Orbital susceptibility of the Fermi seaLet us first re-state the Fermi sea terms of χorb:χsea =µ0e2ℏ2∑n∫BZdk(2π)dnF (εn){ℏemnkΩnk− 2ℏegnxykReMnk + εxnεyn∑m⟨n|hx |m⟩ ⟨m|hy |n⟩+ c.c.(εn − εm)3−∑m,lεxn ⟨n|hy |m⟩ ⟨m|hx |l⟩ ⟨l|hy |n⟩+ (x↔ y) + c.c.(εn − εm)2(εn − εl)+∑m,l,p⟨n|hx |m⟩ ⟨m|hy |l⟩ ⟨l|hx |p⟩ ⟨p|hy |n⟩+ c.c.(εn − εm)(εn − εl)(εn − εp)}≡∑n(χn,Ω + χn,ε + χn,inter) . (4.12)Although the VP state is clearly a phenomenon of interacting bands, we continue to use thenon-interacting band structure for our calculation; as the effect of interactions on orbitalsusceptibility is beyond the scope of this work. Our main goal is to determine whethera paramagnetic in-gap susceptibility, necessarily driven by the Fermi sea contribution ofχinter and/or χε , is possible in the topological moire´ bands.To this end, we plot χsea for the valley ξ = +1 as a function of chemical potentialµ in Fig.4.4, and obtain a somewhat surprising result. The states in the 2nd moire´ band(C = 2), which we showed enters a VP state at 1/2-filling, indeed yield a large paramag-netic contribution to χsea, with χ2,inter and χ2,Ω (defined in Eq.(4.12)) both positive. Atlow temperature, their combined contribution is large enough to overcome the diamag-netism of the lower band and of χ2,ε, generating an overall paramagnetic susceptibilitynear the 2nd band maximum. Fig.4.4(b) simply exaggerates the gap between the 2nd and3rd(C = ±2) moire´ bands, which do not cross but have an indirect gap of −0.08meV in thenon-interacting band structure, to show that a paramagnetic plateau would arise if thisgap were somehow enhanced. By increasing the temperature to T = 1.8K, we find thatthe paramagnetic region is highly sensitive to temperature and gap size; effectively disap-pearing when the gap is determined by the Rashba SOC alone. At the same temperaturebut larger gap, paramagnetism persists but starts to decay at a lower chemical potential,due to filling of diamagnetic states in the C = −2 band. Thus the in-gap paramagnetismbetween the 2nd and 3rd moire´ bands is an interesting exception to the general rule ofdiamagnetic Fermi sea susceptibility. The paramagnetism is driven largely by the χintercontribution, which arises from band coupling induced by the applied field.Interestingly, this result is consistent with an observation by Vignale that the interbandorbital susceptibility of a 2D electron gas in a periodic potential produces paramagnetismat a narrow indirect band gap [41]. The only other feature known to consistently produceparamagnetism is a Van Hove singularity, where the band energy has a saddle point, inwhich case the paramagnetism is driven by the intra-band Landau-Peierls susceptibilitygiven in Eq.(3.48) [39,41,42,45,46].304.3. Orbital susceptibility of the Fermi seaWe note that χ4,Ω and χ4,inter are also positive, although the interband contribution ismuch smaller than that of χ2,inter. This suggests a general rule that whenMbulkz andMedgezhave the same sign, χn,Ω and χn,inter are paramagnetic. The reason for this is obvious inthe case of χn,Ω, which is given by the integral of mn,kΩnk, but is less transparent in thecase of χn,inter.While it is not clear that we would obtain similar results for the orbital susceptibilityin the interacting bands, it is instructive to use the order of magnitude of χsea to estimatethe corresponding coercive field. Assuming a positive in-gap orbital susceptibility χseaµ0 ∼+0.1µB/T and Mz ∼ 1µB in the VP state, we obtain Bc ∼ 10T. This coercive fieldis extremely large compared to the ∼ 10−1T fields with which magnetic switching ofquantized anomalous Hall resistivity has been achieved in TBG heterostructures [5, 54]and MoTe2/WSe2 heterobilayers [10]. Such a large coercive field precludes realization ofmagnetic switching in the spontaneous T -broken moire´ Chern insulator described here,for a very fundamental reason. In small-angle Moire´ heterostructures with lattice constant∼ 10 nm, such as the MoS2 homobilayer studied here with aM = 12.9 nm, a field of ∼ 10Tis large enough that the magnetic flux through one unit cell approaches the flux quantumΦ0 =h2e ∼ 103Tnm2. At this scale of magnetic flux, the combination of periodic potentialand magnetic field splits the electronic states into Landau levels (LL) with magneticsubbands, leading to a fractal energy spectrum known as the Hofstadter butterfly [55,56].Indeed, this long-anticipated phenomenon was observed for the first time experimentallyin a moire´ heterostructure of Bernal-stacked bilayer graphene on hBN substrate, with alattice constant of 15.5±0.9 nm and magnetic fields ∼ 10T [57]. This system also realizesthe integer quantum Hall effect, but the quantization of the Hall conductivity is relatedto the ratio of carrier density to magnetic flux, rather than the intrinsic topology of themoire´ bands which produces the quantum anomalous Hall effect at zero field.The relatively small magnitude of in-gap susceptibility, contributed entirely by χsea,combined with its propensity for diamagnetism suggest that the Fermi sea susceptibilityalone cannot be responsible for magnetic switching of the anomalous Hall resistivity in flatmoire´ bands. However if the chemical potential is tuned close to a band edge, χΩ certainlycontributes a significant paramagnetic susceptibility, which may be further enhanced byχinter. Pie´chon et al. explored two examples of flat bands in two-band systems, andfound in both cases that paramagnetic susceptibility arose due to dominant interbandcontributions [42]. In the case of Mielke’s checkerboard model, the flat band caused adivergent paramagnetic susceptibility. Both Fermi surface and Fermi sea contributionsof χε are suppressed due to the non-dispersive nature of the flat bands, and will haveminimal influence on the overall sign of the susceptibility. From these considerations,we hypothesize that orbital magnetic switching originating in a VP state in moire´ Chernbands is driven by Fermi surface interband orbital susceptibility, as well as χΩ, where thechemical potential must be tuned close to the interacting band edge.314.3. Orbital susceptibility of the Fermi sea26 25 24 23 22 21 20 19µ (meV)0.060.040.020.000.020.040.06χsea/µ0 (µB/T)χΩχinterχεχorb(a) T = 0.2K, no gap26 24 22 20 18µ (meV)0.060.040.020.000.020.040.06χsea/µ0 (µB/T)χΩχinterχεχorb(b) T = 0.2K, 1meVgap26 25 24 23 22 21 20 19µ (meV)0.060.040.020.000.02χsea/µ0 (µB/T)χΩχinterχεχorb(c) T = 1.8K, no gap26 24 22 20 18µ (meV)0.060.040.020.000.020.040.06χsea/µ0 (µB/T)χΩχinterχεχorb(d) T = 1.8K, 1meV gapFigure 4.4: Fermi sea orbital susceptibility χsea vs. µ for non-interacting moire´ bands ofvalley ξ = +1 at different temperatures. Each point incorporates the contribution of allstates which are filled at a given µ, by summing over each k ∈ MBZ and band n. Greyshaded regions indicate the width of each band, with the slight overlap between the 2ndand 3rd band (C = ±2) appearing in darker grey in panel (a) and (c). (a)-(b) and (c)-(d)contain the same data, but we artificially increase the energy of the 3rd and 4th band by1 meV in to highlight the paramagnetic susceptibility at the top of the 2nd band. Thewidth and strength of the paramagnetic region is highly sensitive to temperature and gapsize.32Chapter 5Discussion and outlookWhile ongoing activities in twisted TMDs have focused mainly on VB moire´ bands, ourproposal of moire´ Chern bands generated by SOC opens a new pathway into the largelyunknown territory of CB moire´ physics. When these Chern bands are half-filled, accordingto our predictions in Section 2.3, electron interactions can lead to a spontaneously T -broken VP phase and realize a correlated QAH state, which is signaled by quantized Hallresistance in transport measurements. The Hall resistivity Rxy should be accompaniedby a large orbital magnetization of ∼ 1µB per moire´ unit cell carried by the topologicalflat bands, as demonstrated in Section 4.1. The nontrivial orbital magnetism of the moire´bands provides a mechanism for magnetic switching of the Hall resistivity, as describedby the Ginzburg-Landau theory in Section 4.2.The main original result of this thesis is the multi-band orbital susceptibility formulafor separable Hamiltonians like the BM-type continuum models used to describe moire´systems, developed in Section 3.2.2. We use this formula to study the zero-field suscepti-bility of our twisted homobilayer MoS2 model, focusing on the Fermi sea terms given inEq.(4.12). These terms are easier to compute numerically, and we expected them to be thedominant contribution in the VP state with quantized Hall resistivity, where the chem-ical potential potential is tuned within a correlation-induced gap between Chern bandsin opposite valleys (as illustrated in Fig.2.4). However, we found only a narrow region ofchemical potential at the edge of the 2nd moire band in which χsea yields the paramag-netic response necessary for magnetic switching of Mz and Rxy; and this region shrinksquickly with increasing temperature or inverse gap size. Furthermore, the magnitude ofthe paramagnetic susceptibility implies a coercive field of 10 T is necessary to switch themagnetization, which is so large that Landau level physics would set in.We therefore predict that magnetic/electrical switching in moire´ quantum Hall devicesmust be driven by Fermi surface contributions to orbital susceptibility, especially thoseof χinter and χΩ. Computing these contributions is the natural next step in our study oforbital susceptibility in moire´ flat bands.To our knowledge, only two other groups have studied the orbital susceptibility of moire´systems quantitatively [12,50]. He et al. [12] consider only the Fermi surface contributionof χΩ in order to get the desired paramagnetic response, without including any of thecompeting effects that are likely to produce diamagnetism such as the Landau-Peierls term.The formula used by Guerci et al. [50] appears complete, but they make few connectionsto the phenomenon of magnetic switching and in-gap susceptibility.Although our treatment of the susceptibility is an improvement on these previousapproaches, it is still incomplete. The most significant omission in our calculation ofboth orbital magnetization and susceptibility is that of interactions, which we showed inSection 2.3 play an important role in determining the ground state of twisted MoS2. Shi et33Chapter 5. Discussion and outlookal. showed that the generalization to interacting systems for orbital magnetization in theframework of current- and spin-density functional theory (CS-DFT) is surprisingly simple:the formula for orbital magnetization is identical, but we substitute the interacting bandenergies obtained from CS-DFT [37]. There is little existing literature on the effect ofinteractions on orbital susceptibility, except in a few specific models [58, 59]. We couldalso incorporate effect of spin on the susceptibility, which must be modified from theconventional Pauli susceptibility in systems with SOC [60]. However we do not expectspin susceptibility to play a large role in the C = ±2 moire´ bands, as the integrated spinmagnetic moment is about an order of magnitude smaller than the orbital magnetization.We hope that the theory of orbital susceptibility explored here will be useful for pre-dicting the conditions under which magnetic switching of anomalous Hall resistance canbe realized in moire´ materials, and for identifying more candidate systems for robustmagnetic memory devices. For example, details in the µ dependence of orbital suscep-tibility may hold the key to understanding why QAHE and magnetic switching of Rxyare more robust at 3/4-filling than 1/4-filling of the four degenerate spin/valley flavoursin graphene-based moire´ heterostructures [17], even though both are predicted to havevalley-polarized ground states with non-zero Chern number [61].The theory may also be applicable to more exotic examples of magnetic switching,such as the magnetically tunable superconducting diode effect observed recently in magic-angle TBG Josephson junctions [62]. New methods for measuring orbital susceptibilitydirectly in monolayer materials [63] may help advance our understanding of this complexelectronic property.34Bibliography[1] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, andP. Jarillo-Herrero, “Unconventional superconductivity in magic-angle graphenesuperlattices,” Nature, vol. 556, no. 7699, pp. 43–50, Mar. 2018. [Online]. Available:https://doi.org/10.1038%2Fnature26160[2] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. ToMarken, J. Y. Luo, J. D.Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, R. C. Ashoori, andP. Jarillo-Herrero, “Correlated insulator behaviour at half-filling in magic-anglegraphene superlattices,” Nature, vol. 556, no. 7699, pp. 80–84, Mar. 2018. [Online].Available: https://doi.org/10.1038%2Fnature26154[3] D. Tong, “Lectures on the Quantum Hall Effect,” 2016. [Online]. Available:https://arxiv.org/abs/1606.06687[4] B. Bernevig and T. Hughes, Topological Insulators and Topological Superconductors.Princeton University Press, 2013. [Online]. Available: https://books.google.ca/books?id=wOn7JHSSxrsC[5] M. Serlin, C. L. Tschirhart, H. Polshyn, Y. Zhang, J. Zhu, K. Watanabe,T. Taniguchi, L. Balents, and A. F. Young, “Intrinsic quantized anomalous Halleffect in a moire´ heterostructure,” Science, vol. 367, no. 6480, p. 900–903, Feb. 2020.[Online]. Available: http://dx.doi.org/10.1126/science.aay5533[6] R. Yu, W. Zhang, H.-J. Zhang, S.-C. Zhang, X. Dai, and Z. Fang, “Quantizedanomalous Hall effect in magnetic topological insulators,” Science, vol. 329, no.5987, pp. 61–64, 2010. [Online]. Available: https://www.science.org/doi/abs/10.1126/science.1187485[7] C.-Z. Chang, J. Zhang, X. Feng, J. Shen, Z. Zhang, M. Guo, K. Li, Y. Ou,P. Wei, L.-L. Wang, Z.-Q. Ji, Y. Feng, S. Ji, X. Chen, J. Jia, X. Dai, Z. Fang,S.-C. Zhang, K. He, Y. Wang, L. Lu, X.-C. Ma, and Q.-K. Xue, “Experimentalobservation of the quantum anomalous Hall effect in a magnetic topologicalinsulator,” Science, vol. 340, no. 6129, pp. 167–170, 2013. [Online]. Available:https://www.science.org/doi/abs/10.1126/science.1234414[8] C.-Z. Chang, W. Zhao, D. Y. Kim, H. Zhang, B. A. Assaf, D. Heiman, S.-C.Zhang, C. Liu, M. H. W. Chan, and J. S. Moodera, “High-precision realization ofrobust quantum anomalous Hall state in a hard ferromagnetic topological insulator,”Nature Materials, vol. 14, no. 5, pp. 473–477, Mar. 2015. [Online]. Available:https://doi.org/10.1038%2Fnmat420435Bibliography[9] F. R. Geisenhof, F. Winterer, A. M. Seiler, J. Lenz, T. Xu, F. Zhang, and R. T.Weitz, “Quantum anomalous Hall octet driven by orbital magnetism in bilayergraphene,” Nature, vol. 598, no. 7879, p. 53–58, Oct. 2021. [Online]. Available:http://dx.doi.org/10.1038/s41586-021-03849-w[10] T. Li, S. Jiang, B. Shen, Y. Zhang, L. Li, Z. Tao, T. Devakul, K. Watanabe,T. Taniguchi, L. Fu, J. Shan, and K. F. Mak, “Quantum anomalous Hall effectfrom intertwined moire´ bands,” Nature, vol. 600, no. 7890, p. 641–646, Dec. 2021.[Online]. Available: http://dx.doi.org/10.1038/s41586-021-04171-1[11] K. P. Nuckolls, M. Oh, D. Wong, B. Lian, K. Watanabe, T. Taniguchi, B. A.Bernevig, and A. Yazdani, “Strongly correlated Chern insulators in magic-angletwisted bilayer graphene,” Nature, vol. 588, no. 7839, p. 610–615, Dec. 2020.[Online]. Available: http://dx.doi.org/10.1038/s41586-020-3028-8[12] W.-Y. He, D. Goldhaber-Gordon, and K. T. Law, “Giant orbital magnetoelectriceffect and current-induced magnetization switching in twisted bilayer graphene,”Nature Communications, vol. 11, no. 1, Apr. 2020. [Online]. Available:http://dx.doi.org/10.1038/s41467-020-15473-9[13] J. Liu, Z. Ma, J. Gao, and X. Dai, “Quantum valley Hall effect, orbitalmagnetism, and anomalous Hall effect in twisted multilayer graphene systems,”Physical Review X, vol. 9, no. 3, Aug. 2019. [Online]. Available: http://dx.doi.org/10.1103/PhysRevX.9.031021[14] S.-Y. Li, Y. Zhang, Y.-N. Ren, J. Liu, X. Dai, and L. He, “Experimental evidence fororbital magnetic moments generated by moire´-scale current loops in twisted bilayergraphene,” Physical Review B, vol. 102, p. 121406, Sep. 2020. [Online]. Available:https://link.aps.org/doi/10.1103/PhysRevB.102.121406[15] A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney, K. Watanabe,T. Taniguchi, M. A. Kastner, and D. Goldhaber-Gordon, “Evidence of orbitalferromagnetism in twisted bilayer graphene aligned to hexagonal boron nitride,”Nano Letters, vol. 21, no. 10, p. 4299–4304, May 2021. [Online]. Available:http://dx.doi.org/10.1021/acs.nanolett.1c00696[16] C. L. Tschirhart, M. Serlin, H. Polshyn, A. Shragai, Z. Xia, J. Zhu, Y. Zhang,K. Watanabe, T. Taniguchi, M. E. Huber, and A. F. Young, “Imaging orbitalferromagnetism in a moire´ Chern insulator,” Science, vol. 372, no. 6548, pp. 1323–1327, Jun. 2021. [Online]. Available: https://doi.org/10.1126%2Fscience.abd3190[17] H. Polshyn, J. Zhu, M. A. KuMar, Y. Zhang, F. Yang, C. L. Tschirhart, M. Serlin,K. Watanabe, T. Taniguchi, A. H. MacDonald, and A. F. Young, “Electrical switchingof magnetic order in an orbital Chern insulator,” Nature, vol. 588, no. 7836, pp.66–70, Nov. 2020. [Online]. Available: https://doi.org/10.1038%2Fs41586-020-2963-8[18] R. Bistritzer and A. H. MacDonald, “Moire´ bands in twisted double-layer graphene,”Proceedings of the National Academy of Sciences, vol. 108, no. 30, pp. 12 233–12 237,2011. [Online]. Available: https://www.pnas.org/doi/abs/10.1073/pnas.110817410836Bibliography[19] J. Liu, J. Liu, and X. Dai, “Pseudo Landau level representation of twistedbilayer graphene: Band topology and implications on the correlated insulatingphase,” Physical Review B, vol. 99, p. 155415, Apr. 2019. [Online]. Available:https://link.aps.org/doi/10.1103/PhysRevB.99.155415[20] Z. Song, Z. Wang, W. Shi, G. Li, C. Fang, and B. A. Bernevig,“All magic angles in twisted bilayer graphene are topological,” PhysicalReview Letters., vol. 123, p. 036401, Jul. 2019. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevLett.123.036401[21] N. Bultinck, S. Chatterjee, and M. P. Zaletel, “Mechanism for anomalousHall ferromagnetism in twisted bilayer graphene,” Physical Review Letters., vol.124, p. 166601, Apr. 2020. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevLett.124.166601[22] Y.-H. Zhang, D. Mao, Y. Cao, P. Jarillo-Herrero, and T. Senthil, “Nearly flat Chernbands in moire´ superlattices,” Physical Review B, vol. 99, p. 075127, Feb. 2019.[Online]. Available: https://link.aps.org/doi/10.1103/PhysRevB.99.075127[23] F. Wu, T. Lovorn, E. Tutuc, I. Martin, and A. H. MacDonald, “Topologicalinsulators in twisted transition metal dichalcogenide homobilayers,” PhysicalReview Letters., vol. 122, p. 086402, Feb. 2019. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevLett.122.086402[24] B. T. Zhou, K. Taguchi, Y. Kawaguchi, Y. Tanaka, and K. T. Law, “Spin-orbit coupling induced valley Hall effects in transition-metal dichalcogenides,”Communications Physics, vol. 2, no. 1, Mar. 2019. [Online]. Available:https://doi.org/10.1038%2Fs42005-019-0127-7[25] K. Taguchi, B. T. Zhou, Y. Kawaguchi, Y. Tanaka, and K. T. Law, “Valley Edelsteineffect in monolayer transition-metal dichalcogenides,” Physical Review B, vol. 98, p.035435, Jul. 2018. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevB.98.035435[26] K.-H. Kim and H.-W. Lee, “Berry curvature in monolayer MoS2 with broken mirrorsymmetry,” Physical Review B, vol. 97, p. 235423, Jun. 2018. [Online]. Available:https://link.aps.org/doi/10.1103/PhysRevB.97.235423[27] G.-B. Liu, W.-Y. Shan, Y. Yao, W. Yao, and D. Xiao, “Three-band tight-binding model for monolayers of group-VIB transition metal dichalcogenides,”Physical Review B, vol. 88, p. 085433, Aug. 2013. [Online]. Available:https://link.aps.org/doi/10.1103/PhysRevB.88.085433[28] Z. Bi and L. Fu, “Excitonic density wave and spin-valley superfluid in bilayertransition metal dichalcogenide,” Nature Communications, vol. 12, Jan. 2021.[Online]. Available: https://doi.org/10.1038/s41467-020-20802-z37Bibliography[29] A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney, K. Watanabe, T. Taniguchi, M. A.Kastner, and D. Goldhaber-Gordon, “Emergent ferromagnetism near three-quartersfilling in twisted bilayer graphene,” Science, vol. 365, no. 6453, pp. 605–608, 2019.[Online]. Available: https://www.science.org/doi/abs/10.1126/science.aaw3780[30] H. C. Po, L. Zou, A. Vishwanath, and T. Senthil, “Origin of Mott insulatingbehavior and superconductivity in twisted bilayer graphene,” Physical Review X,vol. 8, p. 031089, Sep. 2018. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevX.8.031089[31] Y.-H. Zhang, D. Mao, and T. Senthil, “Twisted bilayer graphene alignedwith hexagonal boron nitride: Anomalous Hall effect and a lattice model,”Physical Review Researchh, vol. 1, p. 033126, Nov. 2019. [Online]. Available:https://link.aps.org/doi/10.1103/PhysRevResearch.1.033126[32] O. Antebi, A. Stern, and E. Berg, “In-plane orbital magnetization as a probe forsymmetry breaking in strained twisted bilayer graphene,” Physical Review B, vol.105, no. 10, Mar. 2022. [Online]. Available: https://doi.org/10.1103%2Fphysrevb.105.104423[33] D. Xiao, M.-C. Chang, and Q. Niu, “Berry phase effects on electronic properties,”Reviews of Modern Physics, vol. 82, no. 3, p. 1959–2007, Jul. 2010. [Online].Available: http://dx.doi.org/10.1103/RevModPhys.82.1959[34] T. Thonhauser, D. Ceresoli, D. Vanderbilt, and R. Resta, “Orbital magnetizationin periodic insulators,” Physical Review Letters, vol. 95, no. 13, Sep. 2005. [Online].Available: http://dx.doi.org/10.1103/PhysRevLett.95.137205[35] D. Xiao, J. Shi, and Q. Niu, “Berry phase correction to electron density of statesin solids,” Physical Review Letters, vol. 95, no. 13, Sep. 2005. [Online]. Available:http://dx.doi.org/10.1103/PhysRevLett.95.137204[36] D. Ceresoli, T. Thonhauser, D. Vanderbilt, and R. Resta, “Orbital magnetizationin crystalline solids: Multi-band insulators, Chern insulators, and metals,”Physical Review B, vol. 74, no. 2, Jul. 2006. [Online]. Available: http://dx.doi.org/10.1103/PhysRevB.74.024408[37] J. Shi, G. Vignale, D. Xiao, and Q. Niu, “Quantum theory of orbital magnetizationand its generalization to interacting systems,” Physical Review Letters, vol. 99, no. 19,Nov. 2007. [Online]. Available: https://doi.org/10.1103%2Fphysrevlett.99.197202[38] K.-T. Chen and P. A. Lee, “Unified formalism for calculating polarization,magnetization, and more in a periodic insulator,” Physical Review B, vol. 84, p.205137, Nov. 2011. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevB.84.205137[39] A. Raoux, F. Pie´chon, J.-N. Fuchs, and G. Montambaux, “Orbital magnetism incoupled-bands models,” Physical Review B, vol. 91, p. 085120, Feb. 2015. [Online].Available: https://link.aps.org/doi/10.1103/PhysRevB.91.08512038Bibliography[40] E. I. Blount, “Bloch electrons in a magnetic field,” Phys. Rev., vol. 126, pp. 1636–1653,Jun. 1962. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRev.126.1636[41] G. Vignale, “Orbital paramagnetism of electrons in a two-dimensional lattice,”Physical Review Letters., vol. 67, pp. 358–361, Jul. 1991. [Online]. Available:https://link.aps.org/doi/10.1103/PhysRevLett.67.358[42] F. Pie´chon, A. Raoux, J.-N. Fuchs, and G. Montambaux, “Geometric orbitalsusceptibility: Quantum metric without Berry curvature,” Physical Review B,vol. 94, no. 13, Oct. 2016. [Online]. Available: http://dx.doi.org/10.1103/PhysRevB.94.134423[43] Y. Gao, S. A. Yang, and Q. Niu, “Geometrical effects in orbital magneticsusceptibility,” Physical Review B, vol. 91, p. 214405, Jun. 2015. [Online]. Available:https://link.aps.org/doi/10.1103/PhysRevB.91.214405[44] G. Sundaram and Q. Niu, “Wave-packet dynamics in slowly perturbedcrystals: Gradient corrections and Berry-phase effects,” Physical Review B,vol. 59, no. 23, pp. 14 915–14 925, Jun. 1999. [Online]. Available: https://doi.org/10.1103%2Fphysrevb.59.14915[45] A. Gutie´rrez-Rubio, T. Stauber, G. Go´mez-Santos, R. Asgari, and F. Guinea,“Orbital magnetic susceptibility of graphene and MoS2,” Physical Review B, vol. 93,no. 8, Feb. 2016. [Online]. Available: https://doi.org/10.1103%2Fphysrevb.93.085133[46] G. Go´mez-Santos and T. Stauber, “Measurable lattice effects on the charge andmagnetic response in graphene,” Physical Review Letters, vol. 106, no. 4, Jan. 2011.[Online]. Available: https://doi.org/10.1103%2Fphysrevlett.106.045504[47] H. Fukuyama, “Theory of orbital magnetism of Bloch electrons: Coulombinteractions,” Progress of Theoretical Physics, vol. 45, no. 3, pp. 704–729, 03 1971.[Online]. Available: https://doi.org/10.1143/PTP.45.704[48] Z. Wang and P. Zhang, “Orbital magnetization and its effects in spin-chiralferromagnetic kagome´ lattice,” Physical Review B, vol. 76, no. 6, Aug. 2007.[Online]. Available: https://doi.org/10.1103%2Fphysrevb.76.064406[49] R. Peierls, “Zur Theorie des Diamagnetismus von Leitungselektronen,” Zeitschrift furPhysik, vol. 80, no. 11-12, pp. 763–791, Nov. 1933.[50] D. Guerci, P. Simon, and C. Mora, “Moire´ lattice effects on the orbital magneticresponse of twisted bilayer graphene and condon instability,” Physical Review B, vol.103, no. 22, Jun. 2021. [Online]. Available: https://doi.org/10.1103%2Fphysrevb.103.224436[51] Y.-M. Xie, C.-P. Zhang, J.-X. Hu, K. F. Mak, and K. Law, “Valley-polarized quantum anomalous Hall state in moire´ MoTe2/WSe2 heterobilayers,”Physical Review Letters, vol. 128, no. 2, Jan. 2022. [Online]. Available:http://dx.doi.org/10.1103/PhysRevLett.128.02640239Bibliography[52] Y. Wang, J. Kang, and R. M. Fernandes, “Topological and nematic superconductivitymediated by ferro-SU(4) fluctuations in twisted bilayer graphene,” Physical ReviewB, vol. 103, no. 2, Jan. 2021. [Online]. Available: https://doi.org/10.1103%2Fphysrevb.103.024506[53] N. Bultinck, E. Khalaf, S. Liu, S. Chatterjee, A. Vishwanath, and M. P.Zaletel, “Ground state and hidden symmetry of magic-angle graphene at eveninteger filling,” Physical Review X, vol. 10, no. 3, Aug. 2020. [Online]. Available:https://doi.org/10.1103%2Fphysrevx.10.031034[54] C.-C. Tseng, X. Ma, Z. Liu, K. Watanabe, T. Taniguchi, J.-H. Chu, andM. Yankowitz, “Anomalous Hall effect at half filling in twisted bilayer graphene,”2022. [Online]. Available: https://arxiv.org/abs/2202.01734[55] D. R. Hofstadter, “Energy levels and wave functions of Bloch electrons in rationaland irrational magnetic fields,” Physical Review B, vol. 14, pp. 2239–2249, Sep.1976. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevB.14.2239[56] K. Hejazi, C. Liu, and L. Balents, “Landau levels in twisted bilayer graphene andsemiclassical orbits,” Physical Review B, vol. 100, p. 035115, Jul. 2019. [Online].Available: https://link.aps.org/doi/10.1103/PhysRevB.100.035115[57] C. R. Dean, L. Wang, P. Maher, C. Forsythe, F. Ghahari, Y. Gao, J. Katoch,M. Ishigami, P. Moon, M. Koshino, T. Taniguchi, K. Watanabe, K. L. Shepard,J. Hone, and P. Kim, “Hofstadter’s butterfly and the fractal quantum Hall effect inmoire´ superlattices,” Nature, vol. 497, no. 7451, pp. 598–602, May 2013. [Online].Available: https://doi.org/10.1038%2Fnature12186[58] A. Principi, M. Polini, G. Vignale, and M. I. Katsnelson, “Many-body orbitalparamagnetism in doped graphene sheets,” Physical Review Letters., vol. 104,p. 225503, Jun. 2010. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevLett.104.225503[59] H. Matsuura and M. Ogata, “Theory of orbital susceptibility in the tight-binding model: Corrections to the Peierls phase,” Journal of the PhysicalSociety of Japan, vol. 85, no. 7, p. 074709, Jul. 2016. [Online]. Available:https://doi.org/10.7566%2Fjpsj.85.074709[60] R. A. Z˙ak, D. L. Maslov, and D. Loss, “Spin susceptibility of interactingtwo-dimensional electrons in the presence of spin-orbit coupling,” Physical Review B,vol. 82, no. 11, Sep. 2010. [Online]. Available: https://doi.org/10.1103%2Fphysrevb.82.115415[61] M. Xie and A. MacDonald, “Nature of the correlated insulator states in twistedbilayer graphene,” Physical Review Letters, vol. 124, no. 9, Mar. 2020. [Online].Available: https://doi.org/10.1103%2Fphysrevlett.124.09760140[62] J. Diez-Merida, A. Diez-Carlon, S. Y. Yang, Y. M. Xie, X. J. Gao, K. Watanabe,T. Taniguchi, X. Lu, K. T. Law, and D. K. Efetov, “Magnetic Josephson junctionsand superconducting diodes in magic angle twisted bilayer graphene,” 2021. [Online].Available: https://arxiv.org/abs/2110.01067[63] J. V. Bustamante, N. J. Wu, C. Fermon, M. Pannetier-Lecoeur, T. Wakamura,K. Watanabe, T. Taniguchi, T. Pellegrin, A. Bernard, S. Daddinounou,V. Bouchiat, S. Gue´ron, M. Ferrier, G. Montambaux, and H. Bouchiat,“Detection of graphene’s divergent orbital diamagnetism at the Dirac point,”Science, vol. 374, no. 6573, pp. 1399–1402, 2021. [Online]. Available: https://www.science.org/doi/abs/10.1126/science.abf9396[64] K. D. Rasamani, F. Alimohammadi, and Y. Sun, “Interlayer-expanded MoS2,”Materials Today, vol. 20, no. 2, pp. 83–91, 2017. [Online]. Available:https://www.sciencedirect.com/science/article/pii/S1369702116302334[65] F. Zahid, L. Liu, Y. Zhu, J. Wang, and H. Guo, “A generic tight-binding modelfor monolayer, bilayer and bulk MoS2,” AIP Advances, vol. 3, 05 2013. [Online].Available: https://aip.scitation.org/doi/10.1063/1.480493641Appendix ACalculation of effective inter-layertunneling parameter w0In the main text, we mention that an effective inter-layer tunneling parameter w0 isextrapolated using a combined approach of empirical model and Slater-Koster method.Here, we present details on the calculation of w0.The inter-layer tunneling between two MoS2 monolayers with a small twist angle canbe conveniently modeled by the two-center approximation following the general recipeintroduced in previous studies on twisted bilayer graphene [18] and twisted homobilayerMoTe2 [23]. Particularly, electronic states near the conduction band minimum of mono-layer transition-metal dichalcogenides(TMDs) are predominantly from the 4dz2-orbitals ofthe transition metal atoms [27] and the inter-layer tunneling is expected to be mediatedby the σ-bonding between 4dz2-orbitals in each layer. To start with, we consider fermionicoperators that create Bloch states of the form:c†k,1,α =1√N∑Reik·Rc†1,α(R), (A.1)c†k′,2,β =1√N∑R′eik′·(R′+d0)c†2,β(R′ + d0).Here, k,R and k′,R′ denote the momentum and lattice vectors of layer 1 and 2, respec-tively. c†l,α(R) creates a localized Wannier orbital of dz2 character with spin α at site Rin layer l. In general, the inter-layer tunneling Hamiltonian is given by:HT =∑k,k′∑α,βc†k,1,αTαβ(k,k′)ck′,2,β + h.c. (A.2)42Appendix A. Calculation of effective inter-layer tunneling parameter w0where the tunneling matrix T reads:Tαβ(k,k′) = ⟨k, 1, α|HT |k′, 2, β⟩=1N∑R,R′e−ik·Reik′·(R′+d0) ⟨R, 1, α|HT |R′ + d0, 2, β⟩(two-center approximation) = − 1N∑R,R′e−ik·Reik′·(R′+d0)tαβ(R−R′ − d0)= − 1N∑R,R′,qe−ik·Reik′·(R′+d0) tαβ(q)Aeiq·(R−R′−d0)= − 1N∑qtαβ(q)Aei(k′−q)·d0(∑Rei(q−k)·R)(∑R′ei(k′−q)·R′)= −∑G,G′tαβ(k +G)Ωe−iG′·d0δk+G,k′+G′ (A.3)Here, tαβ(q) =∫d2rtαβ(r)e−iq·r is the Fourier transform of tαβ(r), A is the total areaof the system, Ω = A/N is the area of the primitive unit cell. In the last step of thederivation above, we made use of the identity∑R ei(q−k)·R = N∑G δq,k+G.To simplify the expression of Tαβ(k,k′) above, we note that for the relevant electronicstates near the K-points with k ∼ ξK, the dominant tαβ(q)-terms that enter the contin-uum model for valley ξ involve only 3 different reciprocal lattice vectors G = 0, ξG2, ξG3for the top layer (G′ = 0, ξG′2, ξG′3 accordingly for the bottom layer). Thus, the effectiveinter-layer tunneling Hamiltonian can be written as:HξT,eff = −1Ω∑k∑α,βc†k,1,αtαβ(ξK)ck,2,β + c†k,1,αtαβ(ξK + ξG2)e−iξG′2·d0ck+ξGM,2,2,β+ c†k,1,αtαβ(ξK + ξG3)e−iξG′3·d0ck+ξGM,3,2,β + h.c.= −∑k,αβc†k,1,αwαβck,2,β + c†k,1,αwαβe−iξG′2·d0ck+ξGM,2,2,β+ c†k,1,αwαβe−iξG′3·d0ck+ξGM,3,2,β + h.c., (A.4)where wαβ = tαβ(K)/Ω, GM,j = Gj −G′j . Explicitly, the tunneling matrices connectingstates at k on layer 1 to those at k,k + ξGM,2,k + ξGM,3 on layer 2 are listed below:Tαβ(k,k) = −wαβ,Tαβ(k,k + ξGM,2) = −e−iξG′2·d0wαβ,Tαβ(k,k + ξGM,3) = −e−iξG′3·d0wαβ. (A.5)For simplicity, we consider the case with d0 = 0 where metal sites from the two layers arealigned at zero twist angle, and assume wαβ = wδαβ given that the inter-layer tunneling isusually dominated by spin-preserving processes. This leads to the form of the inter-layertunneling Hamiltonian in Eq.3 of the main text.43Appendix A. Calculation of effective inter-layer tunneling parameter w0Figure A.1: Interlayer tunneling strength t(q) as a function of q calculated by numericalFourier transform. Dashed vertical line in black indicates the value w ≈ 10 meV extrapo-lated at q = K where Ka/π = 4/3.To evaluate the effective tunneling strength t(q = K) near the K-points, we use thefollowing simple model for t(r) which incorporates the Slater-Koster σ-bonding betweendz2-orbitals in the two different layers with an empirical exponential decay:t(r) = [n2(r)− (l2(r) +m2(r))/2]2Vddσ(r), (A.6)Vddσ(r) = t0e−λ√d2+|r|2 .Here, d = 0.7 nm is the inter-layer distance between two TMDmonolayers [64], n(r), l(r),m(r)are the directional cosines of the vector r connecting the two position vectors in the twolayers, withn2(r) =d2d2 + |r|2 , l2(r) +m2(r) =|r|2d2 + |r|2 . (A.7)The value of λ is estimated based on the empirical fact that t′ ≈ 0.1t, where t′ andt are the next-nearest-neighbor and neareast-neighbor hopping amplitudes. Thus, t/t′ =eλa ≈ 10, and λa = ln(10) ≈ 2.3. Note that the form of t(r) in Eq.A.6 is a radial functiont(r) ≡ t(r), and its Fourier transform has the general form of:t(q) =∫ ∞0drrt(r)∫ 2π0dθe−iqr cos θ (A.8)= 2π∫ ∞0drrJ0(qr)t(r),where J0(x) is the Bessel function of order zero. With the form of t(r) in Eq.A.6, t(q) can beobtained by numerical integration with the only fitting parameter t0. To fix t0, we note that44Appendix A. Calculation of effective inter-layer tunneling parameter w0states near the Γ-point of the valence band states in MoS2 also originate from 4dz2-orbitals[27], and the inter-layer coupling causes an energy splitting of ∆E ≃ 2t(q = 0) ≈ 800meV at the Γ-point in the valence bands in aligned bilayer MoS2 [65]. This impliest(q = 0) ≈ 400 meV and allows us to fix t0 and obtain the values of t(q) for all q as shownin Fig.A.1. The effective interlayer coupling strength at the K-point in our continuummodel is then extrapolated to be w0 ≈ 10 meV at qa = Ka = 4π/3 (dashed vertical linein black in Fig.A.1).45Appendix BSelf-consistent Hartree-FockcalculationsB.1 Hartree-Fock energy functional and LagrangeequationsWhen the Chern bands in a twisted bilayer MoS2 are half-filled, the low-energy effectiveHamiltonian is written as:Heff = H0,eff +HI,eff (B.1)where the non-interacting Hamiltonian is given byH0,eff =∑k∑ξ=±Eξ(k)c†ξ(k)cξ(k) (B.2)Here, ξ = ± is the valley index, Eξ(k) is the energy of moire´ band from valley ξ. Theeffective interacting Hamiltonian is given by:HI,eff =∑k,k′,qV (q)A∑ξξ′Λξ(k + q,k)Λξ′(k′ − q,k′) (B.3)× c†ξ(k + q)c†ξ′(k′ − q)cξ′(k′)cξ(k),where V (q) =∫d2rV (r)eiq·r is the Fourier transform of the two-body interaction V (r1−r2), A is the total area of the system, Λξ(k + q,k) ≡ ⟨uξ,k+q|uξ,k⟩ are the form factorsarising from projecting the general density-density interaction to the active bands, whichgive the overlap between the cell-periodic parts of Bloch states within the same band fromvalley ξ. Note that due to the Uv(1)-symmetry, wave functions from different valleys donot overlap: Λ+−(k′,k) ≡ ⟨u+,k′ |u−,k⟩ = 0.As mentioned in the main text, the Coulomb energy scale gC ∼ 10 meV in twisted bi-layer MoS2 overhelms the band widthW ∼ 1 meV in the Chern bands found in the contin-uum model. This motivates us to consider possible ground states at 1/2-filling as a generalHartree-Fock state of the form: |Φ⟩ = ∏k[w+(k)eiθ+(k)c†+(k) + w−(k)eiθ−(k)c†−(k)] |Φ0⟩,where wξ(k) and θξ(k) are real numbers and denote the amplitude and phase of the co-efficients of valley ξ at momentum k. |Φ0⟩ is the vaccuum state with all bands below theactive bands being fully filled. The variational parameters can be determined by minimiz-ing the ground state energy subject to the normalization condition w2+(k) + w2−(k) = 1.The energy of |Φ⟩ is given by:EΦ ≡ ⟨Φ|Heff |Φ⟩ = EΦ0 + EH + EΦF , (B.4)46B.1. Hartree-Fock energy functional and Lagrange equationswhere EΦ0 , EH and EΦF refer to the kinetic energy, Hartree energy, and the Fock exchangeenergy of |Φ0⟩, respectively. It is straightforward to show that the Hartree term essentiallyaccounts for the direct coupling: EH ∝ n2e, where ne is the electron density of the system,which is fixed at a given band filling. Thus, we neglect EH from now on and focus on EΦ0and EΦF :EΦ0 =∑kw2+(k)E+(k) + w2−(k)E−(k), (B.5)EΦF = −∑k,qV (q)A∑ξξ′Λξ(k + q,k)Λξ′(k,k + q)wξ(k + q)wξ′(k + q)wξ(k)wξ′(k)× e−iθξ(k+q)eiθξ′ (k+q)e−iθξ′ (k)eiθξ(k).Note that in order to minimize EΦF , the phase factors would adjust themselves such thate−iθξ(k+q)eiθξ(k)Λξ(k + q,k) = |Λξ(k + q,k)|. This reduces EΦF to the following form:EΦF = −∑k,qV (q)A∑ξξ′|Λξ(k + q,k)||Λξ′(k + q,k)|wξ(k + q)wξ′(k + q)wξ(k)wξ′(k) (B.6)where we have used the identity |Λξ′(k,k + q)| = |Λξ′(k + q,k)∗| = |Λξ′(k + q,k)|.To minimize EΦ0 + EΦF subject to the constraint w2+(k) + w2−(k) = 1, we introduce theLangrange multipliers λk such that∂wξ(k)[EΦ0 + EΦF −∑kλk(w2+(k) + w2−(k))]= 0. (B.7)This leads to the following two equations:λkw+(k) =[E+(k)−∑qg+(q,k)w2+(k + q)]w+(k)−[∑qg+−(q,k)w+(k + q))w−(k + q)]w−(k)λkw−(k) =[E−(k)−∑qg−(q,k)w2−(k + q)]w−(k)−[∑qg+−(q,k)w−(k + q)w−(k + q)]w+(k) (B.8)where we introduced the shorthand notation:{gξ(q,k) = V (q)|Λξ(k + q,k)|2/A, (ξ = ±)g+−(q,k) = V (q)|Λ+(k + q,k)||Λ−(k + q,k)|/A.(B.9)As we seek nontrivial solutions of Eq.B.8, it is convenient to introduce the followingtwo quantities:∆ξ(k) =∑qgξ(q,k)w2ξ (k + q), (B.10)∆+−(k) =∑qg+−(q,k)w+(k + q)w−(k + q),47B.2. Mean-field approximationwhich allows us to reformulate Eq.B.8 as an eigenvalue equation:(E+(k)−∆+(k) −∆+−(k)−∆+−(k) E−(k)−∆−(k))(w+(k)w−(k))= λk(w+(k)w−(k)). (B.11)Note that λk has the physical interpretation as the energy contribution from k to thetotal energy of the Hartree-Fock state. This motivates us to focus on the lower branch ofeigenvalues of the 2× 2 symmetric matrix:λ1(k) =12[E+(k) + E−(k)− (∆+(k) + ∆−(k))]−√∆2+−(k) +14[E+(k)− E−(k)− (∆+(k)−∆−(k))]2 (B.12)with the corresponding eigenvector [w+(k), w−(k)]T = [cos(ϕk/2), sin(ϕk/2)]T , wheresin(ϕk) =∆+−(k)√∆2+−(k) +14 [E+(k)− E−(k)− (∆+(k)−∆−(k))]2= 2 cos(ϕk/2) sin(ϕk/2),cos(ϕk) =12(∆+(k)−∆−(k))− (E+(k)− E−(k))√∆2+−(k) +14 [E+(k)− E−(k)− (∆+(k)−∆−(k))]2= cos2(ϕk/2)− sin2(ϕk/2). (B.13)The defining equations (Eq.B.10) require that ∆+(k),∆−(k),∆+−(k) satisfy the fol-lowing self-consistent equations:∆+(k) =1N∑k′J++(k′,k)1 + cos(ϕk′)2, (B.14)∆−(k) =1N∑k′J−−(k′,k)1− cos(ϕk′)2,∆+−(k) =1N∑k′J+−(k′,k)sin(ϕk′)2.Here, Jξξ′(k,k′) are the Fock exchange functions:Jξξ′(k′,k) = Ω−1∑GMV (k′ +GM − k)|Λξ(k′ +GM ,k)||Λξ′(k,k′ +GM )|, (B.15)where Ω is the area of the real-space moire´ unit cell, GM denotes the moire´ reciprocallattice vectors.B.2 Mean-field approximationTo simplify Eq.B.14, we employ the mean-field approximation which replaces the ex-change coupling strengths Jξξ′(k′,k) by their averages: intravalley exchange coupling J0 ≡48B.3. Hartree-Fock ground states⟨J++(k′,k)⟩ = ⟨J−−(k′,k)⟩, and the intervalley exchange coupling J1 ≡ ⟨J+−(k′,k)⟩.This reduces Eq.B.14 to a coupled set of self-consistent equations involving only constantorder parameters:∆+ =J0N∑k1 + cos(ϕk)2, (B.16)∆− =J0N∑k1− cos(ϕk)2,∆+− =J1N∑ksin(ϕk)2,wheresin(ϕk) =∆+−√∆2+− +14 [E+(k)− E−(k)− (∆+ −∆−)]2, (B.17)cos(ϕk) =12(∆+ −∆−)− (E+(k)− E−(k))√∆2+− +14 [E+(k)− E−(k)− (∆+ −∆−)]2.By adding and subtracting the first two equations in Eq.2.5, we end up with a compactform (Eq.5 of the main text):∆+ +∆− = J0, (B.18)∆+ −∆− = J02N∑k(∆+ −∆−)− δ(k)D(k),∆+− =J12N∑k∆+−D(k)Here, δ(k) = E+(k) − E−(k), D(k) =√∆2+− +14 [δ(k)− (∆+ −∆−)]2. Note that time-reversal enforces δ(−k) = E+(−k) − E−(−k) = E−(k) − E+(k) = −δ(k) to be an oddfunction of k.B.3 Hartree-Fock ground states under different regimes ofmean-field exchange coupling (J0, J1)It is worth noting that there exist two possible solutions for the last equation in Eq.B.18,which divides the whole set of solutions to Eq.B.18 into two different classes: (i) ∆+− = 0,and (ii) ∆+− ̸= 0. Particularly, in case (i) where ∆+− = 0, we have D(k) = 12 |δ(k) −(∆+ − ∆−)|, i.e., cos(ϕk) = ±1 for all k (or equivalently, either ∆+ = J0,∆− = 0 or∆+ = 0,∆− = J0). These two degenerate states are the valley-polarized (VP) solutionswhich spontaneously break time-reversal symmetry. The energy of the two degenerate VPstates is given byEV P =∑kE+(k)−NJ0. (B.19)49B.3. Hartree-Fock ground statesIn case (ii) where ∆+− ̸= 0, one can divide both sides of the 3rd equation in Eq.B.18by ∆+−, which simplifies the set of equations as below:∆+ +∆− = J0,∆+ −∆− = J0J1(∆+ −∆−)− J02N∑kδ(k)D(k),1 =J12N∑k1D(k), (B.20)where the three unknown mean-field order parameters ∆+,∆−,∆+− that characterize thisstate need to be solved self-consistently. In particular, we note that ∆+− ̸= 0 implies anontrivial inter-valley order associated with this Hatree-Fock state, in which the micro-scopic valley is no longer conserved (i.e., spontaneous breaking of valley-Uv(1) symmetry).On the other hand, any solution with nonzero ∆M ≡ ∆+ −∆− necessarily lifts the val-ley degeneracy and breaks time-reversal symmetry spontaneously. Thus, we call ∆M thevalley magnetic order. As it turns out, under time-reversal symmetry (e.g., in theabsence of any external magnetic fields), the only Uv(1)-breaking solution must satisfythe condition ∆+ = ∆−, i.e., there is no valley magnetic order ∆M = ∆+ −∆− involved.This solution is known as the (generalized) inter-valley coherence (IVC) state [30](note: w+(k) = w−(−k) ̸= 1/√2 in general except at time-reversal-invariant momenta).According to Eq.B.12, the energy of the IVC state is:EIV C =∑k12[E+(k) + E−(k)− (∆+ +∆−)]−D(k) (B.21)=∑k[E+(k)−D(k)]− N2J0,where we used E−(k) = E+(−k),∆++∆− = J0. The ground state can then be determinedby comparing the two energies EIV C and EV P :∆E = EIV C − EV P = N2J0 −∑kD(k), (B.22)where ∆+,∆−,∆+− involved in D(k) are given by the self-consistent solutions of Eq.B.20.The ground states obtained for each pair of (J0, J1) map out a J0 − J1 mean-field phasediagram.To develop some physical insights into the energetics of the system, we first considersome simple limits. First, we note that the flat band condition implies J0, J1 ≫ |δ(k)|.This combined with the self-consistency condition 1 = J12N∑k1D(k) suggests D(k) ∼ J1/2.Thus, when the intravalley exchange coupling dominates over its intervalley counterpart:J0 ≫ J1, we have ∆E ∼ N(J0 − J1)/2 ≫ 0. In this case, the VP state is favored.By a similar argument, when the intervalley exchange coupling dominates: J1 ≫ J0,∆E ∼ N(J0 − J1)/2 ≪ 0 and the IVC state always wins. In particular, in the entireJ1 ≥ J0 region, using the harmonic mean-arithmetic mean (HM-AM) inequality:N∑Ni=11xi≤∑Ni=1 xiN, (xi > 0,∀i = 1, ..., N) (B.23)50B.4. Physically accessible regimes of (J0, J1)one can show that ∑kD(k) ≥ N2∑k1D(k)=N2J1, (B.24)where we used the last equation in Eq.B.20. This shows that for J1 ≥ J0, ∆E ≤ N(J0 −J1)/2 ≤ 0 always holds, and the system favors the IVC state with nonzero ∆+− and∆M = ∆+ −∆− = 0.Next, we consider the J0 > J1 regime. As we discussed above, at J0 = J1, the systemalways favors the IVC state, while in the limit J0 ≫ J1 the system should favor the VPstate. This suggests that upon increasing the ratio J0/J1, a phase transition must happensomewhere in the J0 > J1 regime. By solving the self-consistent equations (Eq.B.18)numerically, we obtain the J0− J1 phase diagram in Fig.3 of the main text, which revealsthat the VP state is favored almost immediately when the J0 > J1 regime is accessed, anda first-order transition from IVC to VP phase occurs at J0 ≃ J1. This is due to the factthat the energetics of the system is essentially governed by (J0, J1) under the conditionJ0, J1 ≫ |δ(k)|. As a result, ∆E ≃ N(J0 − J1)/2 holds for almost all values of (J0, J1),such that the VP state is favored almost in the entire regime with J0 > J1.B.4 Physically accessible regimes of (J0, J1) underdensity-density repulsionsWe note that the discussions above are based on pure algebraic analysis — it is clearlydesirable to examine the accessible regimes of (J0, J1) on physical grounds. Indeed, givena generic density-density repulsive interaction which respects Uv(1) and T , the relationbetween J0 and J1 can be derived. In particular, the intra-valley exchange couplingstrength isJ0 ≡ ⟨J++(k′,k)⟩=1N2Ω∑k,k′∑GMV (k′ +GM − k)|Λ+(k′ +GM ,k)|2=1N2Ω∑k,k′∑GMV (−k′ −GM + k)|Λ−(−k,−k′ −GM )|2=1N2Ω∑k,k′∑GMV (k′ +GM − k)|Λ−(k′ +GM ,k)|2= ⟨J−−(k′,k)⟩ . (B.25)Note that in the equations above we used the fact that V (q) = V (−q) and the condi-tion Λ+(k′,k) = Λ−(−k,−k′) imposed by time-reversal symmetry. For the inter-valley51B.4. Physically accessible regimes of (J0, J1)exchange coupling strength,J1 ≡ ⟨J+−(k′,k)⟩=1N2Ω∑k,k′∑GMV (k′ +GM − k)|Λ+(k′ +GM ,k)||Λ−(k,k′ +GM )|≤ 1N2Ω∑k,k′∑GMV (k′ +GM − k) |Λ+(k′ +GM ,k)|2 + |Λ+(−k′ −GM ,−k)|22[V (q) = V (−q)] = 12N2Ω[(∑k,k′∑GMV (k′ +GM − k)|Λ+(k′ +GM ,k)|2)+(∑k,k′∑GMV (−k′ −GM + k)|Λ+(−k′ −GM ,−k)|2)]=122 ⟨J++(k′,k)⟩= J0. (B.26)Therefore, we identify the physically accessible points (J0, J1) are located in the J0 ≥J1 regime of the phase diagram. Notably, the relation J0 ≥ J1 derived above holdsregardless of details of the wave functions and the microscopic interaction V (q) as longas the interaction is repulsive in nature. Thus, as shown in Fig.3(a) of the main text, theVP state takes up a large portion of the physical regime with J0 ≥ J1 except a narrowstrip with J0 ≈ J1 where the IVC state is favored.52Appendix CGauge-invariant orbital magneticmoment and Berry curvatureEqs.(3.1) and (3.3) can be expressed in a gauge-invariant form by inserting the projectoronto the complete basis of eigenstates 1 =∑m |umk⟩ ⟨umk|, and using the general relation⟨umk|∂junk⟩ = ⟨umk| ∂jhk |unk⟩(εnk − εmk) (C.1)for m ̸= n. This identity moves the kj-derivatives from the eigenstates, which may haverandom phases when computed numerically, to the hermitian Hamiltonian. First, weexpress the out of plane orbital magnetic moment using the Levi-Civita tensor ϵijkmznk = −ie2ℏϵzij ⟨∂iunk| (hk − εnk) |∂junk⟩ , (C.2)and then insert projectors within the expectation value:ϵzij ⟨∂iunk| (hk − εnk) |∂junk⟩ = ϵzij∑m,m′⟨∂iunk|umk⟩ ⟨umk| (hk − εnk) |um′k⟩ ⟨um′k|∂junk⟩= ϵzij∑m ̸=n⟨unk| ∂ihk |umk⟩ (εmk − εnk) ⟨umk| ∂jhk |unk⟩(εnk − εmk)2 .(C.3)Thus applying this trick leads to the final gauge-invariant expressions:mznk =ie2ℏ∑m ̸=n⟨unk|hxk |umk⟩ ⟨umk|hyk |unk⟩ − ⟨unk|hyk |umk⟩ ⟨umk|hxk |unk⟩(εnk − εmk) , (C.4)andΩznk = i∑m̸=n⟨unk|hxk |umk⟩ ⟨umk|hyk |unk⟩ − ⟨unk|hyk |umk⟩ ⟨umk|hxk |unk⟩(εnk − εmk)2 , (C.5)where we have introduced a shorthand for the derivatives of the Hamiltonian hjk ≡ ∂jhk.53Appendix DDerivation of numerical formulafor orbital susceptibilityD.1 Tr′ Xˆ in the band basisAny single-particle GF has a simple form in the band basis and frequency space:g(k, ω) =∑m|umk⟩ ⟨umk|ω − εmk + iη (D.1)where |umk⟩ are eigenstates of the Hamiltonian hk, εmk the corresponding eigenvalues,and η > 0 a small parameter required to ensure convergence of the Fourier transform. Wewill usually omit η in the following, except when computing the energy integral. UsingEq.(3.26) and Tr′ Xˆ =∑n ⟨unk| Xˆ |unk⟩, we find the following expression:Tr′ Xˆ =∑n,m,l,p⟨n|hx |m⟩ ⟨m|hy |l⟩ ⟨l|hx |p⟩ ⟨p|hy |n⟩(ω − εn)(ω − εm)(ω − εl)(ω − εp) (D.2)To perform the integration over energy, we must split the quotient into different partsby partial fraction decomposition. Any of the εn being equal generally leads to differentdecompositions, therefore we split the trace into four separate sums:Tr{ghxghyghxghy} = S1 + S2 + S3 + S4 , (D.3)where the subscript indicates the highest exponent of the (ω − εn) factors in the denomi-nator.Setting n = m = p = l givesS4 =∑n(⟨n|hx |n⟩ ⟨n|hy |n⟩)2(ω − ϵn)4 =∑n(εxnεyn)2(ω − ϵn)4 . (D.4)Making just one index different (and including all permutations) we obtainS3 = 2∑n̸=m⟨n|hx |n⟩ ⟨n|hy |m⟩ ⟨m|hx |n⟩ ⟨n|hy |n⟩+ c.c.(ω − εn)3(ω − εm)= 2∑n̸=mεxnεyn⟨n|hx |m⟩ ⟨m|hy |n⟩+ c.c.(ω − εn)3(ω − εm) . (D.5)54D.1. Tr′ Xˆ in the band basisThe denominator of S3 can be decomposed as1(ω − εn)3(ω − εm) =1εn − εm1(ω − εn)3 −1(εn − εm)21(ω − εn)2+1(εn − εm)3(1ω − εn −1ω − εm),(D.6)and the ∼ 1ω terms can be combined into one by swapping n↔ m inside the sum.∑n̸=mεxnεyn⟨n|hx |m⟩ ⟨m|hy |n⟩+ c.c.(εn − εm)3(1ω − εn −1ω − εm)=∑n̸=mεxnεyn + εxmεym(εn − εm)3⟨n|hx |m⟩ ⟨m|hy |n⟩+ c.c.ω − εn .(D.7)We define the following expressions:C3,1 =2(εxnεyn + εxmεym)(εn − εm)3 (D.8)C3,2 =2εxnεyn(εn − εm)2 (D.9)C3,3 =2εxnεynεn − εm (D.10)Oxynm = ⟨n|hx |m⟩ ⟨m|hy |n⟩ (D.11)in order to write S3 in the more compact formS3 =∑n̸=m[C3,1ω − εn −C3,2(ω − εn)2 +C3,3(ω − εn)3]2ReOxynm . (D.12)S2 can be broken into two sub-cases, s2 and s′2. s2 is given bys2 =′∑n,m,l⟨n|hx |m⟩ ⟨m|hy |n⟩ ⟨n|hx |l⟩ ⟨l|hy |n⟩+ c.c.(ω − εn)2(ω − εm)(ω − εl)+ 2′∑n,m,lεxn ⟨n|hy |m⟩ ⟨m|hx |l⟩ ⟨l|hy |n⟩+ εyn ⟨n|hx |m⟩ ⟨m|hy |l⟩ ⟨l|hx |n⟩(ω − εn)2(ω − εm)(ω − εl) , (D.13)where we have introduced the notation∑′ indicating that none of the indices of the sumcan be equal to one another, and the denominator can be decomposed as1(ω − εn)2(ω − εm)(ω − εl) =12(εn − εm)(εn − εl)1(ω − εn)2 −1(εn − εm)(εn − εl)21ω − εn+1(εm − εl)(εm − εn)21ω − εm + (m↔ l) . (D.14)55D.2. Energy integralNote that swapping m ↔ l gives the complex conjugate of the expectation values in thenumerator. Now, s′2 is given bys′2 =∑n ̸=m(⟨n|hx |m⟩ ⟨m|hy |n⟩)2(ω − εm)2(ω − εn)2 +∑n̸=mεxnεxm ⟨n|hy |m⟩ ⟨m|hy |n⟩+ εynεym ⟨n|hx |m⟩ ⟨m|hx |n⟩(ω − εm)2(ω − εn)2(D.15)with the denominator decomposing as1(ω − εn)2(ω − εm)2 =1(εn − εm)21(ω − εn)2 −2(εn − εm)31ω − εn + (n↔ m) . (D.16)Likewise, swapping n ↔ m inside the sum in s2 adds a complex conjugate. Now, we seethat setting m = l in the first two terms of Eqs.(D.13-D.14) yields Eqs.(D.15-D.16), thuswe can combine them into a single sum without the m ̸= l restriction.S2 =∑n(∑m,l{OxynmOxynl + εxnOyxynml + εynOxyxnml}××{1(εn − εm)(εn − εl)1(ω − εn)2 −2(εn − εm)(εn − εl)21ω − εn}+∑m ̸=l{OxynmOxynl + εxnOyxynml + εynOxyxnml}{ 2(εm − εl)(εm − εn)21ω − εm}+ c.c.)(D.17)Where Oyxynml = ⟨n|hy |m⟩ ⟨m|hx |l⟩ ⟨l|hy |n⟩ and Oxyxnml = ⟨n|hx |m⟩ ⟨m|hy |l⟩ ⟨l|hx |n⟩. Thethird term of Eq.(D.14) must be treated separately since m = l results in division by zero.Finally, the summand of S1 is identical to that of Eq.(D.2) with all indices different,and the partial fraction decomposition is1(ω − εn)(ω − εm)(ω − εl)(ω − εp) =1(εn − εm)(εn − εl)(εn − εp)1ω − εn+permutations .(D.18)Swapping indices so that all terms are ∼ 1ω−εn we obtainS1 = 2′∑n,m,l,p⟨n|hx |m⟩ ⟨m|hy |l⟩ ⟨l|hx |p⟩ ⟨p|hy |n⟩+ c.c.(εn − εm)(εn − εl)(εn − εp)1ω − εn (D.19)D.2 Energy integralNow, to evaluate the integral we invoke the following identity:Im∫ ∞−∞dωf(ω)ω − εn + iη =∫ ∞−∞dωf(ω) (−πδ(ω − εn)) = −πf(εn) , (D.20)where f(ω) is an arbitrary real function. From Eq.(D.20), we can also show via integrationby parts thatIm∫ ∞−∞dωf(ω)(ω − εn + iη)2 = −πf′(εn) , (D.21)56D.2. Energy integralwhere f ′(ω) is the derivative with respect to ω of f , and f ′(εn) is its value at εn. Eqs.(D.20)and (D.21) can be generalized to any power of the denominator by performing repeatedintegration by parts, yielding the general identity:Im∫ ∞−∞dωf(ω)(ω − εn + iη)p = −π(p− 1)!f(p−1)(εn) (D.22)We now compute the energy integral of each sum comprising Tr Xˆ.S4 and S1 contain just one term and the integrals are straightforward to compute.Im∫ ∞−∞dω nF (ω)S4 = Im∫ ∞−∞dω nF (ω)∑n(εxnεyn)2(ω − ϵn)4= − π3!∑nn′′′F (εn)(εxnεyn)2 (D.23)The contribution of S4 to the overall susceptibility is:χ4 =µ0e212ℏ2∑n∫BZdk(2π)dn′′′F (εn)(εxnεyn)2 (D.24)For S1 the energy integral givesIm∫ ∞−∞dω nF (ω)S1 = −2π′∑n,m,p,lnF (εn)⟨n|hx |m⟩ ⟨m|hy |l⟩ ⟨l|hx |p⟩ ⟨p|hy |n⟩+ c.c.(εn − εm)(εn − εl)(εn − εp)(D.25)and susceptibility contributionχ1 =µ0e2ℏ2∑n∫BZdk(2π)dnF (εn)′∑m,l,p⟨n|hx |m⟩ ⟨m|hy |l⟩ ⟨l|hx |p⟩ ⟨p|hy |n⟩+ c.c.(εn − εm)(εn − εl)(εn − εp) .(D.26)The energy integral of S3 contains three different terms:Im∫ ∞−∞dω nF (ω)S3= −4π∑nεxnεyn[∑mnF (εn)− nF (εm)(εn − εm)3 Re ⟨n|hx |m⟩ ⟨m|hy |n⟩− n′F (εn)gnxyk + n′′F (εn)ℏ2eReMnk],(D.27)so that its contribution to the susceptibility isχ3 =2µ0e2ℏ2∑n∫BZdk(2π)d[nF (εn)∑mεxnεyn + εxmεym(εn − εm)3 Re ⟨n|hx |m⟩ ⟨m|hy |n⟩− n′F (εn)gnxyk + n′′F (εn)ℏ2eReMnk](D.28)57D.2. Energy integralThe energy integral of S2 is given byIm∫ ∞−∞dωnF (ω)S2 =− 2π∑n(∑m,l{OxynmOxynl + εxnOyxynml + εynOxyxnml}××{n′F (εn)2(εn − εm)(εn − εl) −nF (εn)(εn − εm)(εn − εl)2}+∑m ̸=l{OxynmOxynl + εxnOyxynml + εynOxyxnml} nF (εm)(εm − εl)(εm − εn)2 + c.c.)(D.29)Giving the susceptibility contributionχ2 =µ0e2ℏ2∑n∫BZdk(2π)d(∑m,l{OxynmOxynl + εxnOyxynml + εynOxyxnml}××{n′F (εn)2(εn − εm)(εn − εl) −nF (εn)(εn − εm)(εn − εl)2}+∑m̸=l{OxynmOxynl + εxnOyxynml + εynOxyxnml} nF (εm)(εm − εl)(εm − εn)2 + c.c.)(D.30)We now isolate the separable terms from χorb,2, consisting of the OxynmOxynl terms in thefirst line of equation D.30χseporb,2 =µ0e2ℏ2∑n∫BZdk(2π)d[n′F (εn)(12∑m⟨n|hx |m⟩ ⟨m|hy |n⟩εn − εm∑l⟨n|hx |l⟩ ⟨l|hy |n⟩εn − εl + c.c.)︸ ︷︷ ︸≡A1− nF (εn)(∑m⟨n|hx |m⟩ ⟨m|hy |n⟩εn − εm∑l⟨n|hx |l⟩ ⟨l|hy |n⟩(εn − εl)2 + c.c.)︸ ︷︷ ︸≡A2](D.31)A1 = Re∑m̸=n⟨n|hx |m⟩ ⟨m|hy |n⟩εn − εm2(D.32)A2 = 2Re∑m ̸=n⟨n|hx |m⟩ ⟨m|hy |n⟩εn − εm∑l ̸=n⟨n|hx |l⟩ ⟨l|hy |n⟩(εn − εl)2 (D.33)To help us simplify the above, we recall the definition of the quantum geometric tensorintroduced by BerryTnij(k) =∑m ̸=n⟨n|hi |m⟩ ⟨m|hj |n⟩(εn − εm)2 , (D.34)58D.2. Energy integralwhose real part is the quantum metric tensor gnij(k) = ReTnij and whose imaginary partis related to the Berry curvature as:Ωnk = i∑m ̸=n⟨n|hx |m⟩ ⟨m|hy |n⟩ − ⟨n|hy |m⟩ ⟨m|hx |n⟩(εn − εm)2= −2Im∑m ̸=n⟨n|hx |m⟩ ⟨m|hy |n⟩(εn − εm)2= −2ImTnij . (D.35)In analogy with Tnij , we also define the following complex tensor:Mnk =eℏ∑m ̸=n⟨n|hx |m⟩ ⟨m|hy |n⟩εn − εm (D.36)which is related to the orbital magnetic moment asmnk =ie2ℏ∑m ̸=n⟨n|hx |m⟩ ⟨m|hy |n⟩ − ⟨n|hy |m⟩ ⟨m|hx |n⟩εn − εm= − eℏIm∑m ̸=n⟨n|hx |m⟩ ⟨m|hy |n⟩εn − εm= −ImMnk . (D.37)Using these definitions, and the identity Re(XY ) = ReXReY − ImXImY we find that:A1 =ℏ2e2((ReMnk)2 − (ImMnk)2)=ℏ2e2((ReMnk)2 −m2nk), (D.38)A2 =ℏe(2gnxykReMnk −mnkΩnk) , (D.39)which allows us to write the separable part in compact formχsep2 = µ0∑n∫BZdk(2π)d{n′F (εn)[(ReMnk)2 −m2nk]+ nF (εn)[eℏmnkΩnk − 2eℏ gnxyReMnk]}.(D.40)The inseparable part of χ2 isχ2 − χsep2 =µ0e2ℏ2∑n∫BZdk(2π)d{∑m,l(εxnOyxynml + εynOxyxnml)××(n′F (εn)2(εn − εm)(εn − εl) −nF (εn)(εn − εm)(εn − εl)2)+∑m ̸=l(OxynmOxynl + εxnOyxynml + εynOxyxnml) nF (εm)(εm − εl)(εm − εn)2 + c.c.}.(D.41)59D.3. Numerical evaluation of χorbIf we swap indices n ↔ m in the second line and compare to the terms of χ1, we seethat these are a special case of the inseparable term χ1 with different permutations of them, l, p indices equal to one another. Furthermore, if all three are equal we simply get the∼ εxmεym term in χ3. This reduces the total susceptibility toχorb = χsep2 +µ0e2ℏ2∑n∫BZdk(2π)d{n′′′F (εn)(εxnεyn)212− n′F (εn)2εxnεyngnxyk + n′′F (εn)ℏeεxnεynReMnk+ nF (εn)[ ∑m,l,p⟨n|hx |m⟩ ⟨m|hy |l⟩ ⟨l|hx |p⟩ ⟨p|hy |n⟩+ c.c.(εn − εm)(εn − εl)(εn − εp)+∑mεxnεyn(εn − εm)3 (⟨n|hx |m⟩ ⟨m|hy |n⟩+ c.c.)]+∑m,l(εxnOyxynml + εynOxyxnml + c.c.) [ n′F (εn)2(εn − εm)(εn − εl) −nF (εn)(εn − εm)(εn − εl)2]}.(D.42)Finally, rearranging to group terms according to derivatives of the Fermi distributionwe obtain our final expression for the orbital susceptibility of a separable tight bindingHamiltonian:χorb =µ0e2ℏ2∑n∫BZdk(2π)d{n′′′F (εn)(εxnεyn)212+ n′′F (εn)ℏeεxnεynReMnk− n′F (εn)[ℏ2e2(m2nk − (ReMnk)2)+ 2εxnεyngnxyk −∑m,lεxn ⟨n|hy |m⟩ ⟨m|hx |l⟩ ⟨l|hy |n⟩+ (x↔ y)(εn − εm)(εn − εl)]+ nF (εn)[ℏe(mnkΩnk − 2gnxykReMnk) + εxnεyn∑m⟨n|hx |m⟩ ⟨m|hy |n⟩+ c.c.(εn − εm)3−∑m,lεxn ⟨n|hy |m⟩ ⟨m|hx |l⟩ ⟨l|hy |n⟩+ (x↔ y) + c.c.(εn − εm)2(εn − εl)+∑m,l,p⟨n|hx |m⟩ ⟨m|hy |l⟩ ⟨l|hx |p⟩ ⟨p|hy |n⟩+ c.c.(εn − εm)(εn − εl)(εn − εp)]},(D.43)with no restriction on the sums, except that no index is equal to n.For easier comparison to other works which derive general formulas for orbital suscep-tibility in terms of geometric quantities, such as Refs. [42, 43], the terms in the first linecan be integrated by parts to reduce the degree of the nF (εn) derivatives.D.3 Numerical evaluation of χorbThe inseparable χinter contribution presents significant challenges for numerical compu-tation due to the sums over multiple band indices. Calculating the contribution on the60D.3. Numerical evaluation of χorbsecond line of Eq.(3.45) for a model with ∼ 100 bands would require a sum over one millionterms. Instead, we show that all of these terms can be computed efficiently using matrixoperations.We first define the interband velocity matrix Vi, whose elements are expectation valuesof hi:Vinm = ⟨n|hi |m⟩ . (D.44)To illustrate its structure, we write it out explicitly as follows:Vi =εi0 ⟨0|hi |1⟩ . . . ⟨0|hi |N⟩⟨1|hi |0⟩ εi1....... . . ⟨N − 1|hi |N⟩⟨N |hi |0⟩ . . . ⟨N |hi |N − 1⟩ εN . (D.45)We also define the skew-symmetric matrix D, given by:Dnm ={0 if n = m1εn−εm otherwise. (D.46)Using the notation ⊙ for element-wise multiplication (also known as Hadamard product),and juxtaposition for the conventional matrix product, the above objects can be used toexpress, for example, the quantum geometric tensor((Vx ⊙D)(Vy ⊙D))nn = −Tnxy(k) . (D.47)Similarly, we can obtain the interband term involving a sum over three band indices fromthe following matrix product:∑m,l,p ̸=n⟨n|hx |m⟩ ⟨m|hy |l⟩ ⟨l|hx |p⟩ ⟨p|hy |n⟩+ c.c.(εn − εm)(εn − εl)(εn − εp) = 2Re[(((((Vx⊙D)Vy)⊙D)Vx)⊙D)Vy]nn(D.48)We can use a combination of matrix/Hadamard products of Vi and D to compute allcontributions to the susceptibility. This method greatly reduces the time required tocompute the sums over bands, making the susceptibility formula in Eq.(3.37) numericallytractable.61
Thesis/Dissertation
eng
Physics
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Franz, Marcel
Orbital susceptibility and magnetic switching in twisted moiré heterostructures
Text