Quantum Magnetism in Wannier-Obstructed Mott Insulators

We develop a strong coupling approach towards quantum magnetism in Mott insulators for Wannier obstructed bands. Despite the lack of Wannier orbitals, electrons can still singly occupy a set of exponentially-localized but nonorthogonal orbitals to minimize the repulsive interaction energy. We develop a systematic method to establish an effective spin model from the electron Hamiltonian using a diagrammatic approach. The nonorthogonality of the Mott basis gives rise to multiple new channels of spin-exchange (or permutation) interactions beyond Hartree-Fock and superexchange terms. We apply this approach to a Kagome lattice model of interacting electrons in Wannier obstructed bands (including both Chern bands and fragile topological bands). Due to the orbital nonorthogonality, as parameterized by the nearest neighbor orbital overlap $g$, this model exhibits stable ferromagnetism up to a finite bandwidth $W\sim U g$, where $U$ is the interaction strength. This provides an explanation for the experimentally observed robust ferromagnetism in Wannier obstructed bands. The effective spin model constructed through our approach also opens up the possibility for frustrated quantum magnetism around the ferromagnet-antiferromagnet crossover in Wannier obstructed bands.


I. INTRODUCTION
Mott insulators are correlated insulators where electrons singly occupy localized orbitals to avoid the repulsive interaction.In many cases, Mott insulators further develop antiferromagnetic ordering below the charge gap due to the superexchange interaction among lowenergy spin degrees of freedoms.However, in recent twisted bilayer graphene (tBLG) experiments 1,2 , the observation of ferromagnetic hysteresis 3 suggests that the three-quarter-filling Mott insulating state exhibits ferromagnetism.In another experiment on twisted double bilayer graphene (tDBLG) [4][5][6] , an increasing gap under inplane magnetic field also suggests a ferromagnetic phase.Such ferromagnetic Mott state is understood in the flat band limit [7][8][9][10] , where the spin exchange term in Coulomb interaction reduces energy of the spin-polarized state.0][21][22][23][24][25][26] While most analyses consistently point to a robust ferromagnetism near the flat-band limit, it remains challenging to analyze the instability of such ferromagnetic state in competition with the antiferromagnetic superexchange away from the flat-band limit.
One major obstacle to model the competition between exchange and superexchange effects systematically in tBLG has to do with the Wannier obstruction, 27,28 i.e. if Wannier orbitals are constructed with only the relevant bands, they cannot respect all the symmetries. 29,30To construct the Wannier orbitals for tBLG with all symmetries taken into account, the minimal model needs to contain at least ten bands as pointed out by H. C. Po et al. 19,31,32 If we only focus on a few bands within the energy scale of interaction, the Wannier obstruction would prevent us from constructing Wannier orbitals.In lack of Wannier orbitals, it becomes unclear how the electrons should localize to form Mott insulators, which further obscure the derivation of low-energy effective spin (and/or valley) models.The issue of Wannier obstruction has long been identified and studied in quantum Hall systems and other Chern insulators.In the phases with nonzero Chern number, exponentially-localized Wannier orbitals cannot be constructed regardless any symmetry considerations. 33Now with more fragile topological systems being discovered, [34][35][36] the Wannier obstruction becomes a more general issue in the study of strongly correlated electronic systems.
We begin our general discussion by considering an extended Hubbard model on an arbitrary lattice, where c i = (c i,↑ , c i,↓ ) is the electron operator containing spin degrees of freedom and n i = c † i c i is the total electron number on site i.The model may be generalized to include orbital or valley 31 degrees of freedom, but in this work we will only focus on spins for illustration purpose.We assume certain degree of locality in t ij and U ij .Suppose that H t produces several bands, described by the dispersion relations n (k), In many cases, we are only interested in a subset of these bands near the Fermi energy with total bandwidth W , and well separated from other bands by energy gap ∆.Although the full band structure has a tight-binding model description, a subset of these bands may not.H. C. Po et al. 27,28 discussed several scenarios for a subset of bands, which can be trivial, obstructed trivial, fragile topological, or stably topological.As long as these bands are not trivial, they admit the Wannier obstruction, which is an obstruction towards constructing symmetric, exponentially-localized and orthogonal Wannier orbitals by linearly combining Bloch states within the subset of bands.The well known ones are with Chern number where no additional symmetry is required. 33,37,385][36] Finally, for the obstructed trivial bands, the obstruction is only due to the fact that the Wannier center does not reside on a lattice site, which can be resolved by adding empty sites.
The absence of such Wannier orbitals prevents us from writing down an effective tight-binding model targeting only those subset of bands of interest.This is precisely the case for tBLG, where the relevant conduction and valance bands are fragile topological and hence Wannier obstructed by the C 2 T symmetry. 19,31,32eak-coupling approaches have been developed to treat the interaction perturbatively (as U W ), 19,20,23 such that one only need to work with a momentum space description of the effective band structure near the Fermi surface, hence circumventing the Wannier obstruction.However, it remains challenging to understand the strong-coupling physics in Wannier obstructed bands, when the energy scales are arranged in the following hierarchy The band gap ∆, as the leading energy scale, protects a set of Wannier obstructed low-energy bands from mixing with high-energy bands away from the Fermi surface.
To further respect the interaction energy U , electrons should repel each other into real-space localized orbitals by combining states in the low-energy bands.In the standard notion of Mott insulators, electrons are localized in Wannier orbitals with charge fluctuation gapped and spin fluctuation remained active at low energy.Now in the absence of Wannier orbitals, does the many-body Hamiltonian in Eq. ( 1) still admits Mott-like ground states?If so, can we write down the trial wave function to describe the low-lying states?Can we derive an effective model to describe the spin dynamics at low energy?Motivated by these questions, we take a closer look at the requirements of Wannier orbitals, namely symmetry, locality and orthogonality.If we sacrifice one or more of them, the Wannier obstruction can be lifted and it would be possible to construct orbitals to host electrons in a Mott-like state.If we sacrifice the symmetry requirement, we will have to fine tune hopping parameters to fit the band structure.If we sacrifice the locality requirement, we will end up with a non-local hopping model which is hard to deal with.So we decided to explore the possibility of sacrificing orthogonality and working with a set of nonorthogonal Wannier basis.This approach is in analogous to the Maki-Zotos wavefunction 39 and the von Neumann lattice formulation [40][41][42] in quantum Hall systems with nonzero but negligible orbital overlap.
We present the general theory based on nonorthogonal Wannier basis with finite orbital overlap in Sec.II.The nonorthogonality of the orbitals leads to new spin exchange channels and chiral spin exchange channels.The competition among these new channels could lead to a rich magnetic phase diagram.We discuss leading order contributions to the effective spin Hamiltonian in Sec.III and discuss several new channels emerging from the theory of nonorthogonal basis.We also propose an energetic objective function in Sec.IV to construct these nonorthogonal orbitals numerically to facilitate the study of any concrete model.Within this framework we study a toy model proposed in Ref. 28 in Sec.V to demonstrate our framework in Wannier-obstructed bands.In this model, we show that the flat-band ferromagnetism remains stable up to finite band width and a variety of magnetic phases appear around the the ferromagnetantiferromagnet crossover.

II. THEORY OF NONORTHOGONAL BASIS
In this section, we discuss how to project the Hamiltonian in Eq. (1) to a nonorthogonal spin basis of lowenergy states and how to treat perturbative corrections.Let us assume a set of localized and normalized but nonorthogonal orbitals φ I (i) in the real space labeled by the orbital index I, which jointly labels the unit cell and the orbital within the unit cell.We will leave the energetic criterion to optimize these orbitals and their completeness as a set of basis for later discussions in Sec.IV.For now, we assume that electrons will self-organize under repulsive interaction to develop such nonorthogonal localized orbitals.The nonorthogonality implies a nontrivial metric g IJ among these orbitals, Nevertheless, we assume the normalization condition g II = 1 for all I.Given these orbitals, we can define a set of fermion operators a † Iσ that create electrons residing on these orbitals, such that they satisfy the following anticommutation relations {a Iσ , a Jτ } = {a † Iσ , a † Jτ } = 0 and {a Iσ , a † Jτ } = g IJ δ στ .These operators a † Iσ allow us to construct a set of many-body trial states from the vacuum state |vac , In these many-body states, every orbital is singly occupied and the spin configuration is labelled by σ = {σ I }.The nonorthogonality of single particle orbitals also implies the nonorthogonality of these many-body states Ψ σ |Ψ τ = δ στ .In the Mott limit U W , the trial states |Ψ σ span the low energy manifold H.All the charge degrees of freedom are frozen, and the spin degrees of freedom are still allowed to fluctuate, resembling the Mott states.We can then derive the effective theory for these spin degrees of freedom and investigate the resulting phases.

A. Exchange Interactions and Beyond
To formulate an effective spin model, we first introduce the spin Hilbert space H spanned by a fictitious set of orthogonal Ising basis |σ , which allows us to define the spin operator S I = (S x I , S y I , S z I ) in the conventional way , where σ a τ σ denotes the Pauli matrix element.We would like to comment that the Ising states |σ do not directly correspond to physical electronic states, but merely play a bookkeeping role to provide a convenient basis for the purpose of representing spin operators.The relation between different Hilbert spaces is illustrated in Fig. 1.

Full many-body space
(orthogonal) Low-energy space ℋ = span {Ψ σ 〉} (nonorthogonal) 1.The lattice model in Eq. ( 1) is defined in the full many-body Hilbert space spanned by the orthogonal Fock states of electrons, which includes a low-energy subspace H spanned by the nonorthogonal trial states |Ψσ in Eq. (6).A spin Hilbert space H spanned by the orthogonal Ising basis |σ is introduced to represent the effective spin model.The linear map A (and A † ) connects H and H.
We want to project the many-body Hamiltonian Eq. (1) to the low energy subspace H of electrons and then translate it to the spin space H.The complication arises from the nonorthogonality of the many-body trial states |Ψ σ that span the low energy subspace H.Here we present a systematic approach to deal with the nonorthogonality.First we introduce a non-unitary linear transformation A : H → H to map the Ising basis |σ to the many-body trial states |Ψ σ with correspond-ing spin configuration, Next we define the adjoint operator where {| Ψσ } is the dual basis to {|Ψ σ } such that Ψσ |Ψ τ = δ στ .Now we can use the operators A and A † to project the identity operator 1 and Hamiltonian H to the Ising basis, Under the projection, the original eigen problem H|Ψ = E|Ψ becomes a generalized eigen problem by inserting the identity operator AA −1 and multipling A † from the left, where |Φ ≡ A −1 |Ψ is the representation of the manybody eigenstate |Ψ in the Ising basis.Now the generalized eigen problem in Eq. ( 9) is formulated in the spin Hilbert space H with a nice orthogonal basis.
After the projection, we can expand the many-body Gram matrix G and the projected Hamiltonian H as linear combinations of permutation operators χ P in the spin space, where S N denotes the permutation group over all orbitals I, (−) P is the sign of the permutation P, and χ P is the permutation operator that permutes the spins among different orbitals χ P |{σ I } = |{σ P −1 (I) } .For two-spin and three-spin permutations, χ P can be expressed with the familiar Heisenberg term and chiral spin term, To specify the coefficients G P and H P , we introduce the hopping tensor t IJ and the interaction tensor U IJKL , where t ij and U ij are the bare hopping and interaction coefficients in H, as introduced in Eq. (1).Given the tensors g IJ in Eq. ( 4) and t IJ , U IJKL in Eq. ( 12), the coefficients G P , H P ≡ T P + U P in Eq. ( 10) are given by where T P and U P denotes the contribution from the hopping and the interaction terms respectively.
To simplify the notation, we introduce the diagrammatic representation: a closed loop of arrows represents a permutation cycle among the orbitals.If an arrow from I to J is labeled by t or g, it contributes a factor of t IJ or g IJ ; if two arrows, e.g. one from I to J and the other from K to L, are connected and labeled by U , it contributes a factor of U IJKL .For example, Using the diagrammatic representations, one can expand G and H in terms of spin permutations χ P order by order.To the order of two-spin exchange, we have where the permutation operator χ (IJ) can eventually be written in terms of spin operators as in Eq. ( 11).Higher order permutations can be included systematically based on Eq. (10).For all summations appeared in Eq. ( 15), it is assumed that the orbitals labeled by different indices do not coincide.In the orthogonal limit g IJ = δ IJ , all diagrams in Eq. ( 15) that contains g-labeled (red) arrows will vanish, such that G reduces back to the identity operator 1 and H reduces to the following three terms The first diagram is the band energy.The second and third diagrams are respectively the Hartree and the Fock energies between electrons from orbitals I and J (where the Fock interaction is accompanied with the spin exchange χ (IJ) ).For nonorthogonal orbitals, the metric g IJ becomes non-trivial, then non-vanishing terms (containing g-labeled arrows) in Eq. ( 15) suggest contributions beyond Hartree-Fock approximation, which lead to new channels of spin exchange interactions.Furthermore, there exist three-and even more spin interactions that are conventionally only present through the interactionsuppressed super-exchange effects.The various spin permutation interactions competing with each other could result in a frustrated quantum magnet with a rather rich phase diagram.To determine the ground state of the lowenergy spin degrees of freedom in Wannier-obstructed Mott insulators, one will need to solve the generalized eigen problem in Eq. ( 9).

B. Superexchange Interactions from Perturbation
In the above discussion, we project the many-body Hamiltonian H to the low-energy subspace H spanned by the Mott states |Ψ σ to derive the effective spin model in the flat-band limit W U .Away from the flat-band limit, the electrons can virtually hop to the neighboring orbitals and back, which gives rise to the superexchange interactions among the spins.To capture such effect, we should go beyond the low-energy subspace H and consider the perturbation effects in orders of (t/U ).
In the following, we analyze the perturbative corrections to the effective spin model within our framework.In Appendix A, we review the generalized perturbation theory for nonorthogonal basis.We apply the perturbation theory to the spin dynamics in the low energy manifold H by treating the hopping Hamiltonian H t in Eq. ( 2) as a small perturbation. 43Consider the high energy subspace spanned by states |n, α with double occupancy but still the same number of electrons, where n labels the number of doubly occupied orbitals and α labels the configuration.Assuming that the zeroth order energy is completely determined by n and states with different n's are orthogonal, i.e. n, α|m, β = δ mn G nαβ , we get the following energy correction from the second-order perturbation theory Since the dominant contribution of Ψ σ |H t |n, α comes from states with only one doubly occupied site, we can approximate Eq. ( 17) by where U is the energy required to create one double occupancy, and 1 Ψ = σ |Ψ σ Ψσ | is the projection operator to the low energy subspace.Now we can use the projec-tion introduced previously to write down the Hamiltonian H(2) ≡ A † H (2) A in the Ising basis, where H 2 t ≡ A † H 2 t A follows our convention, and in the second term we use where • denotes the composition of permutations (such that χ P•Q = χ P χ Q ), and the coefficients T P and T 2 P are given by and To the order of twospin exchange, we have For all summations appeared in Eq. ( 22), it is assumed that orbitals labeled by different indices do not coincide.
In the orthogonal limit g IJ = δ IJ , all diagrams containing g-labeled arrows will vanish, such that the secondorder perturbation contains only the usual t 2 /U antiferromagnetic superexchange interaction.When we allow non-trivial g IJ , both T P and T 2 P can give rise to new channels in spin interactions.Mediated by g IJ , three-or higher order spin interactions can also arise even at the level of secondorder perturbation in (t/U ).Similar treatment can be generalized to higher order perturbation theory.

III. EFFECTIVE SPIN HAMILTONIAN
In this section, we discuss how to solve the effective spin model constructed in Sec.II.There are two major challenges.First, the model is presented as a generalized eigenvalue problem in Eq. ( 9), which requires to diagonalize a complicated operator G −1 H that is not guaranteed to be short-ranged on the lattice.Second, the summation over all permutation in Eq. ( 10) is hard to track even numerically.Both challenges can be resolved by separating connected permutations and disconnected permutations.A connected permutation is formally defined as a cyclic permutation (or cycle) in group theory, while those that are not cycles are called disconnected in the following text.As demonstrated in Eq. ( 15), diagrams in H can be organized by the connected diagram (containing t or U ), each followed by a series of disconnected diagrams (containing g only).The series of disconnected diagrams is similar to G, which motivates us to factor G out of H.However, residue terms are generated due to over counting diagrams with colliding indices, illustrated as follows where P 0 sums over single-cycle (i.e.connected) permutations and P 1,2,••• sum over those permutations that have non-vanishing index overlap with every other permutation (including P 0 ) in the diagram.The explicit expression and a detailed convergence analysis of the entire series can be found in Appendix B. The rough idea is that the expansion is controlled by the small parameter g ≡ |g IJ | 1 for well-localized orbitals φ I .Since n-spin interaction can only be generated by χ P with length(P) ≥ n, they are suppressed by at least g n−2 .Thus the spin dynamics is still dominated by few-spin interactions as expected.Furthermore, for these few-spin interactions of small n, the contribution from sub-leading terms in the series Eq. ( 24) is further suppressed by ng 2 compared to the leading g n−2 term.Thus, for those dominating few-spin interactions, we only need to consider the leading order connected diagrams: where S * N is the set of single-cycle permutations in S N .With this approximation, the general eigenvalue problem in Eq. ( 9) reduces to the ordinary eigenvalue problem where H c collects the connected pieces in H.The problem reduces to solving the effective spin Hamiltonian H c on a set of orthogonal basis |σ , for which many welldeveloped analytical and numerical tools in quantum magnetism can be applied.Similar treatment applies to the perturbation theory described in Sec.II B, where the G −1 in Eq. ( 19) has cancelled the disconnected pieces in one of the Ht 's.Then the effective spin Hamiltonian H c in Eq. ( 26) gets corrected by H c .In conclusion, despite of the nonorthogonality of the Mott basis and the complication of the generalized eigen problem, we can still work with an effective spin Hamiltonian H c in a ordinary eigen problem to describe the low-energy spin degrees of freedoms approximately.
The full spin-rotation symmetry dictates the spin Hamiltonian H c to take the general form of up to three-spin interactions.Here we keep only the near neighbor interactions and analyze the coupling strengths of the Heisenberg interaction S I • S J and the chiral spin interaction S I • (S J × S K ).In Eq. ( 15), among the terms attached with χ IJ , to the leading order of g, those contribute to the Heisenberg interaction are The second term is the familiar Fock exchange term, which is always positive and thus provides inter-site Hunds coupling.We remark that the Fock term is non-vanishing even in the orthogonal limit due to the density-density overlap between orbitals. 13,17The first term is a new channel arising from nonorthogonality, which can be either ferromagnetic or antiferromagnetic depending on its sign.This channel could potentially provide a stronger antiferromagnetism in the strong coupling (large U ) limit than the usual t 2 /U superexchange antiferromagnetism, 44 which may enhance the magnetic frustration in the spin model.When the competition between ferromagnetism and antiferromagnetism reaches a balance in certain parameter regime, higher-order spin interactions will start to dominate the spin model.For example, there are two terms attached with χ IJK that contribute to the three-spin ring exchange interaction, Both terms give rise to new channels contributing to the chiral spin interaction S I • (S J × S K ).Compared to the usual t 3 /U 2 chiral spin interaction from the 3rd order superexchange channel, 45,46 these nonorthogonality enabled exchange channels could provide stronger chiral spin interaction in the strong coupling (large U ) limit, in favor of the chiral spin liquid ground state. 47e can repeat the analysis for the second-order perturbative correction H (2) .It can also be approximated by the connected part H (2) c as argued previously.If we look at the Heisenberg interaction and the chiral spin interaction, To the leading order in g, we get The Heisenberg term J IJ contains contributions from the standard superexchange channel.The chiral spin term K (2) IJK contains contributions from a new channel as nonorthogonal ring superexchange.
Collecting all contributions from Eq. ( 28), Eq. ( 29) and Eq. ( 31), there are three channels that contribute to the Heisenberg interaction -the gt term from the nonorthogonal exchange, the U term from the conventional exchange, and the t 2 /U term from the superexchange; and there are four channels that contribute to the chiral spin interaction -the g 2 t and gU terms from the nonorthogonal ring exchange in Eq. ( 29), the gt 2 /U term from second-order perturbation theory, and the conventional t 3 /U 2 term from third-order perturbation theory (which will appear in H (3) c ).In the strong coupling (large U ) limit, the novel channels originated from the orbital nonorthogonality typically dominate over the conventional superexchange and ring exchange channels.They are crucial to the analysis of the magnetism in Wannierobstructed Mott insulators.

IV. CONSTRUCTING LOCALIZED ORBITALS
In previous discussions, we have established the lowenergy effective spin model starting from the assumption of electrons localized on a set of nonorthogonal orbitals φ I (i).Now we come back to discuss why such arrangement is favorable and how these orbitals should be determined.In correlated materials, when the interaction energy dominates over the band width U W , it becomes energetically favorable to recombine single-particle states in the energy band to form localized orbitals, and to arrange one electron in each localized orbital to reduce the repulsive interaction.Although the Wannier obstruction prevents us from constructing orthogonal Wannier orbitals, it does not prevent us from constructing nonorthogonal and localized orbitals, on which electrons can reside.The criterion is to minimize the total energy of the system.Thus we start from the trial many-body state |Ψ σ proposed in Eq. ( 6), and minimize its energy Ψ σ |H|Ψ σ so as to optimize the localized orbitals φ I (i) that were used to construct the trial state |Ψ σ .
However, the energy Ψ σ |H|Ψ σ still depends on the spin configurations σ, which makes the objective function undetermined.To proceed, we focus on the "spinindependent part" of the energy, which is naturally the constant piece in the effective spin model H c .It corresponds to the Hartree energy H () in front of the identity operator χ () , as given in Eq. ( 13), Here n (i∇) denotes a real space representation of the band structure n (k).We will assume that the nth band n is well separated from other bands by a large band gap ∆, which is much larger than the interaction strength U .The separating energy scales (∆ U ) allows us to focus on the nth band only and find the optimal orbitals that can minimize the energy H () .9][50] The energy optimization δH () /δφ I = 0 boils down to solving the following Gross-Pitaevskii (GP) equation, The ground state can be found by imaginary time evolution of φ I (i) under the GP Hamiltonian (i.e. the operator on the left-hand-side of Eq. ( 33)).In each step of the evolution, the updated φ I (i) orbital will be broadcasted to other unit cells by translation.During the evolution, we do not impose the orthogonality among Wannier orbitals, so the orbitals we eventually obtain are in general nonorthogonal.The interaction will automatically determine whether or not the optimal orbitals will spontaneously break the point group symmetry.
We make a remark on the completeness of the localized orbitals.As is well known, for a single Chern band (suppose it is isolated for simplicity), it is impossible to choose a set of Bloch states |ψ k such that it is normalized and smooth in the Brillouin zone torus. 33This implies a Wannier obstruction to construct a set of exponentially localized orbitals which are orthonormal, complete and related to each other by translations.What if we relax the orthonormality constraint?Unfortunately, it is still impossible to find exponential localized orbitals which are complete and translation invariant, even if they are allowed to be nonorthogonal.If such orbitals exist, we can Fourier transform to obtain a set of unnormalized but smooth and nowhere vanishing Bloch states.Further normalizing these states leads to normalized and smooth Bloch states, causing a contradiction.It turns out that, under certain assumption which holds for the example we will be considering, it is possible to find a set of exponentially localized orbitals which are complete but nonorthogonal and also breaks the translation symmetry for one orbital.In other words, N − 1 number of orbitals are related to each other by translations but the last orbital takes a different form, where N is the total number of orbitals.We proved this claim in Appendix C. The orbital obtained from the energy optimization procedure described above, if exponentially localized, can be used to generate these N − 1 translation related orbitals, and a different last orbital is needed to complete a basis.We expect a single orbital to have little effect on the overall physics of the system, thus we ignore this subtlety hereafter.In the appendix, we also show that the dual orbital formalism, which has been useful in the study of Hubbard model ferromagnetism in topologically trivial bands 8,9 , is not applicable to a Chern band with the orbitals we constructed.This is one of the reasons that we take a different approach in this work.

V. APPLICATION A. Kagome Lattice Model and Wannier Obstructions
In this section, we apply our theory to an interacting fermion model on a Kagome lattice, whose Hamiltonian H = H t + H U takes the form of Eq. ( 1).The hopping Hamiltonian H t consists of purely imaginary nearest neighbor hopping only, with t ij = i/2 for j → i along the bond direction as specified in Fig. 2  The middle band is strictly flat with zero Chern number C = 0, which is an obstructed trivial band.Its Wannier obstruction is only due to the lack of lattice sites at the hexagon center Wyckoff position, such that the obstruction can be lifted by adding empty sites.The top and bottom bands are Chern bands of Chern number C = ±1 respectively.They combined together to form a set of fragile topological bands, which is Wannier obstructed and does not admit an effective tight-binding model description (despite their total Chern number being zero), 28,51 because the symmetry representations at high-symmetry momentum points do not match those of any atomic insulators of the same lattice symmetry.This situation is analogous to the middle bands in the tBLG around the charge neutrality, which have a fragile topological band structure for each valley. 27,32Placing the tBLG on the aligned hexagonal boron nitride (hBN) substrate further opens up the band gap at charge neutrality.The top and bottom bands are valley Chern bands of opposite Chern numbers, 12 which resembles the Chern bands in the Kagome lattice model as in Fig. 2(b).In the following, we will first set aside the possible connection to tBLG systems and focus on the Kagome lattice model itself.By applying our proposed approach to analyze this toy model, we wish to gain general understanding about Wannier obstructed Mott insulators, which could facilitate future study of correlated insulating phases in Moiré superlattice systems.

B. Nonorthogonal Localized Orbitals in Wannier
Obstructed Bands We follow the method described in Sec.IV to construct the localized but nonorthogonal Wannier orbitals.If we isolate the middle flat band and apply on-site repulsive interaction U ij = U 0 δ ij , the orbital φ I (i) that minimizes the energy H () defined in Eq. ( 32) is found to be strictly localized around the hexagon as shown in Fig. 3(a).This orbital is similar symmetry-wise to the localized orbitals proposed in Ref.52 for a different model, where the strict localization in both cases comes from the destructive interference between orbitals.
If we focus on the bottom Chern band with C = −1 and again apply the on-site repulsion U 0 , we found two degenerated solutions of φ I (i) as shown in Fig. 3(c) and  (d).These orbitals have similar features as the Wannier orbitals in twisted bilayer graphene (tBLG). 29,53They have three major peaks around the triangular lattice site and are equipped with non-zero angular momentum.To make connection to the tBLG system, we follow the notation in Ref. 31 to label these orbitals by the Wyckoff positions of their orbital centers and the angular momenta with respect to their orbital centers.The hexagon and the triangle center Wyckoff positions are denoted by a and b respectively as in Fig. 2(a), and angular momentum 0, ±1, ±2 are denoted by s, p ± and d ± .The orbital types are summarized in the table of Fig. 3.In the Mott limit, electrons will spontaneously choose one of the (b, p − ) orbitals in Fig. 3(c) and (d) to reside, which spontaneously breaks the C 6 rotation symmetry to the C 3 subgroup and results in the C 3 nematic phases.However, if we include longer range interactions, the energetically most favorable orbital can be different.For example, if we apply 3rd neighbor repulsion U ij = U 3 to the bottom Chern band, we obtain a (a, d + ) orbital as shown in Fig. 3(b).This orbital preserves the C 6 rotation symmetry, so the corresponding Mott phases are not nematic.As shown in Fig. 4, all the orbitals we constructed are indeed exponentially localized.As a result, the tensors g IJ , t IJ and U IJKL should all decay exponentially with the inter-orbital distance, so we expect the resulting effective spin model to exhibit well controlled locality.

C. Effective Spin Models and Possible Phases
In the Mott limit, the spin dynamics for any of these orbitals can be described by the effective spin Hamiltonian on a triangular lattice, as shown in Fig. 5.The spin-rotational symmetric Hamiltonian takes the following form where three sites IJK surrounding both up and down triangular plaquettes are arranged in the counterclockwise order.Again we focus on the bottom Chern band.The coupling strengths of the nearest-neighbor Heisenberg interaction J 1 and the chiral spin interactions K / are plotted in Fig. 6 as a function of t/U 0 or t/U 3 .Now we briefly discuss symmetries of this effective Hamiltonian and direct readers to Appendix D for more details.The electron Hamiltonian has translation symmetries T 1 and T 2 along two different directions, six-fold rotation symmetry C 6 , and an anti-unitary mirror symmetry M. The optimal localized orbitals φ I can spontaneously break some of the symmetries.For example, the (b, p − ) orbital breaks C 6 to C 3 (see Fig. 5).Given the symmetries of the orbitals, we can infer the symmetry of tensors g IJ , t IJ and U IJKL .It turns out that between nearest neighboring sites I and J, t IJ (g IJ ) can be generated by a single parameter t ≡ |t IJ | (g ≡ |g IJ |) given T 1,2 , C 3 and M. Then the chiral spin interaction K /K result from t IJ and g IJ must be opposite on neighboring triangular plaquettes, i.e.K = −K .The interaction tensor U IJKL breaks this pattern, but K /K result from U IJKL is much smaller than others in this model, so we still have a good approximate symmetry K −K .If we further have C 6 symmetry like in the case of (a, d + ) orbital, t IJ and g IJ will be restricted to real numbers and furthermore the chiral spin interaction is restricted to be uniform K = K .This symmetry analysis is in agreement with a previous study on a similar model. 17r the (b, p − ) orbital favored by the onsite interaction U 0 , the Heisenberg interaction changes from ferromagnetic (FM) to anti-ferromagnetic (AFM) as t increases (Fig. 6 (a)).In this case, the nonorthogonality enabled channel t IJ g JI favoring AFM.Thus, AFM interaction starts to dominate after this new channel takes over at large t.The same physics happens for the (a, d + ) orbital favored by the 3rd neighbor interaction U 3 , while the transition occurs at a much smaller t (Fig. 6 (b)).Close to the transition, there can be intermediate phases, which requires to take the chiral spin interaction into account.It turns out that near the FM-AFM transition regime, the chiral spin interaction is dominated by the tg 2 channel in Eq. ( 29) (see Appendix D).The (a, d + ) orbital respects the C 6 symmetry, which constrains the chiral spin interaction to be vanishing small (Fig. 6 (b)).On the other hand, the (b, p − ) orbital breaks the C 6 symmetry, so K and K are approximately related by the staggered pattern (Fig. 6 (a)).FIG. 6.The coupling strengths of the nearest-neighbor Heisenberg interaction J1 and the chiral spin interaction K , K for (a) the (b, p−) orbital in Fig. 3(c) favored by the U0 interaction, and for (b) the (a, d+) orbital favored by the U3 interaction (where K = K = K).FM/AFM stands for the (anti-)ferromagnetic phase; canted stands for the canted 120 • AFM configuration shown in the inset of Fig. 7. Now we compute the coupling strengths perturbatively with two control parameters t and g to get insight for their behavior in more general orbitals.We approximate the orbitals by only major and secondary peaks, and the orbital wave function is determined by the angular momentum and the amplitude ratio between major and secondary peaks tunable by g.Most importantly, the nearest orbital overlap g parameterizes the nonorthogonality of the orbitals, and controls the competition between different magnetic phases.Since g only slightly depends on the magnitude of U , we approximate it by a constant in the following analysis.For the (b, p − ) orbital, to the leading order in the hopping t and the orbital overlap g, we have where new channels in the theory of nonorthogonal basis contribute to those terms that contain g.The Heisenberg coupling J 1 changes sign at (U 0 /t) * = 1/g as shown in Fig. 6(a).In the flat band limit U 0 /t 1/g, the FM exchange dominates.In Appendix E, we provide a rigorous proof of the ferromagnetism for Wannier obstructed Mott insulators in the flat band limit.As the band dispersion gets larger (but still on the strong coupling side) 1/g U 0 /t 1, the AFM interaction takes over.Due to the geometric frustration on the triangular lattice (especially when higher order spin interactions are also taken into account), several candidate orders may compete for the ground state, which we will leave for later discussion.
From the above analysis, we expect the FM phase to become unstable toward AFM-like phases around U 0 /t 1/g.However, around this point, the exchange interaction J 1 tends to vanish, so we need to consider higher order interactions, e.g. the chiral spin interaction, Again, terms that contain g arise from nonorthogonality enabled channels.The chiral spin interaction dominates the spin model |K , | > |J 1 | in a narrow window of |U 0 /t − 1/g| < 13/8.A similar analysis can be repeated for the (a, d + ) orbital under third neighbor repulsion U 3 .
Combining these information, we get a schematic phase diagram in Fig. 7.
As a reminder, the coefficients in Eq. ( 35) and Eq. ( 36) are specific to the localized orbital in our model.However, we expect the Heisenberg interaction J 1 to be quadratic in t/U and g for generic systems, with different order-one coefficients, such that the FM-AFM cross over happens at t/U ∼ 1/g scale.A similar analysis can be generalized to the chiral spin interaction which is cubic in t/U and g.The spin model in Eq. ( 34) can give rise to different phases.We first consider the case when the on-site interaction U 0 dominates, which favors the (b, p − ) orbital.The (b, p − ) orbital breaks the C 6 symmetry to C 3 spontaneously, resulting in nematic phases.In this case, the chiral spin interaction is approximately staggered K −K .We present a classical picture of possible phases in the following.Away from the U 0 /t ∼ 1/g transition regime, the nearest Heisenberg interaction J 1 dominates, which leads to a Heisenberg FM state when J 1 < 0 and a 120 • AFM state when J 1 > 0. At the classical level, the in-plane 120 • AFM state can be tuned towards the z-axis FM state by canting the 120 • spin configuration in the x-y plane toward the z-axis, as illustrated in the inset of Fig. 7.The canted 120 • AFM configuration happens to have opposite spin chirality between up and down triangles, which is indeed favored by the staggered chiral spin interaction.We denote this intermediate state as the canted state in Fig. 6 and Fig. 7, which carries both non-zero magnetization and staggered scalar spin chirality.This spin configuration has been previously studied as an umbrella-type noncoplanar phase in Ref.54 within a Kondo system under different settings.Since both the canted state and the 120 • AFM state spontaneously break the spin U(1) symmetry, they actually belong to the same phase.In contrast, there has to be a phase transition between the canted AFM phase and the Heisenberg FM phase which preserves the spin U(1) symmetry (Fig. 7).However, this classical picture might be modified under quantum fluctuations and longer-range geometric frustrations.We then comment on the other case when the longer-range interaction becomes important.For example, the third-neighboring interaction U 3 favors the (a, d + ) orbital, which preserves all the lattice symmetries.In this case, the chiral spin term is parametrically small.The sign change of J 1 still happens when U 3 /t is around the order of 1/g, which drives the transition between FM and AFM phases.Such transition is likely first order (Fig. 7).Finally, when U/t ∼ 1, the system is no longer captured by the strong coupling theory presented in this work.We simply denote the weak coupling phase as the Fermi liquid (FL) phase, whose instability should be further analyzed using weak coupling approaches.
Let us further remark on some previous studies on the spin model in Eq. ( 34) regarding more exotic phases due to quantum fluctuation.Though it is well believed that uniform flux K K can drive the system toward a chiral spin liquid (CSL) phase 55,56 , it is not clear what happens in the staggered limit K −K .Some numerical study 57 and parton construction 58 on related models suggest a possibility of gapless spin liquid.When we further add a next-nearest-neighbor AFM coupling J 2 to the model [59][60][61][62] , it can drive the system toward a Dirac spin liquid (DSL) phase before the system fully develops a stripe order.These are all possible phases of the effective spin model we construct through nonorthogonal projection and perturbation.We will leave these rich possibilities for future numerical investigations.

VI. CONCLUSION
In this work, we present a different approach to understand Mott physics in Wannier obstructed systems including Chern insulators and fragile topological systems like twisted bilayer graphene.To get around the obstruction, we sacrifice the orthogonality of Wannier basis and develop a method to construct the trial nonorthogonal Wannier orbitals by numerically optimizing the Hartree energy of the system.In the Mott limit, we fill these trial orbitals with one electron per orbital.To study the low energy spin dynamics, we systematically project the Hamiltonian to a nonorthogonal spin basis and further study perturbative corrections.This new procedure concerning nonorthogonal Wannier basis gives rise to new channels to spin interactions.For example, at the level of direct projection, we find new channels that contribute to ferromagnetic, antiferromagnetic, or spin liquid phases.We demonstrate our approach with a toy model that carries Chern bands and fragile topological bands.In this model, new channels widen the antiferromagnetic phase and enhance chiral spin interactions that may lead to rich magnetic phases.
Our result may shed light on the magnetisms in Moiré superlattice systems, which often host Wannier obstructed bands.For example, the two middle bands near the charge neutrality in twisted bilayer graphene (tBLG) are identified to be fragile topological bands. 19,31,32ligning the tBLG with hexagonal boron nitride (hBN) substrate, the fragile topological bands further develop into separate Chern bands within each valley, which is analogous to the Chern bands in our toy model.The observation of ferromagnetic hysteresis 3 suggests that the three-quarter-filling insulating state in such system exhibits ferromagnetism.As the tBLG band structure only preserves the C 3 rotation symmetry within each valley, the scenario is similar to the U 0 dominated case in our toy model, where the localized (b, p − ) orbital exhibits the famous fidget spinner structure. 19,29,53r analysis shows a non-vanishing inter-site ferromagnetic coupling from the Fock term due to finite overlap g between nonorthognal orbitals even in the limit when the band width W approaches zero.This ferromagnetism becomes unstable when the band width increases up to U g set by the interaction strength U and the orbital nonorthogonality g, which is in consistent with previous studies on narrow Chern bands. 13When the system is close to the ferromagnet-antiferromagnet crossover, our approach provides a systematic framework to write down an explicit effective spin model that enables further numerical investigation of intermediate phases.Meanwhile, new channels from our framework also open up the possibility of new and richer magnetic phases close to the crossover.
intersection P 0 ∩ P 1 = ∅.Therefore, we would like to claim that H c dominates over others.To this end, we perform a numerical test of it on a finite-size system.We define the deviation between G −1 H and H c as The relation between D and typical overlapping weight g on a finite-size triangular lattice is plotted in Fig. 8.We artificially rescale the weight of each orbital g = αg , which is equivalent to transferring weights onto the Wyckoff position a and b.They do not belong to the kagome lattice thus have no contribution on permutations.We find that D(α = 1) < 0.01, meaning a high accuracy of approximating G −1 H with H c .In conclusion, within numeric capability, H c is a promising starting point to investigate the physics of the model.
Next, we want to address the issue related to the convergence of the whole series.We take the FM phase to examine this convergence.The energy of FM phase H c FM is calculated using χ FM = 1.Generically, there will be a g l total weight in front of each term with l = |P| being the length of the connected permutation.When l increases, the number of graphs also grows, but we have little knowledge of the speed.We need to determine which factor will dominate in the thermodynamical limit.To this end, we first transform the energy into the language of random walk.Since each term only encounters one H P , and each individual H P only differs from G P by a local factor t P(I)I or U P(I)IP(J)J , we are allowed to regard them as g IJ correspondingly without changing the convergence of the series.Then we have where k i is the overcounting factor for each intersection (e.g.k i = 2 for site i being visited twice).This is because whenever there is an intersection, simple random walk (RW) will have multiple choice to go through, which overcounts the H c FM .In fact, for backtracking process , there should be no discounting factor since there is only one way for RW to act.But neglecting backtracking does not influence the results of RW much, especially when the number of neighbors is large.Generally, the random walk problem can be written as f (g, p) = N 2 l l(−g) l p nc 0|RW(l)|0 , (B5) Depending on the optimal localized orbital achieved by minimizing the energy H () , part of the lattice symmetry can be spontaneously broken.For example, the (b, p − ) orbital in Fig. 3 breaks the six-fold rotation symmetry C 6 to three-fold C 3 , because there are two Wyckoff positions b in each unit cell and choosing one of them to occupy will necessary break the lattice symmetry.On the other hand the (a, d + ) orbital respects the C 6 symmetry, as it transforms under the C 6 symmetry as an irreducible representation.Under symmetry action, the localized orbital φ I transforms as Here G(I) denotes the new orbital index that I transforms to under the symmetry group element G.To preserve the translation symmetry, the orbitals follows the arrangement of the unit cells, and form a triangular lattice.The orbital (or unit cell) index I can be considered as the site index on the triangular lattice.In the Mott state, the electron spin degrees of freedom will ret on these sites.Given the symmetry transformations of the orbitals φ I , we can infer the symmetry transformations of the tensors g IJ , t IJ and U IJKL , which are defined via Note that t IJ has identical symmetry property as g IJ .We can use symmetry transformations to bring one tensor element at a particular link or plaquette to elsewhere through out the triangular lattice.For example, between the nearest neighboring sites I, J, required by the symmetries T 1 , T 2 , C 3 , M, g IJ = g * JI = g, t IJ = t * JI = t, (D6) noncoplanar phase by canting the 120 o Neel configuration in the x-y plane toward z direction, from which the chirality is automatically staggered.Further decreasing the interaction strength goes beyond our strong coupling framework, and the electrons fail to arrange into localized orbitals under weak repel interaction.In the main text, we present the localized orbitals of the Chern band.Now, we show the ones for fragile topological bands (combining top and bottom bands).To fill both of the bands, we need to fill two electrons per unit cell.Then the orbital index I labels both the unit cell and the orbital in the unit cell.In the case of two electrons per unit cell, the strong repulsion between two electrons force them to form nematic orbitals localized on a single site as shown in (a).The model was introduced by Ref.28 to demonstrate fragile topological insulators.This lattice model preserves translation symmetry and six-fold rotation symmetry C 6 (about the hexagon center).The single particle spectrum consists of three bands fully gapped from each other as shown in Fig.2(b).

FIG. 2 .
FIG. 2. (a) Kagome lattice model with imaginary hopping.Bound directions are specified by arrows.The gray hexagon marks out the unit cell.The Wyckoff positions a and b are respectively the hexagon and triangle centers.(b) The band structure of the Kagome lattice model, with the band Chern number C labeled.The inset shows the Brillouin zone and high-symmetry momentum points.

FIG. 4 .
FIG.3.Localized orbitals found by minimizing objective energy H () (the constant piece of the effective spin Hamiltonian Hc).The phase of the orbital wave function is specified by the colorbar.

FIG. 5 .
FIG. 5.The triangular lattice formed by the (b, p−) orbitals (in the bottom Chern band) arranged on the Kagome lattice.The effective spin model contains the Heisenberg interaction J1 across the bonds and the chiral spin interaction K / around the up/down triangles.It respects the translation T1,2, three-fold rotation C3 and anti-unitary mirror M symmetries.

FIG. 7 .
FIG. 7. Schematic phase diagram.The inset shows the spin configuration of the canted 120 • AFM order.

FIG. 8 .
FIG.8.In (a), we plot D against rescaled weight g = αg for two different nonorthogonal Wannier orbitals (Fig.3) on a six-site system.For α = 1 we have the original configuration.In (b), we plot RG flow of Eq.(B8), in which red points indicate fixed points.

2 l
B3) where |.|represent the length of connected permutation and the factor |P| for H comes from the fact that we can choose any bond to be t P(I)I .The 0|W |0 means the total number of graphs starting from origin and ending at it under rule W .And W = 1-loop(|P|) means we can only take a closed connected permutations with no self-intersections.Applying (−) n (−) n i=0 Pi = (−) n (−) n i=0 (|Pi|−1) = (−) −1+ n i=0 |Pi| , we have H c FM = N

FIG. 9 .
FIG. 9. (a) Symmetries of the lattice model.Consider (b, p−) orbitals that breaks the C6 symmetry, the remaining symmetries fix the pattern of (b) gIJ (same as tIJ ) and (c) UIJJK .
D3)For unitary symmetries G = T 1 , T 2 , C 3 , C 6 , they transform asG : g IJ → g G(I)G(J) , t IJ → t G(I)G(J) , U IJKL → U G(I)G(J)G(K)G(L) ,(D4)such that tensor elements related by the symmetry should simply be equal to each other.Only anti-unitary symmetry M relates them by additional complex conjugateM : g IJ → g * M(I)M(J) , t IJ → t * M(I)M(J) , U IJKL → U * M(I)M(J)M(K)M(L) ,(D5)

FIG. 10 .
FIG.10.Phase transition for the stably topological (bottom) band with respect to on-site interaction U0 (a) and 3rd-neighbor interaction U3 (b).In (b), we multiply the CSO by 10 to fit the graph.The CSO is staggered for (a) (not shown here).

Fig. 11 .
FIG. 11.Nonorthogonal Wannier orbitals for fragile topological band (combining top and bottom bands).(a) and (b) are a possible set of orbitals in a single unit cell.