Interaction-driven topological and nematic phases on the Lieb lattice

We show that topological states are often developed in two dimensional semimetals with quadratic band crossing points (BCPs) by electron-electron interactions. To illustrate this, we construct a concrete model with the BCP on an extended Lieb lattice and investigate the interaction-driven topological instabilities. We find that the BCP is marginally unstable against infinitesimal repulsions. Depending on the interaction strengths, topological quantum anomalous/spin Hall, charge nematic, and nematic-spin-nematic phases, develop separately. Possible physical realizations of quadratic BCPs are provided.


I. INTRODUCTION
The search of new topological states of matter has never been stopped since the discovery of the quantum Hall state in 1980's. 1 In particular, in recent years the study of "topological insulators" (TI) becomes one of the most active fields in condensed-matter physics, not only for its importance to fundamental physics but also for its potential application in spintronics or thermoelectrics. 2,3 This new insulating phase is distinguished from the conventional one by a non-trivial Z 2 topological invariant and robust gapless edge states in two dimensions (2D) [4][5][6][7] or surface states in three dimensions (3D) 8,9 , against moderate perturbations which preserve time reversal symmetry (TRS).
Generally, such topological insulating state can occur in a system with strong spin-orbit coupling that explicitly breaks spin rotational symmetry (SRS), resulting in the band inversion phenomenon. Typical materials which exhibit TI phase are found in, for instance, the HgTe/CdTe quantum wells (2D), Bi x Sb 1−x , Bi 2 Se 3 , Bi 2 Te 3 (3D), and so on. [12][13][14][15][16] However, an alternative route to TI is possible and it is associated with the concept of topological Mott insulator, first studied in Ref. [11] in strongly correlated systems. There are two key and generic ingredients in this approach. First, the Fermi surface of the non-interacting system should shrink to discrete points (2D) or lines (3D), and hence it is semi-metallic; second, there exists a suitable repulsive interaction, which can dynamically generate spin-orbit coupling through spontaneously broken SRS, as first discovered by Wu and Zhang. 10 A few pioneering examples along this line of thought have been discussed in various lattice geome-tries, e.g., honeycomb, 11,17,18 checkerboard, 19 kagome, 20 and diamond lattices. Here we show that topological states can be generally realized in a system with quadratic band crossing points (BCPs), which are symmetry protected at noninteracting level. Near such kind of BCP, instability towards phases with broken symmetry symmetries is inevitable even if there is only weak interaction between electrons. To demonstrate this, we construct a concrete model with such a BCP on an extended 2D Lieb lattice. There are several reasons for us to choose the Lieb lattice. First, the Lieb lattice has three sites per unit cell, as shown in Fig. 1. With only nearest-neighbor (NN) hoppings, there is a dispersionless (flat) band in the middle of the band structure. There are three scenarios to alter this non-generic crossing between this flat band and the two other bands, depending on the values of the onsite potentials or the presence of further range hoppings: (i) The flat band touches upper and lower linearly dispersing bands at one point when all lattice sites have equal onsite potentials. (ii) The flat band can be isolated. For instance, one can add intrinsic spin-orbit couplings as discussed by Weeks and Franz. 23 (iii) When the onsite potential on A sublattice is not equal to that of B/C sublattices, the flat band only touches one of the other two bands, which becomes quadratically dispersing, instead of linearly dispersing. 24 Thus, these choices could variegate our results. Second, a nearly flat band can have non-zero Chern number (see Sec. IV). Finally, the 2D Lieb lattice has been the most important building block in many 3D perovskite materials featured with complex phase diagram and strongly electron-electron correlations. Thus study of the model can be viewed as a preliminary investigation of the TI phase especially in layered perovskites composed of weakly coupled 2D planes with Lieb lattice structure (e.g., the well-known high-T c cuprates).
In this paper, we start with the construction of the explicit model and reveal the topological nature of the BCP at non-interacting level. We then examine the consequence of such topological BCP, i.e., with a symmetry protected Berry flux 2π, under the presence of shortrange repulsive interactions. We investigate various symmetry breaking instabilities at BCP within self-consistent mean-field approximation. Note that we mainly focus on type (iii) band structure, namely, only two bands touch together, and compare it with the case of type (i) when necessary. In principle, for the BCP there are two ways to gain energy: One is to open a gap at BCP, and the other one is to split into two Dirac points (each with Berry flux π), but with the price of broken C 4 symmetry. To justify this speculation, we show phase diagrams for spinless/spinful fermions at 1/3 or 2/3 filling, according to the position of the BCP in the band structure. In fact, at both fillings the phase diagrams are qualitatively similar with subtle differences due to particle-hole asymmetry introduced by the interactions. In the spinless case, the leading order under "weak" next nearest-neighbor (NNN) repulsion is the quantum anomalous Hall (QAH) insulating state (TRS broken). For "strong" NNN repulsion, the ground state evolves into insulating nematic state (C 4 symmetry broken down to C 2 ). In addition, for intermediate strength, there exists a narrow coexistence region between these two orders. In the spinful case, the phase diagrams are more complicated. Besides the phases we find in the spinless case, there are also a quantum spin Hall (QSH) insulating state and a nematicspin-nematic semi-metallic phase with Dirac nodes. [25][26][27] Thus, we clearly demonstrate that, in principle, correlated systems with Lieb lattice structure can be a potential host to TIs. This paper is organized as follows. In Section II, we define the model and demonstrate the topological nature of its BCP from both momentum and real-space point of view. Next, the consequences of introducing shortrange repulsions are discussed for spinless fermions in Section III A and for spinful fermions in Section III B, respectively. Finally, we discuss some issues and make conclusions in Sec IV.

II. THE LATTICE MODEL
We begin with the tight-binding model for noninteracting fermions, where c † i creates a fermion on site i of the 2D Lieb lattice, the unit cell of which is given by A, B, C sites shown in Fig. 1. For simplicity, we take the hopping amplitudes, t ij = t, between NN sites ij , and t ij = 0 otherwise for the moment. The effect of adding longer-range hopping amplitudes (but small in magnitude) will be discussed later when appropriate. Note that the C 4 point group symmetry dictates that the onsite potentials on the B and C are equal, ε B = ε C . Although ε (A,B,C) are generically non-zero, only their relative values are essential to determine the symmetry of the lattice, and hence, the band structure. Therefore, hereafter we set the units of energy t ≡ 1, the lattice constant a ≡ 1, and, without loss of generality, ε B = ε C = 0.
An important feature of this model is that the presence or absence of ε A can change electronic properties dramatically. When ε A = 0, the flat band touches two linearly dispersing bands at M point in the FBZ [type (i)], where the linear bands meet as if there were a "Dirac point". However, the touching point in fact has completely different structure. It becomes clear once we expand H 0 (k) around M point with k = M + p, |p| ≪ 1. To the first order in p, H 0 (k) can be written as where the Fermi velocity v F = t, B = (p x , p y , 0), and the (pseudo) spin-1 matrices are defined as obeying Lie algebra of SU(2), i.e., [L i , L j ] = iǫ ijk L k , 33 instead of a Clifford algebra as in the case of graphene. This is the fundamental reason why there is no Dirac point and hence no fermion doubling problem 34 on the Lieb lattice. Viewing the low-energy effective Hamiltonian H 0 (B) as a "spin" L in an external "magnetic field" B, its eigenvalues can be easily read out as ǫ(B) = v F |p|l p , where l p = 0, ±1 are the quantized angular momenta along the axis parallel to B in three dimensions. When ε A = 0, however, the spin-1 structure mentioned above is no longer valid. The flat band now touches only one dispersive (massive) energy band either above or below at M point, depending on the sign of ε A [see Fig. 2 To make this structure transparent, we again expand H 0 (k) around M with small p. Assuming |p| ≪ |ε A /t| and ε A < 0 at 2/3 filling, we then integrate out the contribution from basis A (due to almost fully filled A sublattice) and obtain a low-energy effective two-band Hamiltonian, where m 0 = −ε A /v 2 F . In the last equality, we express H ef f 0 in terms of the identity and Pauli matri- Interestingly, if we further allow small third-neighbor hoppings t ′′ > 0 (but forbid to hop whenever there is a site in the middle of the path), the flat band becomes slightly dispersive and the effective Hamiltonian changes to without removing the band crossing point (BCP) [see Fig. 2 (b)]. Such point at p = 0 is the so-called quadratic band crossing point (QBCP), which has been studied recently by several research groups. 19,[28][29][30] One of the key features for QBCP in 2D is that its density of states (DOS) is non-zero at the crossing point, in sharp contrast to the case of Dirac points. This will lead to essential difference when responding to the weak interactions present in the system. In the following, we will mainly focus on the ε A = 0 case and show that the BCP in our model is not only topologically non-trivial, but also make the system be a potential host to topological phases under weak repulsive interactions. The band touching phenomenon on the Lieb lattice is quite generic and stable for non-interacting fermions. Such stability deserves a full analysis here. We shall provide two different approaches to show it: One is based on momentum-space topology, and the other one is based on real-space topology.
From the first point of view, the BCP actually forms a topological defect in the momentum space, similar to a vortex in a 2D superconductor, here with a winding number ±2. To see this, let us rewrite Eq. (6) as . This effective Hamiltonian now represents a spin-1/2 particle sitting in a magnetic field B ′ , which has a vortex structure at p = 0, as shown in Fig. 3(a). For comparison, recall that for ε A = 0 we instead have a spin-1 particle in an external field B [Eq. (4)], whose structure is shown in Fig. 3(b). The winding number W then can be easily extracted from the figures that in the former case, W = 2; in the latter case, W = 1. However, somewhat counter intuitively, both cases are associated with the same Berry phase of the BCP, 31 which can be calculated precisely by where Γ is a contour in the p space enclosing the touching point, n denotes any one of the involved bands, and |u np represents the Bloch wave function for nth band. A simple argument solves this puzzle. The line integral along any loop enclosing 0 in p space given above is known to be 1/2 (1) times the solid angle subtended by B ′ (p) [B(p)] from the origin for a spin 1/2 (1) particle. Thus, B n = 2πW/2 = 2π in the former case, which is just equivalent to B n = 2πW = 2π in the latter one. In fact, when B n = 0, any infinitesimal mixing (perturbation) between bands would lift the degeneracy. With non-vanishing B n = ±2π, we confirm that the BCP on the Lieb lattice at non-interacting level is topologically stable (i.e., not opening a gap) as long as the spinless system preserves both TRS and C 4 point group symmetry. Note that C 4 symmetry in our model is quite essential, as a similar QBCP happens in the A-B stacking bilayer graphene (with C 3 symmetry) while it can easily decay into Dirac BCPs and thus is topologically unstable. 19,22 An alternative point of view for the protection of such BCP on the Lieb lattice can be associated with certain topological structure present in the real space, or more specifically, with the existence of the eigenstates which are extended along non-contractible loops winding around the whole lattice with periodic boundary conditions (i.e., a torus). 32 To demonstrate this feature, we first take the merit of the flat band, which allows us to construct its corresponding localized, one-particle eigenstates of H 0 (Wannier states). Taking R to be the coordinate of the central site of the shaded plaquette shown in  , we find that the creation operator for the localized eigenstate at R can be written as The key reason for these states being localized is rooted on the fact that all A-sites have vanishing amplitudes and remain zero after the action of H 0 on them due to destructive interference. In localized-state language, the existence of the BCP in our model with ε A = 0 is equivalent to state that the dimension of the space expanded by independent localized eigenstates with zero energy has a dimension which is one larger than the number of unit cells, N . The extra state cannot come from the flat band, but from one of the dispersive bands. The plaquette states we constructed in Eq. (8) seem to form N linearly independent states with zero energy. For our model with periodic boundary conditions, however, the following relation, reduces the naive counting by one and hence only N − 1 states are independent. The missing two states, in fact, are accounted for by two non-contractible loops around the whole lattice (torus), as illustrated in Figs. 5(a) and (b). When H 0 acts on these states, the destructive interference again guarantees the zero eigenvalue. Now, in total we have N + 1 independent states. Therefore, provided not destroying the flat band, such band touching phenomenon is protected by the topological character of the lattice.

III. INTERACTION DRIVEN INSTABILITIES
The existence of such symmetry/topology protected BCP on the Lieb lattice at the non-interacting level motivates us to further ask if it is stable in the presence of repulsive interactions. To see this, we will first examine whether generic short-range repulsions are relevant to this BCP from perturbative renormalization group (RG) analysis, and next, if the interactions are relevant, we will investigate possible consequences of such instability, i.e. symmetry breaking phases, at mean-field (MF) level.
To perform RG analysis, we consider a continuum, spin-1/2 Hamiltonian, which can be obtained by projecting the original three-band interacting model onto an effective two-band theory near the BCP in the continuum and |ε A | ≫ t limit, where the free part H 0 reads as The subscript σ denotes spin polarization and the fermion field Ψ † σ = (ψ † 1σ , ψ † 2σ ) with 1,2 representing orbital (i.e., two touching bands) degrees of freedom. Note that in the momentum space the expression for H ef f 0 is given by Eq. (6), and is independent of σ. For simplicity, we take small t ′′ = −t 2 /ε A such that d I = 0, making the effective theory particle-hole symmetric. In fact, non-vanishing d I would not change our main conclusion, provided m −1 0 > t ′′ . The projected interacting part includes only 1) intra-orbital and 2) inter-orbital contact interactions, , where u and g σσ ′ are intra-orbital and inter-orbital coupling parameters, respectively. For the chemical potential µ = 0, the non-interacting system H ef f 0 leads to one Fermi point at p = 0 with nonvanishing DOS, instead of a Fermi surface. Setting the dimension [p] = 1 and understanding the dynamical critical exponent z = 2 due to quadratic dispersion, it is straightforward to see that the dimension [ψ aσ (r)] = 1; in the interacting part the coupling constants [u] = [g σσ ′ ] = 0, implying that they are superficially marginal interactions. However, as shown in Appendix A, we find that they are generically marginally relevant and bring the system to the strong coupling regime. More explicitly, up to oneloop order, the coupled RG equations for the coupling parameters are where we have used the fact that g ↑↑ = g ↓↓ and g ↑↓ = g ↓↑ . l denotes the momentum rescaling p → pe −l and the coefficients α = 1 In fact, no new fixed point (FP) is produced in this set of RG equations, except for the non-interacting one at which u and g σσ ′ vanish. Moreover, we find that given generic bare coupling parameters (u, g σσ ′ > 0), at least one of them diverges first when reaching low enough energy scale. This indicates that the non-interacting FP is an unstable FP, which can drive the system to strong coupling regime in the presence of short-range repulsions. In addition, it is worth mentioning that by setting g ↑↓ = u = 0, we reduce the RG equation back to the spinless case, with α > 0, which is consistent with the work done in Ref. 19 and, importantly, it implies that short-range repulsions are again marginally relevant.

A. Spinless fermions
From above we know that a QBCP is generally unstable against weak repulsive interactions. We now discuss its consequence on the Lieb lattice and explore possible symmetry breaking phases driven by interactions at mean-field level. We warm up with the spinless case to gain some physical insights before including spin degrees of freedom.
The lattice model we study is given by Eq. (1), with short-range interacting terms, where V 1 and V 2 are repulsive coupling constants for NN and NNN interactions, respectively. n i = c † i c i is the number operator on the site i. The chemical potential is suitably chosen to keep the system at 2/3 (1/3) filling for ε A < 0 (ε A > 0). We proceed by treating H int in the MF approximation, including both the on-site and bond MF decoupling particle-hole channels, where φ ij = φ * ji = c † i c j represent certain current/bond order with i, j belonging to NN and NNN bonds. Note that in this work only translation-invariant MF ansatz is considered. The repulsive interactions can produce the following possible phases: (i) Nematic state. This is a phase associated with broken C 4 symmetry down to C 2 . It behaves like an anisotropic metal (one QBCP splits into two Dirac points) or an insulator (two Dirac points meet at zone boundary and end up with a gap), depending on the strength of repulsive interaction. There are two types within this phase. Type I is "site" nematic with order parameter, η = 1 , where δ ′ = ±x/2 ±ŷ/2 denoting four NNN bonds of B-site and A, B, C are sublattice indices. Without loss of generality, at 2/3 filling we can set the charge density, c † Ai c Ai = 2 3 + ρ, c † Bi c Bi = 2 3 − ρ 2 + η, and c † Ci c Ci = 2 3 − ρ 2 − η. As ε A > 0, at 1/3 filling we simply replace 2/3 by 1/3 in the above expression. Note that the use of the parameter ρ is to take into account the renormalization of the onsite potentials due to interactions. The nonzero expectation value of it does not break any symmetry of the model. On the other hand, type II is "bond" nematic with order parameter either in the form of This topological state is known to be characterized by quantized Hall conductance without Landau levels (or equivalently, by non-zero Chern number) and has topology-protected, gapless chiral edge modes. 35 We compute the Chern number for each band within this state and find that 1) for |ε A | > 0, the previous two touching bands now carry Chern numbers ±1 separately. In particular, one of two bands is (nearly) dispersionless. The third one simply carries zero; 2) for ε A = 0, the middle flat band carries zero Chern number, while the upper and lower bands carry ±1, respectively.
The other possible current loop states, however, are not topological insulating. For case (b) (Varma Θ I loop state 36 ), it is semi-metallic with order parameter given by, simultaneously, , it is insulating with broken inversion symmetry (IS) as well but is invariant under the combined TRS and IS. Thus, there is no Hall or uniform Kerr response by noticing that for any given momentum k it changes sign under TRS or IS. 29 This order can be described as, Finally, case (d) is a semi-metallic state with broken IS and hence no net Hall current in it. Its order is described as Φ ′′′ In momentum space, the mean-field Hamiltonian at 2/3 filling can be now written in the matrix form, with the fermion spinor, Ψ † k = (c † Ak , c † Bk , c † Ck ), and where δt (δt ′ ) represents a renormalization due to NN(NNN) repulsions. The Hamiltonian matrix H k reads , and ν x,y = −2t ′′ cos k x,y ; Γ parameters are given by Thus, the mean-field free energy can be expressed as where β = 1/k B T and E k are eigenvalues of H k . The ground state with given coupling parameters can then be determined by minimizing the free energy with respect to each order parameter, yielding a set of coupled gap equations. Notice that at 1/3 filling we follow the same procedure and only the diagonal part of H k and E 0 need to be changed accordingly due to the shift of the average charge density. We numerically solve the coupled gap equations selfconsistently and obtain the zero-temperature V 1 -V 2 phase diagrams at both 1/3 (ε A = 4t) and 2/3 (ε A = −4t) fillings, as seen in Figs. 7(a) and (b). Note that in this study, only weak short-range repulsions, i.e., V 1 , V 2 ≤ t are considered. At 1/3 filling, we find that there are three phases in the absence of V 1 : QAH phase, coexisting QAH+nematic phase, and nematic phase. Beginning with V 2 ≪ t, the infinitesimal instability of QBCP leads to QAH phase by the second order phase transition, with where N 0 denotes the finite DOS at QBCP and Λ is an energy cutoff; on the other hand, for t V 2 > V 2c ∼ 0.22t the ground state breaks C 4 symmetry spontaneously down to C 2 and exhibits insulating nematic phase with a gap ∆ N ∼ η. In this phase, we find that the site-nematic order (η) is the dominant one and a small component of the bond-nematic order (Q 1 ) accompanies with it. In fact, we notice that the bond-nematic order cannot be induced by V 1 itself. Finally, with intermediate value of V 2 , there exists a narrow window for the coexistence of both QAH and nematic orders. In addition, we also notice that there is no room for the current loop states other than QAH state for systems with relatively large |ε A | and weak repulsions. At 2/3 filling, the phase diagram is qualitatively similar to that at 1/3 filling. However, there are a few remarks worth mentioning here: 1) Although the noninteracting energy spectrum for both fillings can be related by translating "particle" into "hole" language, which causes t → −t, t ′′ → −t ′′ , ε A → −ε A , and µ → −µ. The interactions given in the present form ruin such relation and hence the two phase diagrams must be different. 37 2) The response to the weak V 1 repulsion is different for the two cases. At 1/3 filling with ε A ≫ V 1 > 0, since the charge density on A sublattice is much less than 1/3, the stronger V 1 is, the more energy cost the system suffers from V 2 . Thus, the region of the QAH phase reduces and the ground state evolves into the nematic phase faster than the case with V 1 = 0. In contrast, at 2/3 filling with |ε A | ≫ V 1 > 0, the presence of small V 1 can only reduce the charge density slightly on A sublattice, and therefore its effect on the energy cost from V 2 is relatively small. 3) When ε A , t ′′ → 0, the spin-1 structure near band touching point, as discussed in the previous section, is recovered. Our mean-field study shows that the infinitesimal instability (near 1/3 or 2/3 filling) is absent due to the vanishing DOS of the disper-sive bands and the semi-metallic phase is robust until V 2 reaches certain critical value. Moreover, for V 2 ≥ V 2c , we find that the QAH phase only survives in a negligible window of V 2 , and the nematic phase becomes the dominant one in the phase diagram (not shown). This result is similar to the work done by Q. Liu et al. 20 on the 2/3-filled kagome lattice with Dirac BCPs.

B. Spinful fermions
We now take the spin degrees of freedom into account. The model Hamiltonian again consists of the free and interacting parts, i.e, H = H 0 + H int . The free part is again given by Eq. (1) with extra spin index σ in the fermion creation/annihilation operators; the interacting terms now contain In addition to V 1 (V 2 ) denoting the coupling constant of NN (NNN) repulsion, we also consider the repulsive Hubbard (U i ) terms. For simplicity, we will assume uniform onsite repulsions, i.e., U i = U . Different from the spinless case, there are not only spin-singlet order parameters (as we had before), but also spin-triplet order parameters within mean-field approximation. The possible phases under translation-invariant ansatz are classified below: (i) Charge nematic state (CN). This phase is associated with spontaneously rotational (C 4 ) symmetry breaking. One can either have the site-nematic or bond-nematic state, whose order parameter is a spin-singlet and simply the same as that for spinless fermions with additional summation on spin σ times a normalization factor 1/2. Note that for the site-nematic case, the driving force is now from both U and V 2 terms, combined together to give out an effective NNN repulsion, V ′ 2 = 2(V 2 − U 8 ), playing similar role of V 2 in the spinless model.
(ii) Nematic-spin-nematic state (NSN). This phase breaks C 4 symmetry in the spin sector, not in the charge sector. Consequently, it turns the (spin) doubly degenerate QBCP into four Dirac points (two pairs with opposite spin polarizations), and C 4 symmetry of the band structure remains intact. Similar to its charge counterpart in (i), there are two types: One is the site-NSN with a spin-triplet order parameter, η t = 1 16 αi,σ s σσ ′ c αi,σ ′ with α = A, B, C and s, the Pauli matrices. Note that this phase can occur simply due to the presence of the Hubbard term, which provides a spin-triplet channel, − 2 3 U ( S i ) 2 with S i representing usual spin operator. The other one is the bond-NSN, which is described by Q t (iii) Charge current-loop state with broken TRS. As mentioned in the spinless model, among all the current patterns the case (a) in Fig. 6 is of most interest. This state, characterized by spontaneously broken TRS and parity symmetry, exhibits QAH effect with a spin-singlet order parameter, . Mainly driven by V 2 terms, fermions with opposite spin polarizations flow in the same way and hence provide the same flux pattern, penetrating the whole lattice. In addition, we investigate the possibility for the other currentloop states [e.g., from Fig. 6(b)-(d)] and find that none of them is stabilized by the presence of short-range repulsions in our MF study. Therefore, we will only consider case (a) hereafter.
(iv) Spin current-loop state with TRS. The key difference of this phase from its charge counterpart (QAH) is TRS unbroken. One can view it as a combined doublelayer QAH system: Fermions of opposite spin polarizations, residing in different layers, producing just opposite flux patterns separately. Thus, for the whole system TRS is preserved. In fact, this is known as quantum spin Hall (QSH) phase with spin-triplet order parameter described by Φ t . Two remarks deserve mentioning here. First, both QAH and QSH are topological phases, characterized by nontrivial topological index with robust edge states. However, the former one acquires non-vanishing Chern number, while the latter one has zero Chern number due to TRS. Thus, a new topological index, called Z 2 index ν, needs to be introduced. 4,5,9 As detailed in Appendix B, the QSH phase on the Lieb lattice indeed acquires non-trivial ν = 1. Second, it is straightforward to see that at MF level, the energy spectra for both QAH and QSH are the same. As a result, they have equal energy gain from V 2 repulsion and hence one cannot distinguish them in the MF phase diagram. If there were an extra NNN exchange coupling J 2 present in the system, the QAH would be favored for J 2 > 0; reversely, the QSH would be favored for J 2 < 0 due to its spin-triplet nature.
Under our assumption of translational invariance within our mean-field study, we do not consider any charge or spin density wave order. However, it is still worth pointing out that if |ε A | → ∞ and hence makes the A sublattice be effectively decoupled from rest of the lattice sites, at large U the (0,0) antiferromagnetic order could be realized at 1/3 (2/3) filling with ε A > 0 (ε A < 0), as guaranteed by the Lieb's theorem. 38 Following the same procedure as in the spinless case, we decouple the interacting terms within MF approximation and obtain the MF free energy in similar form of Eq. (22), where we have used the MF ansatz, n Ai,σ = 2 3 + ρ + 1 2 σS Az , n Bi,σ = 2 3 − ρ 2 + η + 1 2 σ(S Bz + 2η t z ), and n Ci,σ = 2 3 − ρ 2 − η + 1 2 σ(S Bz − 2η t z ) for the 2/3-filled lattice. Note that, for simplicity, we have assumed that spin points to z-direction after SU(2) symmetry breaking. Such treatment will be applied to other spin-triplet order parameters as well. E 0 in the spinful case becomes whereε A1 ,ε B1 ,ε B2 ,ε C1 , andε C2 are defined as before. E k are eigenvalues of the MF Hamiltonian, which is now a 6 × 6 matrix. By minimizing free energy with respect to various order parameters, the T = 0 mean-field U -V 2 phase diagram at 2/3 filling (ε A = −4t) for spinful fermions is shown in Fig. 8. We first notice that in the absence of U there again exist three phases: QAH/QSH phase, coexisting QAH/QSH+CN phase, and CN phase from weak to strong V 2 repulsion. In particular, the fact that topological QAH/QSH phase can arise from infinitesimal instability of QBCP further justifies the interaction-driven scenario as a promising way for producing TI. In the presence of U , however, NSN phase begins to compete with the topological phase and clearly dominates over QAH/QSH whenever U ≫ V 2 > 0. On the other hand, as V ′ 2 = 2(V 2 − U 8 ) 0.22t, the insulating CN phase takes over the phase diagram and this is consistent with the result shown in the spinless model. Two remarks are worth mentioning here. The first one is about the effect of NN repulsion V 1 . As in the spinless case, at 2/3 filling the phase diagram is not sensitive to the presence of V 1 . Especially, V 1 itself does not lead to any order. How-ever, when V 1 becomes stronger, it is quite possible that the system might gain certain energy by opening a gap due to translation symmetry breaking, and hence beyond our current consideration. The second point is about the system at 1/3 filling with ε A = 4t. In fact, the phase diagram is qualitatively similar to that of 2/3 filling and therefore we omit it without further discussion.

IV. DISCUSSION AND CONCLUSION
The model we have solved on the Lieb lattice demonstrates that the TI (QAH/QSH state) can be induced by appropriate interactions through spontaneously symmetry breaking mechanism, which dynamically generates spin-orbit couplings necessary for a topological insulator. It is then natural to ask how it can be realized experimentally.
We notice that in the interaction-driven scenario there are two key conditions that a system has to fulfill: 1) The band structure should contain suitable band crossing point at which two touching bands have opposite curvatures. 2) The system should have weak (or no) spin-orbit coupling [spin SU(2) symmetry is preserved] and NNN repulsions need to be more significant than the other shortrange repulsions. While the condition 2) is tricky and we have to reserve it for future investigation, we would like to comment on some possible routes for condition 1) below.
First of all, the most promising candidate, we believe, is from cold atom system. As discussed by Goldman et al. in Ref. 33, the Lieb lattice may be constructed as an optical lattice created by properly arranged laser beams. In particular, the spin-orbit interaction, which might be an issue in traditional materials, now becomes irrelevant. Another potential way for realizing Lieb lattice may come from layered perovskites. A well-known example is the CuO 2 plane in high-T c cuprates such as La 1−x Sr x CuO 4 or YBa 2 Cu 3 O 7 , whose electronic structure might be captured by a three-band model with p x , p y and d x 2 −y 2 orbitals. Here, for illustration purpose we take a typical three-band (Emery) model, written in hole language, having the same form of Eq. (1). 39,40 The corresponding model parameters are given by Hybertsen et al. 42 : t = 1.5, t ′ = 0.65, ε A = 0, and ε B,C = 3.6, where t ′ denotes NNN hopping. The band structure in the FBZ along high symmetry lines is presented in Fig. 9(a). As one can see, there is indeed a BCP at M, but with "wrong" curvatures for two touching bands. As a result, cuprates may be impractical to produce topological phase as we desire. To overcome this issue, we offer two speculative suggestions. The first way out is to add decorated elements with suitable orbital nature (e.g., p-orbital) between atoms on B and C sites that could change the sign of t ′ . This change leads to the band structure shown in Fig. 9(b) with "right" curvatures now. Keeping the same orbital characters, the second way is to search for a new system among perovskite materials, whose model Hamil-tonian is similar to that given by Sun et al.'s recent work in Appendix F. 43 The key feature of such three-band model is that the dominating transfer integrals are now associated with distance a after hopping. For instance, if, by certain geometric reason, the relevant orbital on B (C) sublattice becomes p y (p x ), instead of p x (p y ) as in the cuprates, the above consideration could be plausible. Before concluding our work, we would like make a brief remark for the possible new physics brought by the (nearly) flat band appearing in the Lieb lattice. Consider ε A > 0 and only NN hoppings in our spinless model. Now, the middle band is completely flat and touches with the lower dispersive band at M in FBZ. At 1/3 filling if we turn on V 2 , the ground state would enter QAH state with an energy gap ∆ QAH opened at M. As we have mentioned in Sec. III A, the middle and lower bands acquire Chern numbers ±1, while the upper one has zero Chern number in this phase. Moreover, the flatness of the middle band is approximately determined by the ratio of band-gap to bandwidth ∼ O(ǫ A /∆ QAH ). In other words, with appropriately chosen ǫ A /∆ QAH , one can produce a nearly flat band with non-trivial Chern number. This result is quite significant since it may provide an opportunity to realize fractional QAH state (or fractional Z 2 TI in the spinful case) when such band is partially filled. [44][45][46][47][48][49] In conclusion, we have studied interacting spinless/spinful fermions on the (extended) Lieb lattice and have explored the possibilities of various spontaneously broken symmetries associated with a BCP in the band structure. Due to the topological nature of the BCP, namely, with Berry phase 2π, we have seen that in the T = 0 phase diagram the system can exhibit topological QAH/QSH, nematic, and nematic-spin-nematic phases, depending on the strengths of short-range density-density interactions. In particular, the existence of TI phase firmly justifies the interaction-driven scenario. Moreover, for a quadratic BCP (as ǫ A , t ′′ = 0), only weak interaction is necessary for inducing the TI phase, which is in sharp contrast to the systems with Dirac points. In addition, in our model there exists a flat band, which is interesting on its own right and we argue that in principle, one can obtain a nearly flat topological band without external magnetic field, a starting point to realize exotic correlated phases of matter. Still, there are many open issues and they deserve further investigation. For instance, one could consider the effect of the chemical potential away from the BCP or the instability to superconductivity from both repulsions and attractions, with special focus on the possibility of any topological nature. Finally, we hope our work can further stimulate people to search for a new family of TIs in the perovskite-related materials, where, due to the presence of transition elements, the strong correlations may not be ignored.
In this section, we derive RG equations given in Eq. (13) in path integral formulation and show the RG flows for two typical cases. We begin with defining the action and the shorthand notations below. At zero temperature, the full action is given by where the free action S 0 reads in momentum space. The Grassmann variables,Ψ σ = (ψ 1σ ,ψ 2σ ) with 1,2 (labeled by 'a' for later use) representing orbital degrees of freedom. As mentioned in the main text, for simplicity, we consider the inverse of the non-interacting Green's function [Ĝ σ 0 (ω, p)] −1 = iωI − d x σ x −d z σ z , which is suitable for a particle-hole symmetric QBCP with d x = p x p y /m 0 , d z = (p 2 x − p 2 y )t/(2m 0 ), andt = 2 and is independent of spin polarization σ. We will fix 1/(2m 0 ) as our unit in the analysis. The action of interactions can be written as, dξuψ a↑ (4)ψ a↓ (3)ψ a↓ (2)ψ a↑ (1), (2π) 2 and the shorthand notation b = (ω b , p b ). Note that "4=1+2-3" is understood by energy and momentum conservation and we ignore the momentum dependence of u and g σσ ′ , which turns out to be irrelevant in RG sense. In addition, it is worth mentioning here that the present RG analysis is much simpler than that of usual Fermi liquids in the sense that there is no Fermi surface (p F = 0). This fact thus waives the complexity brought by the Fermi surface, a situation similar to φ 4 -theory.
a, σ a, σ Let us now sketch how to obtain one-loop RG equations, following the standard perturbative RG procedure by Shankar. 50 Defining a momentum cutoff Λ and ψ < (ψ > ) as slow (fast) modes with 0 < p < Λ/s (Λ/s < p < Λ), the key formula we use to derive the renormalized action by integrating out the fast modes is given by Note that in deriving last equality we have used the cumulant expansion. Z 0 denotes the non-interacting partition function and 0,> represents average with respect to the fast modes of action S 0 . The collected exponent of the exponential terms now results in the renormalized action S ′ .
Before examining the renormalization of coupling parameters u and g σσ ′ , we first notice that, to one-loop level, the quartic terms in S int do induce quadratic terms from the tadpole diagram shown in Fig. 10, which are momentum independent. This would indicate there is no non-interacting fixed point in our system. However, explicit calculations show that all contributions from such type of diagram are zero and hence our starting FP survives without flowing (t, ψ aσ are not renormalized), in contrast to the case of the Luttinger liquid.
Next, we turn our attention to the diagrams that renormalize coupling parameters. As shown in Figs. 11(a)-(c), they are, respectively, total contributions to the renormalized u, g σσ , and g σσ . Obviously, they have standard structures such as 'BCS', 'ZS' (zero sound), and 'ZS ′ ' used in Shankar's seminal work. 50 Note that we have taken the convention that1 = 2,↑ =↓ and vice versa. The essential fermion loop integrals are listed below (without vertex): andt, α(t), γ(t) are defined in the main text [just below Eq. (13)]. The superscript, σ ′ , of G 0 can be either σ orσ. The presence of dl is due to the approximation we have used, ln s = ln(1 + dl) ≈ dl in the limit s → 1. In order to evaluate these integrals, we have set all external ω = p = 0 and made q x = q cos θ, q y = q sin θ. Combining all the Feynman diagrams in Fig. 11 and keeping in mind the order of fermion operators, the straightforward algebra gives rise a set of RG equations, as we desire, Eq. (13).
Finally, we investigate the RG flow beginning with weak coupling regime for two typical cases: (a) u > g σσ = g σσ ; (b) u < g σσ = g σσ . As one can see in Fig. 12, in either case, they all flow to strong coupling regime and hence the non-interacting fixed point is actually unstable against short-range repulsions.  The existence of the QSH state can be justified by a non-trivial Z 2 -valued invariant ν = 1. This Z 2 index represents the topological nature of the QSH state and leads to topologically protected gapless edge states when the 2D system is made open with boundaries. Here, we present a method, first invented by Fu and Kane 9 , to derive Z 2 index explicitly. This method essentially takes the advantage of inversion symmetry such that the index can be related to the parity eigenvalues ξ 2m (Γ i ) of 2m occupied energy bands at four TRIM: Γ i = Γ, X, Y, and M.
For the purpose of instruction, we focus on 1/3 filling with ε A,B,C = 0 and hence only the lower band is relevant. As shown in Sec III, the presence of V 2 can dynamically generate the "spin-orbit" interaction, which is equivalent to add the following term into Eq. (2).
where w k ≡ 2 j=1 (−1) j+1 cos(k·a ′ j ) with a ′ 1 = a 1 + a 2 , and a ′ 2 = a 2 −a 1 . The +(−) sign refers to spin up (down) fermions. The eigenvalues of H(k) = H 0 (k) + H SO (k) are ǫ ± (k) = ± b k + 4λ 2 SO w 2 k and ǫ 0 (k) = 0 and their corresponding eigenstates can be written in the form, where the expressions of the components q lk (l = A, B, C) and the normalization factor G nk for each band n(= 0, ±) are given in Table I.
At TRIM, we have eigenstates |u −Γ = ( 1  Table I, respectively. Now, let us make a gauge transformation on eigenstates such that only one coordinate R would be assigned to each unit cell (instead of R + a i for each atom within the cell), and the Bloch Hamiltonian H(k) now obeys H(k + G) = H(k) and, with the inversion operator 2 ) T . Picking up any one of A-sites as our inversion center, we can determine the parity eigenvalue for each TRIM by examining the wave functions (eigenstates) in real space. Since moving a unit cell by R is simply multiplying our eigenstate at Γ i by a factor e iΓi·R , we now see the parity eigenvalues for Γ i =(Γ, X, Y, M), are, respectively, P =(+, +, +, −). The Z 2 index ν is then determined by where δ i = N m=1 ξ 2m (Γ i ), and N = 1 indicating either spin-up or spin-down band in our example here. Therefore, ν = 1 suggests that the ground state of our system at 1/3 filling has topological nature. We summarize the parity eigenvalue of each band in Table II. Finally, for ε A = 0, one can follow similar procedure and the conclusion of ν = 1 is unchanged. Γi n = + n = − n = 0