Interaction-driven topological and nematic phases on the Lieb lattice

Wei-Feng Tsai; Chen Fang; Hong Yao; Jiangping Hu

doi:10.1088/1367-2630/17/5/055016

1. Introduction

The search for new topological states of matter has not stopped since the discovery of the quantum Hall state in the 1980s [1]. In particular, in recent years, the study of 'topological insulators' (TI) has become 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₂ topological invariant and robust gapless edge states in two dimensions (2D) [4–7] or surface states in three dimensions (3D) [8, 9], against moderate perturbations that preserve time reversal symmetry (TRS).

Generally, such topological insulating states 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 that exhibit the TI phase are found in, for instance, the HgTe/CdTe quantum wells (2D), Bi_x ${\rm S}{{{\rm b}}_{1-x}}$ , Bi₂Se₃, Bi₂Te₃ (3D), and so on [11–15]. However, an alternative route to TI is possible and it is associated with the concept of the topological Mott insulator, first studied in [16] 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 geometries, e.g., honeycomb, [16–22] checkerboard, [23] kagome, [24, 25] diamond lattices [27], and in the low-energy continuum theory [26].

Here we show that topological states can be generally realized in a system with quadratic band crossing points (BCPs), which are symmetry protected at the non-interacting level. Near such kinds of BCP, instability toward phases with broken 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 figure 1. With only nearest-neighbor (NN) hoppings, there is a dispersionless (flat) band in the middle of the band structure. The three-band touching point is the result of a combined effect of (a) crystalline symmetry group of the Lieb lattice ( ${{D}_{4h}}$ -group), (b) spin rotation symmetry (or the absence of spin-orbital coupling), and (c) equal intra-sublattice hopping amplitudes and vanishing hopping amplitudes between B and C sublattices. Depending on the values of the onsite potentials or the presence of further range hoppings, the band crossing feature between this flat band and the two other bands may include: (i) The flat band touches upper and lower linearly dispersing bands at one point when (a, b, c) are all satisfied. (ii) The flat band can be isolated. For instance, one can add intrinsic spin–orbit couplings, i.e., breaking (b), as discussed by Weeks and Franz [29]. (iii) When the onsite potential on the A sublattice is not equal to that of B/C sublattices [i.e., breaking (c)], the flat band only touches one of the other two bands, which becomes quadratically dispersing, instead of linearly dispersing [30]. Thus, these choices could variegate our results. Second, a nearly flat band has effectively large correlation effects due to the small bandwidth, leading to fractional Chern insulating phases if the band has a non-zero Chern number (see section 4). Finally, the 2D Lieb lattice has been the most important building block in many 3D perovskite materials featured with complex phase diagram and strong 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).

**Figure 1.** Schematic plot of the 2D (extended) Lieb lattice. The dotted curve indicates the unit cell, which contains three sublattice sites (A, B, C).
Download figure:
Standard image High-resolution image

In this paper, we start with the construction of the explicit model and reveal the topological nature of the BCP at the non-interacting level. We then examine the consequence of such topological BCP, i.e., with a symmetry protected quadratic dispersion, under the presence of short-range repulsive interactions. We investigate various symmetry breaking instabilities at BCPs within a self-consistent mean-field(MF) 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 open a gap and gain energy: one is to open a full gap at BCP, and the other one is to split the BCP into two Dirac points (each with Berry flux π), but at the price of broken C₄ 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₄ symmetry broken down to C₂). 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 nematic-spin-nematic semi-metallic phase with Dirac nodes [31–33]. Thus, we clearly demonstrate that, in principle, correlated systems with Lieb lattice structure can be host to various non-trivial phases including TIs.

This paper is organized as follows. In section 2, we define the model and demonstrate the topological nature of its BCP from both momentum and real-space points of view. Next, the consequences of introducing short-range repulsions are discussed for spinless fermions in section 3.1 and for spinful fermions in section 3.2, respectively. Finally, we discuss some issues and make conclusions in section 4.

2. The lattice model

We begin with the tight-binding model for non-interacting fermions,

$\begin{eqnarray}&&{{\mathcal{H}}_{0}}=-\mathop{\sum }\limits_{ij}{{t}_{ij}}c_{i}^{\dagger }{{c}_{j}}+\mathop{\sum }\limits_{i}{{\varepsilon }_{i}}c_{i}^{\dagger }{{c}_{i}},\end{eqnarray} \tag{ 1 }$

where $c_{i}^{\dagger }$ 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 figure 1. For simplicity, we take the hopping amplitudes, ${{t}_{ij}}=t$ , between NN sites $\langle ij\rangle$ , 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₄ point group symmetry dictates that the onsite potentials on B and C are equal, ${{\varepsilon }_{B}}={{\varepsilon }_{C}}$ . Although ${{\varepsilon }_{(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\equiv 1$ , the lattice constant $a\equiv 1$ , and, without loss of generality, ${{\varepsilon }_{B}}={{\varepsilon }_{C}}=0$ .

2.1. Band structure

The band structure of equation (1) can be obtained by transforming ${{\mathcal{H}}_{0}}$ into momentum space,

$\begin{eqnarray}&&{{\mathcal{H}}_{0}}=\mathop{\sum }\limits_{{\bf k}}\psi _{{\bf k}}^{\dagger }{{H}_{0}}({\bf k}){{\psi }_{{\bf k}}},\end{eqnarray} \tag{ 2 }$

where the fermion spinor, $\psi _{{\bf k}}^{\dagger }=(c_{A{\bf k}}^{\dagger },c_{B{\bf k}}^{\dagger },c_{C{\bf k}}^{\dagger })$ , with sublattice (basis) index $A,B,C$ and ${\bf k}=({{k}_{x}},{{k}_{y}})$ . Defining the displacement vectors, ${{{\bf a}}_{1}}=(1/2,0)$ and ${{{\bf a}}_{2}}=(0,1/2)$ , ${{H}_{0}}({\bf k})$ is of the form

$\begin{eqnarray}{{H}_{0}}({\bf k})=\left( \begin{array}{ccccccccccccccc} {{\varepsilon }_{A}} & -2t{\rm cos} \left( {\bf k}\cdot {{{\bf a}}_{1}} \right) & -2t{\rm cos} \left( {\bf k}\cdot {{{\bf a}}_{2}} \right) \\ {} & 0 & 0 \\ {} & {} & 0 \\ \end{array} \right),\end{eqnarray} \tag{ 3 }$

where the lower triangular matrix is understood to be filled for keeping the whole matrix hermitian. In this notation, the first Brillouin zone (FBZ) is a square with four time reversal invariant momenta (TRIM): ${\boldsymbol{ \Gamma }} =(0,0),{\bf X}=(\pi ,0),{\bf M}=(\pi ,\pi )$ , and ${\bf Y}=(0,\pi )$ [see the inset of figure 2(a)]. The energy spectrum consists of two dispersive bands, ${{\epsilon }_{\pm }}({\bf k})=\frac{1}{2}({{\varepsilon }_{A}}\pm \sqrt{\varepsilon _{A}^{2}+4{{b}_{{\bf k}}}})$ with ${{b}_{{\bf k}}}=\sum _{i=1}^{2}{{[2t{\rm cos} ({\bf k}\cdot {{{\bf a}}_{i}})]}^{2}}$ , and one dispersionless flat band, ${{\epsilon }_{0}}({\bf k})=0$ .

**Figure 2.** Band structure of the model along the high symmetry lines in the FBZ for (a) ${{\varepsilon }_{A}}=0$ and for (b) ${{\varepsilon }_{A}}=-4t$ . In case (b), when ${{t}^{{\prime\prime} }}=0.1t$ (with the constraint mentioned in the text) is present, the BCP at ${\bf M}$ becomes standard QBCP and the spectrum is changed slightly as shown by dashed lines. Inset of (a): FBZ of the Lieb lattice. There are four time reversal invariant momenta: ${\boldsymbol{ \Gamma }} ,{\bf X},{\bf M}$ , and ${\bf Y}$ .
Download figure:
Standard image High-resolution image

**Figure 2.** Band structure of the model along the high symmetry lines in the FBZ for (a) ${{\varepsilon }_{A}}=0$ and for (b) ${{\varepsilon }_{A}}=-4t$ . In case (b), when ${{t}^{{\prime\prime} }}=0.1t$ (with the constraint mentioned in the text) is present, the BCP at ${\bf M}$ becomes standard QBCP and the spectrum is changed slightly as shown by dashed lines. Inset of (a): FBZ of the Lieb lattice. There are four time reversal invariant momenta: ${\boldsymbol{ \Gamma }} ,{\bf X},{\bf M}$ , and ${\bf Y}$ .
Download figure:
Standard image High-resolution image

An important feature of this model is that the presence or absence of ${{\varepsilon }_{A}}$ can change electronic properties dramatically. When ${{\varepsilon }_{A}}=0$ , the flat band touches two linearly dispersing bands at ${\bf M}$ point in the FBZ [type (i)], where the linear bands meet as if there was a 'Dirac point'. However, the touching point in fact has completely different structure. It becomes clear once we expand ${{H}_{0}}({\bf k})$ around ${\bf M}$ point with ${\bf k}={\bf M}+{\bf p},|{\bf p}|\ll \;$ 1. To the first order in ${\bf p}$ , ${{H}_{0}}({\bf k})$ can be written as

$\begin{eqnarray}&&{{H}_{0}}({\bf B})\sim {{v}_{F}}{\bf L}\cdot {\bf B},\end{eqnarray} \tag{ 4 }$

where the Fermi velocity ${{v}_{F}}=t$ , ${\bf B}=({{p}_{x}},{{p}_{y}},0)$ , and the (pseudo) spin-1 matrices are defined as

$\begin{eqnarray}{{L}_{x}}=\left( \begin{array}{ccccccccccccccc} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0 \\ \end{array} \right),{{L}_{y}}=\left( \begin{array}{ccccccccccccccc} 0 & 0 & 1 \\ 0 & 0 & 0 \\ 1 & 0 & 0 \\ \end{array} \right),{{L}_{z}}=\left( \begin{array}{ccccccccccccccc} 0 & 0 & 0 \\ 0 & 0 & -i \\ 0 & i & 0 \\ \end{array} \right),\end{eqnarray} \tag{ 5 }$

obeying Lie algebra of SU(2), i.e., $[{{L}_{i}},{{L}_{j}}]=i{{\epsilon }_{ijk}}{{L}_{k}}$ [39, 40], 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 [41] on the Lieb lattice. Viewing the low-energy effective Hamiltonian ${{H}_{0}}({\bf B})$ as a 'spin' ${\bf L}$ in an external 'magnetic field' ${\bf B}$ , its eigenvalues can be easily read out as $\epsilon ({\bf B})={{v}_{F}}|{\bf p}|{{l}_{{\bf p}}}$ , where ${{l}_{{\bf p}}}=0,\pm 1$ are the quantized angular momenta along the axis parallel to ${\bf B}$ in three dimensions.

When ${{\varepsilon }_{A}}\ne 0$ , however, the spin-1 structure mentioned previously is no longer valid. The flat band now touches only one dispersive (massive) energy band either above or below at ${\bf M}$ point, depending on the sign of ${{\varepsilon }_{A}}$ [see figure 2(b)]. To make this structure transparent, we again expand ${{H}_{0}}({\bf k})$ around ${\bf M}$ with small ${\bf p}$ . Assuming $|{\bf p}|\ll |{{\varepsilon }_{A}}/t|$ and ${{\varepsilon }_{A}}\lt 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,

$\begin{eqnarray}H_{0}^{eff}({\bf p})\sim \frac{1}{{{m}_{0}}}\left( \begin{array}{ccccccccccccccc} p_{x}^{2} & {{p}_{x}}{{p}_{y}} \\ {{p}_{x}}{{p}_{y}} & p_{y}^{2} \\ \end{array} \right)={{d}_{I}}I+{{d}_{x}}{{\sigma }_{x}}+{{d}_{z}}{{\sigma }_{z}},\end{eqnarray} \tag{ 6 }$

where ${{m}_{0}}=-{{\varepsilon }_{A}}/v_{F}^{2}$ . In the last equality, we express H₀^eff in terms of the identity and Pauli matrices with ${{d}_{I}}=\frac{1}{2}(p_{x}^{2}+p_{y}^{2})m_{0}^{-1}$ , ${{d}_{x}}={{p}_{x}}{{p}_{y}}/{{m}_{0}}$ , and ${{d}_{z}}=\frac{1}{2}(p_{x}^{2}-p_{y}^{2})m_{0}^{-1}$ . Interestingly, if we further allow small, lattice symmetry unbroken third-neighbor hoppings ${{t}^{{\prime\prime} }}\gt 0$ (but forbid to hop whenever here is a site in the middle of the path), the flat band becomes slightly dispersive and the effective Hamiltonian changes to ${{d}_{I}}=\frac{1}{2}(p_{x}^{2}+p_{y}^{2})(m_{0}^{-1}-{{t}^{{\prime\prime} }})$ , ${{d}_{x}}={{p}_{x}}{{p}_{y}}/{{m}_{0}}$ , and ${{d}_{z}}=\frac{1}{2}(p_{x}^{2}-p_{y}^{2})(m_{0}^{-1}+{{t}^{{\prime\prime} }})$ without removing BCP [see figure 2(b)]. Such point at ${\bf p}=0$ is the so-called quadratic BCP (QBCP), which has been studied recently by several research groups. [23, 34–36] One of the key features for a 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 ${{\varepsilon }_{A}}\ne 0$ case and show that the BCP in our model is not only topologically non-trivial, but also makes the system be a potential host to a topological phase under weak repulsive interactions.

2.2. Topological nature of the band touching

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 equation (6) as $H_{0}^{eff}({\bf p})={{d}_{I}}I+{{{\bf B}}^{\prime }}({\bf p})\cdot {\boldsymbol{ \sigma }}$ , where the 'magnetic field' ${{{\bf B}}^{\prime }}({\bf p})=({{d}_{x}},0,{{d}_{z}})$ . This effective Hamiltonian now represents a spin-1/2 particle sitting in a magnetic field ${{{\bf B}}^{\prime }}$ , which has a vortex structure at ${\bf p}=0$ , as shown in figure 3(a). For comparison, recall that for ${{\varepsilon }_{A}}=0$ we instead have a spin-1 particle in an external field ${\bf B}$ [equation (4)], whose structure is shown in figure 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 counterintuitively, both cases are associated with the same Berry phase of the BCP, [37] which can be calculated precisely by

$\begin{eqnarray}&&{{B}^{n}}=i{{\oint }_{\Gamma }}{\rm d}{\bf p}\cdot \langle {{u}_{n{\bf p}}}|{{\nabla }_{{\bf p}}}|{{u}_{n{\bf p}}}\rangle ,\end{eqnarray} \tag{ 7 }$

where Γ is a contour in the ${\bf p}$ space enclosing the touching point, n denotes any one of the involved bands, and $|{{u}_{n{\bf p}}}\rangle$ represents the Bloch wave function for nth band. A simple argument solves this puzzle. The line integral along any loop enclosing 0 in ${\bf p}$ space given previously is known to be 1/2 (1) times the solid angle subtended by ${{{\bf B}}^{\prime }}({\bf p})$ [ ${\bf B}({\bf p})$ ] from the origin for a spin-1/2 (1) particle. Thus, ${{B}^{n}}=2\pi W/2=2\pi$ in the former case, which is just equivalent to ${{B}^{n}}=2\pi W=2\pi$ in the latter one. In fact, when Bⁿ = 0, any infinitesimal mixing (perturbation) between bands would lift the degeneracy. With non-vanishing ${{B}^{n}}=\pm 2\pi$ , we confirm that the BCP on the Lieb lattice at the non-interacting level is topologically stable (i.e., not opening a gap) as long as the spinless system preserves both TRS and C₄ point group symmetry. Note that C₄ symmetry in our model is quite essential, as a similar QBCP happens in the A-B stacking bilayer graphene (with C₃ symmetry) while it can easily decay into Dirac BCPs and thus is topologically unstable [23, 28].

**Figure 3.** Vortex structure of the BCP in ${\bf p}$ space: (a) the planar vector field ${{{\bf B}}^{\prime }}$ ; (b) the planar vector field ${\bf B}$ .
Download figure:
Standard image High-resolution image

**Figure 3.** Vortex structure of the BCP in ${\bf p}$ space: (a) the planar vector field ${{{\bf B}}^{\prime }}$ ; (b) the planar vector field ${\bf B}$ .
Download figure:
Standard image High-resolution image

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 that are extended along non-contractible loops winding around the whole lattice with periodic boundary conditions (i.e., a torus) [38]. 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 ${{\mathcal{H}}_{0}}$ (Wannier states). Taking ${\bf R}$ to be the coordinate of the central site of the shaded plaquette shown in figure 4, we find that the creation operator for the localized eigenstate at ${\bf R}$ can be written as

$\begin{eqnarray}&&\mathcal{A}_{{\bf R}}^{\dagger }=\frac{1}{2\sqrt{6}}[\mathop{\sum }\limits_{j=1}^{4}{{(-1)}^{j}}(2c_{{\bf R}+{{{\bf b}}_{j}}}^{\dagger }-c_{{\bf R}+{{{\bf b}}_{2j+3}}}^{\dagger }-c_{{\bf R}+{{{\bf b}}_{2j+4}}}^{\dagger })],\end{eqnarray} \tag{ 8 }$

where ${{{\bf b}}_{1}}={{{\bf a}}_{1}},{{{\bf b}}_{2}}={{{\bf a}}_{2}},{{{\bf b}}_{3}}=-{{{\bf b}}_{1}},{{{\bf b}}_{4}}=-{{{\bf b}}_{2}}$ , ${{{\bf b}}_{5}}=(1,-1/2),{{{\bf b}}_{6}}=(1,1/2),{{{\bf b}}_{7}}=(1/2,1),{{{\bf b}}_{8}}=(-1/2,1)$ , ${{{\bf b}}_{9}}=-{{{\bf b}}_{5}},{{{\bf b}}_{10}}=-{{{\bf b}}_{6}},{{{\bf b}}_{11}}=-{{{\bf b}}_{7}},{{{\bf b}}_{12}}=-{{{\bf b}}_{8}}$ . The key reason for these states being localized is rooted in the fact that all A-sites have vanishing amplitudes and remain zero after the action of ${{\mathcal{H}}_{0}}$ on them due to destructive interference.

**Figure 4.** Schematic plot of the localized eigenstate at ${\bf R}$ on the Lieb lattice ( ${{\varepsilon }_{A}}\ne 0$ ). Only those sites in the shaded area have non-zero weights, whose magnitudes are denoted by numbers (up to a normalization factor). The weights on A sublattice (blue dots) are all zero.
Download figure:
Standard image High-resolution image

**Figure 4.** Schematic plot of the localized eigenstate at ${\bf R}$ on the Lieb lattice ( ${{\varepsilon }_{A}}\ne 0$ ). Only those sites in the shaded area have non-zero weights, whose magnitudes are denoted by numbers (up to a normalization factor). The weights on A sublattice (blue dots) are all zero.
Download figure:
Standard image High-resolution image

In localized-state language, the existence of the BCP in our model with ${{\varepsilon }_{A}}\ne 0$ is equivalent to state that the dimension of the space expanded by independent localized eigenstates with zero energy has a dimension that 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 equation (8) seem to form N linearly independent states with zero energy. For our model with periodic boundary conditions, however, the following relation,

$\begin{eqnarray}&&\mathcal{A}_{{\bf q}=(\pi ,\pi )}^{\dagger }=\mathop{\sum }\limits_{{\bf R}}{{e}^{i{\bf q}\cdot {\bf R}}}\mathcal{A}_{{\bf R}}^{\dagger }=0,\end{eqnarray} \tag{ 9 }$

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 figures 5(a) and (b). When ${{\mathcal{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.

**Figure 5.** The two non-contractible loop states winding the lattice along periodic x direction (a) and y direction (b). Only those sites in the shaded area have non-zero weights (denoted by numbers in the plots). The weights on A sublattice (blue dots) are all zero.
Download figure:
Standard image High-resolution image

3. Interaction-driven instabilities

The existence of such symmetry/topology-protected BCPs 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 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 $|{{\varepsilon }_{A}}|\gg t$ limit,

$\begin{eqnarray}&&{{\mathcal{H}}^{eff}}=\mathcal{H}_{0}^{eff}+\mathcal{H}_{\operatorname{int}}^{eff},\end{eqnarray} \tag{ 10 }$

where the free part ${{\mathcal{H}}_{0}}$ reads as

$\begin{eqnarray}&&\mathcal{H}_{0}^{eff}=\mathop{\sum }\limits_{\sigma }\int d{\bf r}\Psi _{\sigma }^{\dagger }({\bf r})H_{0}^{eff}{{\Psi }_{\sigma }}({\bf r}).\end{eqnarray} \tag{ 11 }$

The subscript σ denotes spin polarization and the fermion field $\Psi _{\sigma }^{\dagger }=(\psi _{1\sigma }^{\dagger },\psi _{2\sigma }^{\dagger })$ with 1, 2 represents orbital (i.e., two touching bands) degrees of freedom. Note that in the momentum space the expression for H^eff₀ is given by equation (6), and is independent of σ. For simplicity, we take small ${{t}^{{\prime\prime} }}=-{{t}^{2}}/{{\varepsilon }_{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_{0}^{-1}\gt {{t}^{{\prime\prime} }}$ . The projected interacting part includes only (1) intra-orbital and (2) inter-orbital contact interactions,

$\begin{eqnarray}\begin{array}{ccccccccccccccc} \mathcal{H}_{\operatorname{int}}^{eff} & = & \mathcal{H}_{1}^{eff}+\mathcal{H}_{2}^{eff}, \\ \mathcal{H}_{1}^{eff} & = & \mathop{\sum }\limits_{a=1}^{2}u\int {\rm d}{\bf r}\psi _{a\uparrow }^{\dagger }({\bf r})\psi _{a\downarrow }^{\dagger }({\bf r}){{\psi }_{a\downarrow }}({\bf r}){{\psi }_{a\uparrow }}({\bf r}), \\ \mathcal{H}_{2}^{eff} & = & \mathop{\sum }\limits_{\sigma ,{{\sigma }^{\prime }}}{{g}_{\sigma {{\sigma }^{\prime }}}}\int {\rm d}{\bf r}\psi _{1\sigma }^{\dagger }({\bf r})\psi _{2{{\sigma }^{\prime }}}^{\dagger }({\bf r}){{\psi }_{2{{\sigma }^{\prime }}}}({\bf r}){{\psi }_{1\sigma }}({\bf r}), \\ \end{array}\end{eqnarray} \tag{ 12 }$

where u and ${{g}_{\sigma {{\sigma }^{\prime }}}}$ are intra-orbital and inter-orbital coupling parameters, respectively.

For the chemical potential $\mu =0$ , the non-interacting system $\mathcal{H}_{0}^{eff}$ leads to one Fermi point at ${\bf p}=0$ with non-vanishing 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 $[{{\psi }_{a\sigma }}({\bf r})]=1$ ; in the interacting part the coupling constants $[u]=[{{g}_{\sigma {{\sigma }^{\prime }}}}]=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 one-loop order, the coupled RG equations for the coupling parameters are

$\begin{eqnarray}\begin{array}{ccccccccccccccc} \frac{du}{dl} & = & (2\gamma -\alpha ){{u}^{2}}+2\gamma {{g}_{\uparrow \uparrow }}u-2\gamma {{g}_{\uparrow \uparrow }}{{g}_{\uparrow \downarrow }}, \\ \frac{d{{g}_{\uparrow \uparrow }}}{dl} & = & \alpha g_{\uparrow \uparrow }^{2}-2\gamma {{g}_{\uparrow \downarrow }}u+\gamma g_{\uparrow \downarrow }^{2}, \\ \frac{d{{g}_{\uparrow \downarrow }}}{dl} & = & (\alpha -2\gamma )g_{\uparrow \downarrow }^{2}-2\gamma {{g}_{\uparrow \uparrow }}u+2\gamma {{g}_{\uparrow \uparrow }}{{g}_{\uparrow \downarrow }}, \\ \end{array}\end{eqnarray} \tag{ 13 }$

where we have used the fact that ${{g}_{\uparrow \uparrow }}={{g}_{\downarrow \downarrow }}$ and ${{g}_{\uparrow \downarrow }}={{g}_{\downarrow \uparrow }}$ . l denotes the momentum rescaling $p\to p{{e}^{-l}}$ and the coefficients $\alpha =\frac{1}{2{{\pi }^{2}}|\bar{t}|}K(\sqrt{1-{{{\bar{t}}}^{-2}}})$ and $\gamma =\frac{{{{\bar{t}}}^{2}}E(1-{{{\bar{t}}}^{-2}})-K(1-{{{\bar{t}}}^{-2}})+|\bar{t}|(E(1-{{{\bar{t}}}^{2}})-K(1-{{{\bar{t}}}^{2}}))}{8{{\pi }^{2}}(-1+{{{\bar{t}}}^{2}})|\bar{t}|}$ , with K(x) [E(x)], the elliptic function of first (second) kind and $\bar{t}=(m_{0}^{-1}+{{t}^{{\prime\prime} }})/m_{0}^{-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}_{\sigma {{\sigma }^{\prime }}}}$ vanish. Moreover, we find that given generic bare coupling parameters ( $u,{{g}_{\sigma {{\sigma }^{\prime }}}}\gt 0$ ), at least one of them diverges first when reaching a low enough energy scale. This indicates that the non-interacting FP is an unstable FP, which can drive the system to the strong coupling regime in the presence of short-range repulsions. In addition, it is worth mentioning that by setting ${{g}_{\uparrow \downarrow }}=u=0$ , we reduce the RG equation back to the spinless case,

$\begin{eqnarray}&&\frac{d{{g}_{\uparrow \uparrow }}}{dl}=\alpha g_{\uparrow \uparrow }^{2}\end{eqnarray} \tag{ 14 }$

with $\alpha \gt 0$ , which is consistent with the work done in [23] and, importantly, it implies that short-range repulsions are again marginally relevant.

3.1. 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 MF 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 equation (1), with short-range interacting terms,

$\begin{eqnarray}&&{{\mathcal{H}}_{\operatorname{int}}}={{V}_{1}}\mathop{\sum }\limits_{\langle ij\rangle }{{n}_{i}}{{n}_{j}}+{{V}_{2}}\mathop{\sum }\limits_{\langle \langle ij\rangle \rangle }{{n}_{i}}{{n}_{j}},\end{eqnarray} \tag{ 15 }$

where V₁ and V₂ are repulsive coupling constants for NN and NNN interactions, respectively. ${{n}_{i}}=c_{i}^{\dagger }{{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 ${{\varepsilon }_{A}}\lt 0$ ( ${{\varepsilon }_{A}}\gt 0$ ). We proceed by treating ${{\mathcal{H}}_{\operatorname{int}}}$ in the MF approximation, including both the onsite and bond MF decoupling particle-hole channels,

$\begin{eqnarray}&&{{n}_{i}}{{n}_{j}}\to {{n}_{i}}\langle {{n}_{j}}\rangle +{{n}_{j}}\langle {{n}_{i}}\rangle -\langle {{n}_{i}}\rangle \langle {{n}_{j}}\rangle ,\end{eqnarray} \tag{ 16 }$

$\begin{eqnarray}&&{{n}_{i}}{{n}_{j}}\to -{{\phi }_{ij}}c_{j}^{\dagger }{{c}_{i}}-\phi _{ij}^{*}c_{i}^{\dagger }{{c}_{j}}+{{\phi }_{ij}}\phi _{ij}^{*},\end{eqnarray} \tag{ 17 }$

where ${{\phi }_{ij}}=\phi _{ji}^{*}=\langle c_{i}^{\dagger }{{c}_{j}}\rangle$ 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₄ symmetry down to C₂. In particular, it does not break the translational symmetry by any lattice vectors, and thus is in contrast to the conventional charge density wave (CDW) order, which enlarges the unit cell due to translational symmetry breaking [see figure 6]. This phase 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, $\eta =\frac{1}{8}{{\sum }_{{{\delta }^{\prime }}}}(\langle c_{Bi}^{\dagger }{{c}_{Bi}}\rangle -\langle c_{Ci+{{\delta }^{\prime }}}^{\dagger }{{c}_{Ci+{{\delta }^{\prime }}}}\rangle )$ , where ${{\delta }^{\prime }}=\pm \hat{x}/2\pm \hat{y}/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, $\langle c_{Ai}^{\dagger }{{c}_{Ai}}\rangle$ $\;=\;\frac{2}{3}+\rho$ , $\langle c_{Bi}^{\dagger }{{c}_{Bi}}\rangle$ $\;=\;\frac{2}{3}-\frac{\rho }{2}+\eta$ , and $\langle c_{Ci}^{\dagger }{{c}_{Ci}}\rangle$ $\;=\;\frac{2}{3}-\frac{\rho }{2}-\eta$ .

**Figure 6.** The schematic plots for (a) site charge nematic state and (b) site nematic-spin-nematic state. Since the electronic information of A sublattice sites (blue dots) is irrelevant here, only those of B and C sublattice sites are compared. Each open circle represents the local charge density, whose magnitude depends on the size of it. The arrow at each site represents a net spin polarization.
Download figure:
Standard image High-resolution image

As ${{\varepsilon }_{A}}\gt 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 non-zero 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 ${{Q}_{1}}=\frac{1}{4}\operatorname{Re}[{{\sum }_{\delta =\pm \hat{x}}}\langle c_{Ai}^{\dagger }{{c}_{Bi+\delta }}\rangle -{{\sum }_{\delta =\pm \hat{y}}}\langle c_{Ai}^{\dagger }{{c}_{Ci+\delta }}\rangle ]$ , or in the form of ${{Q}_{2}}=\frac{1}{4}{{\sum }_{{{\delta }^{\prime }}}}{{D}_{{{\delta }^{\prime }}}}\operatorname{Re}\langle c_{Bi}^{\dagger }{{c}_{Ci+{{\delta }^{\prime }}}}\rangle$ for ${{D}_{{{\delta }^{\prime }}=\pm (\hat{x}/2-\hat{y}/2)}}=1$ and ${{D}_{{{\delta }^{\prime }}=\pm (\hat{x}/2+\hat{y}/2)}}=-1$ . The subscript of the order parameters indicates their origin of either V₁ or V₂ repulsion.

(ii) Current-loop state. This type of phase is featured by spontaneous TRS breaking. The most probable current patterns that preserve translational invariance with no (charge) source and drain present on the lattice sites are shown in figures 7(a)–(d). Each state basically comes from the non-vanishing imaginary part of certain bond orders in the MF decouplings and may behave Hall-insulating [7(a)], semi-metallic [7(b), 7(d)], or insulating [7(c)]. In particular, the most significant one is case (a), which exhibits QAH effect with order parameter, ${{\Phi }_{2}}=\frac{1}{4}{\rm Im}[{{\sum }_{{{\delta }^{\prime }}}}{{D}_{{{\delta }^{\prime }}}}\langle c_{Bi}^{\dagger }{{c}_{Ci+{{\delta }^{\prime }}}}\rangle ]$ . 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 [42]. We compute the Chern number for each band within this state and find that (1) for $|{{\varepsilon }_{A}}|\gt 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 ${{\varepsilon }_{A}}=0$ , the middle flat band carries zero Chern number, while the upper and lower bands carry ±1, respectively.

**Figure 7.** The current loop ground states for (a) ${{\Phi }_{2}}\ne 0$ and others are zero (QAH state); (b) ${{\Phi }_{1}},\Phi _{2}^{\prime }\ne 0$ and others are zero (Varma ${{\Theta }_{I}}$ loop)(c) $\Phi _{1}^{\prime },\Phi _{2}^{{\prime\prime} }=\Phi _{2}^{\prime\prime\prime}\ne 0$ (Varma ${{\Theta }_{II}}$ loop); (d) $\Phi _{1}^{\prime }=\Phi _{1}^{{\prime\prime} },\Phi _{2}^{\prime\prime\prime}\ne 0$ .
Download figure:
Standard image High-resolution image

**Figure 7.** The current loop ground states for (a) ${{\Phi }_{2}}\ne 0$ and others are zero (QAH state); (b) ${{\Phi }_{1}},\Phi _{2}^{\prime }\ne 0$ and others are zero (Varma ${{\Theta }_{I}}$ loop)(c) $\Phi _{1}^{\prime },\Phi _{2}^{{\prime\prime} }=\Phi _{2}^{\prime\prime\prime}\ne 0$ (Varma ${{\Theta }_{II}}$ loop); (d) $\Phi _{1}^{\prime }=\Phi _{1}^{{\prime\prime} },\Phi _{2}^{\prime\prime\prime}\ne 0$ .
Download figure:
Standard image High-resolution image

The other possible current-loop states, however, are not topological insulating. For case (b) (Varma ${{\Theta }_{I}}$ loop state [43]), it is semi-metallic with the order parameter given by, simultaneously, ${{\Phi }_{1}}=\frac{1}{4}{\rm Im}[{{\sum }_{\delta =\pm \hat{x}/2}}\langle c_{Ai}^{\dagger }{{c}_{Bi+\delta }}\rangle -{{\sum }_{\delta =\pm \hat{y}/2}}\langle c_{Ai}^{\dagger }{{c}_{Ci+\delta }}\rangle ]$ and $\Phi _{2}^{\prime }=\frac{1}{4}{{\sum }_{{{\delta }^{\prime }}}}{\rm Im}\langle c_{Bi}^{\dagger }{{c}_{Ci+{{\delta }^{\prime }}}}\rangle$ . For case (c) (Varma ${{\Theta }_{II}}$ loop state [43]), 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 ${\bf k}$ it changes sign under TRS or IS [35]. This order can be described as $\Phi _{1}^{\prime }\ne 0,\Phi _{2}^{{\prime\prime} }=\Phi _{2}^{\prime\prime\prime}\ne 0$ , where $\Phi _{1}^{\prime }=\frac{1}{4}{\rm Im}[{{\sum }_{\delta =\pm \hat{x}/2}}(2\delta \cdot \hat{x})\langle c_{Ai}^{\dagger }{{c}_{Bi+\delta }}\rangle +{{\sum }_{\delta =\pm \hat{y}/2}}(2\delta \cdot \hat{y})\langle c_{Ai}^{\dagger }{{c}_{Ci+\delta }}\rangle ]$ , $\Phi _{2}^{{\prime\prime} }=\frac{1}{4}{\rm Im}[\sum \delta ^{\prime} (2{{\delta }^{\prime }}\cdot \hat{x}){{D}_{{{\delta }^{\prime }}}}\langle c_{Bi}^{\dagger }{{c}_{Ci+{{\delta }^{\prime }}}}\rangle ]$ , and $\Phi _{2}^{\prime\prime\prime}=\frac{1}{4}{\rm Im}[{{\sum }_{{{\delta }^{\prime }}}}(2{{\delta }^{\prime }}\cdot \hat{x})\langle c_{Bi}^{\dagger }{{c}_{Ci+{{\delta }^{\prime }}}}\rangle ]$ . Finally, case (d) is a semi-metallic state with broken IS and hence has no net Hall current in it. Its order is described as $\Phi _{2}^{\prime\prime\prime}\ne 0,\Phi _{1}^{\prime }=\Phi _{1}^{{\prime\prime} }\ne 0$ , where $\Phi _{1}^{{\prime\prime} }=\frac{1}{4}{\rm Im}[{{\sum }_{\delta =\pm \hat{x}/2}}(2\delta \cdot \hat{x})\langle c_{Ai}^{\dagger }{{c}_{Bi+\delta }}\rangle -{{\sum }_{\delta =\pm \hat{y}/2}}(2\delta \cdot \hat{y})\langle c_{Ai}^{\dagger }{{c}_{Ci+\delta }}\rangle ]$ .

In momentum space, the MF Hamiltonian at 2/3 filling can be now written in the matrix form,

$\begin{eqnarray}&&{{\mathcal{H}}_{MF}}=\mathop{\sum }\limits_{{\bf k}}\Psi _{{\bf k}}^{\dagger }({{H}_{{\bf k}}}-\mu I){{\Psi }_{{\bf k}}}+{{E}_{0}},\end{eqnarray} \tag{ 18 }$

with the fermion spinor, $\Psi _{{\bf k}}^{\dagger }=(c_{A{\bf k}}^{\dagger },c_{B{\bf k}}^{\dagger },c_{C{\bf k}}^{\dagger })$ , and

$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{E}_{0}} & = & -\frac{N{{{\bar{\varepsilon }}}_{A1}}({{{\bar{\varepsilon }}}_{B1}}+{{{\bar{\varepsilon }}}_{C1}})}{4{{V}_{1}}}-\frac{N{{{\bar{\varepsilon }}}_{B2}}{{{\bar{\varepsilon }}}_{C2}}}{4{{V}_{2}}} \\ {} & {} & +4N{{V}_{2}}(\delta {{t}^{{\prime} 2}}+Q_{2}^{2}+\Phi _{2}^{2}+\Phi _{2}^{{\prime} 2}+\Phi _{2}^{{\prime\prime} 2}+\Phi _{2}^{^{\prime\prime\prime}2}) \\ {} & {} & +4N{{V}_{1}}(\delta {{t}^{2}}+Q_{1}^{2}+\Phi _{1}^{2}+\Phi _{1}^{{\prime} 2}+\Phi _{1}^{{\prime\prime} 2}), \\ \end{array}\end{eqnarray} \tag{ 19 }$

where $\delta t$ ( $\delta {{t}^{\prime }}$ ) represents a renormalization due to NN(NNN) repulsions. The Hamiltonian matrix ${{H}_{{\bf k}}}$ reads

$\begin{eqnarray}{{H}_{{\bf k}}}=\left( \begin{array}{ccccccccccccccc} {{{\bar{\varepsilon }}}_{A1}}+{{\varepsilon }_{A}} & {{\Gamma }_{x}} & {{\Gamma }_{y}} \\ {} & {{{\bar{\varepsilon }}}_{B1}}+{{{\bar{\varepsilon }}}_{B2}}+{{\nu }_{y}} & {{\Gamma }_{xy}} \\ {} & {} & {{{\bar{\varepsilon }}}_{C1}}+{{{\bar{\varepsilon }}}_{C2}}+{{\nu }_{x}}, \\ \end{array} \right),\end{eqnarray} \tag{ 20 }$

where ${{\bar{\varepsilon }}_{A1}}={{V}_{1}}(8-6\rho )/3$ , ${{\bar{\varepsilon }}_{B1}}=2{{V}_{1}}(2+3\rho )/3$ , ${{\bar{\varepsilon }}_{B2}}=2{{V}_{2}}(4-3\rho -6\eta )/3$ , ${{\bar{\varepsilon }}_{C1}}=2{{V}_{1}}(2+3\rho )/3$ , ${{\bar{\varepsilon }}_{C2}}=2{{V}_{2}}(4-3\rho +6\eta )/3$ , and ${{\nu }_{x,y}}=-2{{t}^{{\prime\prime} }}{\rm cos} {{k}_{x,y}}$ ; Γ parameters are given by

$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{\Gamma }_{x,y}} & = & (-2t\mp 2{{V}_{1}}({{Q}_{1}}-i{{\Phi }_{1}})){\rm cos} \frac{{{k}_{x,y}}}{2} \\ {} & {} & -2{{V}_{1}}(\Phi _{1}^{\prime }\pm \Phi _{1}^{{\prime\prime} }){\rm sin} \frac{{{k}_{x,y}}}{2}, \\ {{\Gamma }_{xy}} & = & 4i{{V}_{2}}\Phi _{2}^{\prime }{\rm cos} \frac{{{k}_{x}}}{2}{\rm cos} \frac{{{k}_{y}}}{2}+4{{V}_{2}}\Phi _{2}^{\prime\prime\prime}{\rm sin} \frac{{{k}_{x}}}{2}{\rm cos} \frac{{{k}_{y}}}{2} \\ {} & {} & -4{{V}_{2}}({{Q}_{2}}-i{{\Phi }_{2}}){\rm sin} \frac{{{k}_{x}}}{2}{\rm sin} \frac{{{k}_{y}}}{2} \\ {} & {} & -4{{V}_{2}}\Phi _{2}^{{\prime\prime} }{\rm cos} \frac{{{k}_{x}}}{2}{\rm sin} \frac{{{k}_{y}}}{2}. \\ \end{array}\end{eqnarray} \tag{ 21 }$

Thus, the MF free energy can be expressed as

$\begin{eqnarray}&&F=-\frac{1}{\beta N}\mathop{\sum }\limits_{{\bf k}}{\rm ln} (1+{{e}^{-\beta ({{E}_{{\bf k}}}-\mu )}})+{{E}_{0}},\end{eqnarray} \tag{ 22 }$

where $\beta =1/{{k}_{B}}T$ and ${{E}_{{\bf k}}}$ are eigenvalues of ${{H}_{{\bf 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}_{{\bf k}}}$ and E₀ need to be changed accordingly due to the shift of the average charge density.

We numerically solve the coupled gap equations self-consistently and obtain the zero-temperature V₁-V₂ phase diagrams at both 1/3 ( ${{\varepsilon }_{A}}=4t$ ) and 2/3 ( ${{\varepsilon }_{A}}=-4t$ ) fillings, as seen in figures 8(a) and (b). Note that in this study, only weak short-range repulsions, i.e., ${{V}_{1}},{{V}_{2}}\leqslant t$ are considered. At 1/3 filling, we find that there are three phases in the absence of V₁: QAH phase, coexisting QAH+nematic phase, and nematic phase. Beginning with ${{V}_{2}}\ll t$ , the infinitesimal instability of QBCP leads to QAH phase by the second-order phase transition, with a T = 0 gap, ${{\Delta }_{QAH}}=({{V}_{2}}{{\Phi }_{2}})\propto \Lambda {\rm exp} (\frac{-1}{{{N}_{0}}{{V}_{2}}})$ , where N₀ denotes the finite DOS at QBCP and Λ is an energy cutoff; on the other hand, for $t\gtrsim {{V}_{2}}\gt {{V}_{2c}}\sim 0.22t$ the ground state breaks C₄ symmetry spontaneously down to C₂ and exhibits insulating nematic phase with a gap ${{\Delta }_{N}}\sim \eta$ . In this phase, we find that the site-nematic order (η) is the dominant one and a small component of the bond-nematic order (Q₁) accompanies with it. In fact, we notice that the bond-nematic order cannot be induced by V₁ itself.

**Figure 8.** Schematic zero-temperature MF phase diagrams for spinless fermions on the (a) 1/3 filled lattice with ${{\varepsilon }_{A}}/t=4$ , ${{t}^{{\prime\prime} }}/t=-0.1$ and on the (b) 2/3 filled lattice with ${{\varepsilon }_{A}}/t=-4$ , ${{t}^{{\prime\prime} }}/t=0.1$ . The shaded area represents the coexistence region.
Download figure:
Standard image High-resolution image

**Figure 8.** Schematic zero-temperature MF phase diagrams for spinless fermions on the (a) 1/3 filled lattice with ${{\varepsilon }_{A}}/t=4$ , ${{t}^{{\prime\prime} }}/t=-0.1$ and on the (b) 2/3 filled lattice with ${{\varepsilon }_{A}}/t=-4$ , ${{t}^{{\prime\prime} }}/t=0.1$ . The shaded area represents the coexistence region.
Download figure:
Standard image High-resolution image

Finally, with an intermediate value of V₂, there exists a narrow window for the coexistence of both QAH and nematic orders. This can be further seen in figure 9, showing the magnitude of both order parameters as a function of V₂ with fixed V₁. Since the bulk energy gap never really closes as V₂ increases, it suggests that the quantum phase transition between QAH and nematic phases is not continuous, lacking a quantum critical point. In addition, we also notice that there is no room for the current-loop states other than QAH state for systems with relatively large $|{{\varepsilon }_{A}}|$ and weak repulsions.

**Figure 9.** The magnitude of both QAH (blue, solid) and nematic (green, dashed) order parameters as a function of V₂ at 1/3-filled lattice with ${{\varepsilon }_{A}}/t=4$ , ${{t}^{{\prime\prime} }}/t=-0.1$ , and fixed ${{V}_{1}}/t=0.2$ .
Download figure:
Standard image High-resolution image

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 non-interacting energy spectrum for both fillings can be related by translating 'particle' into 'hole' language, which causes $t\to -t$ , ${{t}^{{\prime\prime} }}\to -{{t}^{{\prime\prime} }}$ , ${{\varepsilon }_{A}}\to -{{\varepsilon }_{A}}$ , and $\mu \to -\mu$ . The interactions given in the present form ruin such relation and hence the two phase diagrams must be different⁶ . 2) Notice that since in our consideration $|{{\varepsilon }_{A}}|$ is the largest energy scale among others, it is easy to realize that the charge density at A-site is $\delta n$ at 1/3 filling and $1-\delta n$ at 2/3 filling, where $\delta n$ denotes small density fluctuation. Such fact makes the NN repulsion (V₁ term) almost a constant depending on the total number of fermions, leading to V₁-insignificant phase diagrams in both cases. However, a close study in energetics (assuming the system is in QAH phase) can show that V₁ enters the dynamics through the first order of $\delta n$ for 1/3 filling while through the second order of $\delta n$ for 2/3 filling. Therefore, the phase diagram is relatively insensitive to V₁ in the 2/3-filling case. (3) When ${{\varepsilon }_{A}},{{t}^{{\prime\prime} }}\to 0$ , the spin-1 structure near band touching point, as discussed in the previous section, is recovered. Our MF study shows that the infinitesimal instability (near 1/3 or 2/3 filling) is absent due to the vanishing DOS of the dispersive bands and the semi-metallic phase is robust until V₂ reaches certain critical value. Moreover, for ${{V}_{2}}\geqslant {{V}_{2c}}$ , we find that the QAH phase only survives in a negligible window of V₂, and the nematic phase becomes the dominant one in the phase diagram (not shown). This result is similar to the work done by Liu et al [24] on the 2/3-filled kagome lattice with Dirac BCPs.

3.2. 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., $\mathcal{H}={{\mathcal{H}}_{0}}+{{\mathcal{H}}_{\operatorname{int}}}$ . The free part is again given by equation (1) with extra spin index σ in the fermion creation/annihilation operators; the interacting terms now contain

$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{\mathcal{H}}_{\operatorname{int}}} & = & \mathop{\sum }\limits_{\langle ij\rangle }\mathop{\sum }\limits_{\sigma ,{{\sigma }^{\prime }}}{{V}_{1}}{{n}_{i,\sigma }}{{n}_{j,{{\sigma }^{\prime }}}}+\mathop{\sum }\limits_{\langle \langle ij\rangle \rangle }\mathop{\sum }\limits_{\sigma ,{{\sigma }^{\prime }}}{{V}_{2}}{{n}_{i,\sigma }}{{n}_{j,{{\sigma }^{\prime }}}} \\ {} & {} & +\mathop{\sum }\limits_{i}{{U}_{i}}{{n}_{i,\uparrow }}{{n}_{i,\downarrow }}. \\ \end{array}\end{eqnarray} \tag{ 23 }$

In addition to V₁ (V₂) 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 MF approximation. The possible phases under translation-invariant ansatz are classified below⁷ :

(i)
Charge nematic state (CN). This phase is associated with spontaneously rotational (C₄) 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₂ terms, combined together to give out an effective NNN repulsion, $V_{2}^{\prime }=2({{V}_{2}}-\frac{U}{8})$ , playing similar role of V₂ in the spinless model.
(ii)
Nematic-spin-nematic state (NSN). This phase breaks C₄ 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₄ 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, ${{\vec{\eta }}^{t}}$ = $\frac{1}{16}{{\sum }_{{{\delta }^{\prime }}}}({{{\bf S}}_{Bi}}-{{{\bf S}}_{Ci+{{\delta }^{\prime }}}})$ , where ${{{\bf S}}_{\alpha }}$ = $\langle c_{\alpha i,\sigma }^{\dagger }{{{\bf s}}_{\sigma {{\sigma }^{\prime }}}}{{c}_{\alpha i,{{\sigma }^{\prime }}}}\rangle$ with $\alpha =A,B,C$ and ${\bf 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, $-\frac{2}{3}U{{({{\vec{S}}_{i}})}^{2}}$ with ${{\vec{S}}_{i}}$ representing usual spin operator. The other one is the bond-NSN, which is described by $\vec{Q}_{1}^{t}$ = $\frac{1}{8}\operatorname{Re}$ [ ${{\sum }_{\delta =\pm \hat{x}/2}}\langle c_{Ai,\sigma }^{\dagger }{{{\bf s}}_{\sigma {{\sigma }^{\prime }}}}{{c}_{Bi+\delta ,{{\sigma }^{\prime }}}}\rangle -{{\sum }_{\delta =\pm \hat{y}/2}}\langle c_{Ai,\sigma }^{\dagger }{{{\bf s}}_{\sigma {{\sigma }^{\prime }}}}{{c}_{Ci+\delta ,{{\sigma }^{\prime }}}}\rangle ]$ or $\vec{Q}_{2}^{t}=\frac{1}{8}{{\sum }_{{{\delta }^{\prime }}}}{{D}_{{{\delta }^{\prime }}}}\operatorname{Re}\langle c_{Bi,\sigma }^{\dagger }{{{\bf s}}_{\sigma {{\sigma }^{\prime }}}}{{c}_{Ci+\delta ,{{\sigma }^{\prime }}}}\rangle$ (breaking C₄ along a diagonal direction of the lattice).
(iii)
Charge current-loop state with broken TRS. As mentioned in the spinless model, among all the current patterns the case (a) in figure 7 is of most interest. This state, characterized by spontaneously broken TRS and parity symmetry, exhibits QAH effect with a spin-singlet order parameter, ${{\Phi }_{2}}=\frac{1}{8}{\rm Im}[{{\sum }_{{{\delta }^{\prime }},\sigma }}{{D}_{{{\delta }^{\prime }}}}\langle c_{Bi,\sigma }^{\dagger }{{c}_{Ci+{{\delta }^{\prime }},\sigma }}\rangle ]$ . Mainly driven by V₂ 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 current-loop states [e.g., from figures 7(b)–(d)] and find that none of them are 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 double-layer 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, or equivalently, 2D Z₂-non-trivial TI, with spin-triplet order parameter described by $\vec{\Phi }_{2}^{t}=\frac{1}{8}{{\sum }_{{{\delta }^{\prime }}}}{{D}_{{{\delta }^{\prime }}}}{\rm Im}[\langle c_{Bi,\sigma }^{\dagger }{{{\bf s}}_{\sigma {{\sigma }^{\prime }}}}{{c}_{Ci+{{\delta }^{\prime }},{{\sigma }^{\prime }}}}\rangle ]$ . Two remarks deserve mentioning here. First, both QAH and QSH are topological phases, characterized by a non-trivial topological index with robust edge states. However, the former one acquires a non-vanishing Chern number, while the latter one has zero Chern number due to TRS. Thus, a new topological index, called Z₂ index ν, needs to be introduced [4, 5, 9]. As detailed in appendix B, the QSH phase on the Lieb lattice indeed acquires non-trivial $\nu =1$ . Second, it is straightforward to see that at the MF level, the energy spectra (not shown here) for both QAH and QSH are the same. As a result, they have equal energy gain from V₂ repulsion and hence one cannot distinguish them in the MF phase diagram. If there were an extra NNN exchange coupling J₂ present in the system, the QAH would be favored for ${{J}_{2}}\gt 0$ ; reversely, the QSH would be favored for ${{J}_{2}}\lt 0$ due to its spin-triplet nature.

Under our assumption of translational invariance within our MF study, we do not consider any charge or spin density wave order. However, it is still worth pointing out that if $|{{\varepsilon }_{A}}|\to \infty$ 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 ${{\varepsilon }_{A}}\gt 0$ ( ${{\varepsilon }_{A}}\lt 0$ ), as guaranteed by the Lieb's theorem [44].

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 equation (22), where we have used the MF ansatz, $\langle {{n}_{Ai,\sigma }}\rangle =\frac{2}{3}+\rho +\frac{1}{2}\sigma {{S}_{Az}}$ , $\langle {{n}_{Bi,\sigma }}\rangle =\frac{2}{3}-\frac{\rho }{2}+\eta +\frac{1}{2}\sigma ({{S}_{Bz}}+2\eta _{z}^{t})$ , and $\langle {{n}_{Ci,\sigma }}\rangle =\frac{2}{3}-\frac{\rho }{2}-\eta +\frac{1}{2}\sigma ({{S}_{Bz}}-2\eta _{z}^{t})$ 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₀ in the spinful case becomes

$\begin{eqnarray}\begin{array}{ccccccccccccccc} \frac{{{E}_{0}}}{N} & = & -\frac{{{{\bar{\varepsilon }}}_{A1}}({{{\bar{\varepsilon }}}_{B1}}+{{{\bar{\varepsilon }}}_{C1}})}{{{V}_{1}}}-\frac{{{{\bar{\varepsilon }}}_{B2}}{{{\bar{\varepsilon }}}_{C2}}}{{{V}_{2}}}+8{{V}_{1}}(\delta {{t}^{2}}+Q_{1}^{2}+Q_{1z}^{t2}) \\ {} & {} & +8{{V}_{2}}(\delta {{t}^{{\prime} 2}}+Q_{2}^{2}+\Phi _{2}^{2}+Q_{2z}^{t2}+\Phi _{2z}^{t2}) \\ {} & {} & +\frac{U}{4}(S_{Az}^{2}+S_{Bz}^{2}+S_{Cz}^{2}+8\eta _{z}^{t2}) \\ {} & {} & -U[\frac{{{({{{\bar{\varepsilon }}}_{B1}}+{{{\bar{\varepsilon }}}_{C1}})}^{2}}}{16V_{1}^{2}}+\frac{\bar{\varepsilon }_{B2}^{2}+\bar{\varepsilon }_{C2}^{2}}{16V_{2}^{2}}], \\ \end{array}\end{eqnarray} \tag{ 24 }$

where ${{\bar{\varepsilon }}_{A1}},{{\bar{\varepsilon }}_{B1}},{{\bar{\varepsilon }}_{B2}},{{\bar{\varepsilon }}_{C1}}$ , and ${{\bar{\varepsilon }}_{C2}}$ are defined as before. ${{E}_{{\bf 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 MF U-V₂ phase diagram at 2/3 filling ( ${{\varepsilon }_{A}}=-4t$ ) for spinful fermions is shown in figure 10. 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₂ repulsion. In particular, the fact that the 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, the NSN phase begins to compete with the topological phase and clearly dominates over QAH/QSH whenever $U\gg {{V}_{2}}\gt 0$ . On the other hand, as $V_{2}^{\prime }=2({{V}_{2}}-\frac{U}{8})\gtrsim 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₁. As in the spinless case, at 2/3 filling the phase diagram is not sensitive to the presence of V₁. Especially, V₁ itself does not lead to any order. However, when V₁ becomes stronger, it is quite possible that the system might gain certain energy by opening a gap due to translation symmetry breaking, and hence is beyond our current consideration. The second point is about the system at 1/3 filling with ${{\varepsilon }_{A}}=4t$ . In fact, the phase diagram is qualitatively similar to that of 2/3 filling and therefore we omit it without further discussion.

**Figure 10.** Schematic zero-temperature mean-field phase diagram for spinful fermions at 2/3-filling with ${{\varepsilon }_{A}}/t=-4$ , ${{t}^{{\prime\prime} }}/t=0.1$ , and ${{V}_{1}}=0$ . The short-hand notations represent: NSN = nematic-spin-nematic; QAH = quantum anomalous Hall; CN = charge nematic. The shaded area denotes the coexistence region.
Download figure:
Standard image High-resolution image

4. 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 a 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 a suitable BCP 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 short-range 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 [39], 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 a Lieb lattice may come from layered perovskites. A well-known example is the CuO₂ plane in high-T_c cuprates such as La $_{1-x}$ Sr_xCuO₄ or YBa₂Cu₃O₇, 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 purposes, we take a typical three-band (Emery) model, written in hole language, having the same form of equation (1) [45, 46]. The corresponding model parameters are given by Hybertsen et al [48]: t = 1.5, ${{t}^{\prime }}=0.65$ , ${{\varepsilon }_{A}}=0$ , and ${{\varepsilon }_{B,C}}=3.6$ , where ${{t}^{\prime }}$ denotes NNN hopping. The band structure in the FBZ along high symmetry lines is presented in figure 11(a). As one can see, there is indeed a BCP at ${\bf 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}^{\prime }}$ . This change leads to the band structure shown in figure 11(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 Hamiltonian is similar to that given by Sun et al's recent work in appendix F [49]. 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.

**Figure 11.** (a) Energy spectrum of the model for cuprates along the high symmetry lines in the BZ with parameters mentioned in the text. (b) The same model as used in (a), but with reversed sign of ${{t}^{\prime }}$ .
Download figure:
Standard image High-resolution image

**Figure 11.** (a) Energy spectrum of the model for cuprates along the high symmetry lines in the BZ with parameters mentioned in the text. (b) The same model as used in (a), but with reversed sign of ${{t}^{\prime }}$ .
Download figure:
Standard image High-resolution image

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 ${{\varepsilon }_{A}}\gt 0$ and only NN hoppings in our spinless model. Now, the middle band is completely flat and touches with the lower dispersive band at ${\bf M}$ in FBZ. At 1/3 filling if we turn on V₂, the ground state would enter QAH state with an energy gap ${{\Delta }_{QAH}}$ opened at ${\bf M}$ . As we have mentioned in section 3.1, 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 $\sim \mathcal{O}({{\epsilon }_{A}}/{{\Delta }_{QAH}})$ . In other words, with appropriately chosen ${{\epsilon }_{A}}/{{\Delta }_{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 a fractional QAH state (or fractional Z₂ TI in the spinful case) when such band is partially filled [50–55].

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\pi$ , we have seen that in the T = 0 phase diagram the system can exhibit topological QAH/QSH, nematic, and NSN phases, depending on the strengths of short-range density–density interactions. In particular, the existence of the TI phase firmly justifies the interaction-driven scenario. Moreover, for a quadratic BCP (as ${{\epsilon }_{A}},{{t}^{{\prime\prime} }}\ne 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 nearly flat band, which is interesting in 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.

Acknowledgments

We thank L Fu, H Lin, C-K Lu, Y Ran, K Sun, and C Xu for helpful discussions and, especially, for K Sun's early collaboration in this work. Wei-Feng Tsai would like to thank for the hospitality of IoP, Chinese Academy of Sciences during his visit, where part of the present work was done. This work was supported in part by the MOST in Taiwan with Grant No. 103-2112-M-110-008-MY3 (WFT) and by the NSFC under Grant No.11474175 at Tsinghua (HY). Note added. After this work is completed, we were informed by S Kivelson that there is a partly related paper by M H Fischer and E-A Kim (arXiv:1106.6060).

Appendix A.: Perturbative RG analysis for a QBCP

In this section, we derive RG equations given in equation (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

$\begin{eqnarray}&&S={{S}_{0}}+{{S}_{\operatorname{int}}},\end{eqnarray} \tag{ A.1 }$

where the free action S₀ reads

$\begin{eqnarray}&&{{S}_{0}}=\mathop{\sum }\limits_{\sigma }\int \frac{{\rm d}\omega }{2\pi }\int \frac{{{{\rm d}}^{2}}p}{{{(2\pi )}^{2}}}{{\bar{\Psi }}_{\sigma }}(\omega ,{\bf p}){{[\hat{G}_{0}^{\sigma }(\omega ,{\bf p})]}^{-1}}{{\Psi }_{\sigma }}(\omega ,{\bf p}),\end{eqnarray} \tag{ A.2 }$

in momentum space. The Grassmann variables, ${{\bar{\Psi }}_{\sigma }}=({{\bar{\psi }}_{1\sigma }},{{\bar{\psi }}_{2\sigma }})$ 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 ${{[\hat{G}_{0}^{\sigma }(\omega ,{\bf p})]}^{-1}}=i\omega I-{{d}_{x}}{{\sigma }_{x}}-{{d}_{z}}{{\sigma }_{z}}$ , which is suitable for a particle-hole symmetric QBCP with ${{d}_{x}}={{p}_{x}}{{p}_{y}}/{{m}_{0}}$ , ${{d}_{z}}=(p_{x}^{2}-p_{y}^{2})\bar{t}/(2{{m}_{0}})$ , and $\bar{t}=2$ and is independent of spin polarization σ. We will fix $1/(2{{m}_{0}})$ as our unit in the analysis. The action of interactions can be written as,

$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{S}_{\operatorname{int}}} & = & {{S}_{1}}+{{S}_{2}}, \\ {{S}_{1}} & = & -\mathop{\sum }\limits_{a=1}^{2}\int {\rm d}\xi u{{{\bar{\psi }}}_{a\uparrow }}(4){{{\bar{\psi }}}_{a\downarrow }}(3){{\psi }_{a\downarrow }}(2){{\psi }_{a\uparrow }}(1), \\ {{S}_{2}} & = & -\mathop{\sum }\limits_{\sigma {{\sigma }^{\prime }}}\int {\rm d}\xi {{g}_{\sigma {{\sigma }^{\prime }}}}{{{\bar{\psi }}}_{1\sigma }}(4){{{\bar{\psi }}}_{2{{\sigma }^{\prime }}}}(3){{\psi }_{2{{\sigma }^{\prime }}}}(2){{\psi }_{1\sigma }}(1), \\ \end{array}\end{eqnarray} \tag{ A.3 }$

where $\int d\xi =\prod _{b=1}^{3}\int _{-\infty }^{\infty }\frac{d{{\omega }_{b}}}{2\pi }\int \frac{{{d}^{2}}{{p}_{b}}}{{{(2\pi )}^{2}}}$ and the shorthand notation $b=({{\omega }_{b}},{{{\bf 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}_{\sigma {{\sigma }^{\prime }}}}$ , 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 ( ${{{\bf p}}_{F}}=0$ ). This fact thus waives the complexity brought by the Fermi surface, a situation similar to ${{\phi }^{4}}$ -theory.

Let us now sketch how to obtain one-loop RG equations, following the standard perturbative RG procedure by Shankar [56]. Defining a momentum cutoff Λ and ${{\psi }_{\lt }}$ ( ${{\psi }_{\gt }}$ ) as slow (fast) modes with $0\lt p\lt \Lambda /s$ ( $\Lambda /s\lt p\lt \Lambda$ ), the key formula we use to derive the renormalized action by integrating out the fast modes is given by

$\begin{eqnarray}\begin{array}{ccccccccccccccc} Z & = & \int \mathcal{D}[{{{\bar{\psi }}}_{\gt }}{{\psi }_{\gt }};{{{\bar{\psi }}}_{\lt }}{{\psi }_{\lt }}]{{e}^{S[{{{\bar{\psi }}}_{\lt }}+{{{\bar{\psi }}}_{\gt }},{{\psi }_{\lt }}+{{\psi }_{\gt }}]}} \\ {} & = & \int \mathcal{D}[{{{\bar{\psi }}}_{\lt }}{{\psi }_{\lt }}]{{e}^{{{S}^{\prime }}[{{{\bar{\psi }}}_{\lt }}{{\psi }_{\lt }}]}}, \\ \end{array}\end{eqnarray} \tag{ A.4 }$

where

$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{e}^{{{S}^{\prime }}[{{{\bar{\psi }}}_{\lt }},{{\psi }_{\lt }}]}} & = & {{Z}_{0}}(\Lambda ,\Lambda /s){{e}^{{{S}_{0}}[{{{\bar{\psi }}}_{\lt }},{{\psi }_{\lt }}]}}{{\langle {{e}^{{{S}_{\operatorname{int}}}}}\rangle }_{0,\gt }} \\ {} & = & {{Z}_{0}}(\Lambda ,\Lambda /s){{e}^{{{S}_{0}}[{{{\bar{\psi }}}_{\lt }},{{\psi }_{\lt }}]}} \\ {} & {} & \times {{e}^{[{{\langle {{S}_{\operatorname{int}}}\rangle }_{0,\gt }}+\frac{1}{2}({{\langle S_{\operatorname{int}}^{2}\rangle }_{0,\gt }}-\langle {{S}_{\operatorname{int}}}\rangle _{0,\gt }^{2})+\cdots ]}}. \\ \end{array}\end{eqnarray} \tag{ A.5 }$

Note that in deriving the last equality we have used the cumulant expansion. Z₀ denotes the non-interacting partition function and ${{\langle \quad \rangle }_{0,\gt }}$ represents average with respect to the fast modes of action S₀. The collected exponent of the exponential terms now results in the renormalized action ${{S}^{\prime }}$ .

Before examining the renormalization of coupling parameters u and ${{g}_{\sigma {{\sigma }^{\prime }}}}$ , we first notice that, to the one-loop level, the quartic terms in S_int do induce quadratic terms from the tadpole diagram shown in figure 12, which are momentum independent. This would indicate there is no non-interacting FP in our system. However, explicit calculations show that all contributions from such types of diagram are zero and hence our starting FP survives without flowing ( $\bar{t}$ , ${{\psi }_{a\sigma }}$ are not renormalized), in contrast to the case of the Luttinger liquid.

**Figure 12.** The tadpole diagram, induced by the quartic terms, can in principle renormalize the free action. The arrowed lines on the left represent the intra-orbital (solid) and inter-orbital (dashed) non-interacting Green's functions, respectively.
Download figure:
Standard image High-resolution image

Next, we turn our attention to the diagrams that renormalize coupling parameters. As shown in figures 13(a)–(c), they are, respectively, total contributions to the renormalized u, ${{g}_{\sigma \sigma }}$ , and ${{g}_{\sigma \bar{\sigma }}}$ . Obviously, they have standard structures such as 'BCS', 'ZS' (zero sound), and 'ZS $^{\prime }$ ' used in Shankar's seminal work [56]. Note that we have taken the convention that $\bar{1}=2,\bar{\uparrow }=\downarrow$ and vice versa. The essential fermion loop integrals are listed below (without vertex):

$\begin{eqnarray}\begin{array}{ccccccccccccccc} {\rm BCS}1 & : & {{\int }_{q}}G_{0aa}^{\sigma }(q)G_{0\bar{a}\bar{a}}^{{{\sigma }^{\prime }}}(-q)={\rm d}l\cdot \gamma (\bar{t}), \\ {\rm BCS}2 & : & {{\int }_{q}}G_{0a\bar{a}}^{\sigma }(q)G_{0\bar{a}a}^{{{\sigma }^{\prime }}}(-q)={\rm d}l\cdot \gamma (\bar{t}), \\ {\rm BCS}3 & : & {{\int }_{q}}G_{0aa}^{\sigma }(q)G_{0aa}^{{\bar{\sigma }}}(-q)={\rm d}l\cdot [\alpha (\bar{t})-\gamma (\bar{t})], \\ {\rm ZS}1 & : & {{\int }_{q}}G_{0aa}^{\sigma }(q)G_{0\bar{a}\bar{a}}^{{{\sigma }^{\prime }}}(q)=-{\rm d}l\cdot [\alpha (\bar{t})-\gamma (\bar{t})], \\ {\rm ZS}2 & : & {{\int }_{q}}G_{0a\bar{a}}^{\sigma }(q)G_{0\bar{a}a}^{\sigma }(q)={\rm d}l\cdot \gamma (\bar{t}), \\ {\rm ZS}3 & : & {{\int }_{q}}G_{0aa}^{\sigma }(q)G_{0aa}^{{{\sigma }^{\prime }}}(q)=-{\rm d}l\cdot \gamma (\bar{t}), \\ \end{array}\end{eqnarray} \tag{ A.6 }$

where ${{\int }_{q}}\equiv \int _{-\infty }^{\infty }\frac{{\rm d}\omega }{2\pi }\int _{0}^{2\pi }\frac{{\rm d}\theta }{2\pi }\int \Lambda /{{s}^{\Lambda }}\frac{{\rm d}q}{2\pi }$ and $\bar{t}$ , $\alpha (\bar{t})$ , $\gamma (\bar{t})$ are defined in the main text [just below equation (13)]. The superscript, ${{\sigma }^{\prime }}$ , of G₀ can be either σ or $\bar{\sigma }$ . The presence of dl is due to the approximation we have used, ${\rm ln} s={\rm ln} (1+dl)\approx dl$ in the limit $s\to 1$ . In order to evaluate these integrals, we have set all external $\omega =p=0$ and made ${{q}_{x}}=q{\rm cos} \theta ,{{q}_{y}}=q{\rm sin} \theta$ . Combining all the Feynman diagrams in figure 13 and keeping in mind the order of fermion operators, the straightforward algebra gives rise a set of RG equations, as we desire, equation (13).

**Figure 13.** Diagrams to one-loop level that contribute to renormalize (a) u, (b) ${{g}_{\sigma \sigma }}$ , and (c) ${{g}_{\sigma \bar{\sigma }}}$ . Different types of vertices are shown on the left. The labels for external legs in each one-loop diagram are the same as the vertex on the left and hence are omitted.
Download figure:
Standard image High-resolution image

Finally, we investigate the RG flow beginning with weak coupling regime for two typical cases: (a) $u\gt {{g}_{\sigma \sigma }}={{g}_{\sigma \bar{\sigma }}}$ ; (b) $u\lt {{g}_{\sigma \sigma }}={{g}_{\sigma \bar{\sigma }}}$ . As one can see in figure 14, in either case, they all flow to strong coupling regime and hence the non-interacting FP is actually unstable against short-range repulsions.

**Figure 14.** RG flows for u(l) (solid, blue), ${{g}_{\uparrow \uparrow }}(l)$ (dotted,black), and ${{g}_{\uparrow \downarrow }}(l)$ (dashed, red) with interaction profile (a) $u(0)=0.5,{{g}_{\uparrow \uparrow }}(0)={{g}_{\uparrow \downarrow }}(0)=0.05$ ; (b) $u(0)=0.05,{{g}_{\uparrow \uparrow }}(0)={{g}_{\uparrow \downarrow }}(0)=0.5$ . In both cases, $\bar{t}=2$ .
Download figure:
Standard image High-resolution image

**Figure 14.** RG flows for u(l) (solid, blue), ${{g}_{\uparrow \uparrow }}(l)$ (dotted,black), and ${{g}_{\uparrow \downarrow }}(l)$ (dashed, red) with interaction profile (a) $u(0)=0.5,{{g}_{\uparrow \uparrow }}(0)={{g}_{\uparrow \downarrow }}(0)=0.05$ ; (b) $u(0)=0.05,{{g}_{\uparrow \uparrow }}(0)={{g}_{\uparrow \downarrow }}(0)=0.5$ . In both cases, $\bar{t}=2$ .
Download figure:
Standard image High-resolution image

Appendix B.: The proof of non-trivial Z₂ index for the quantum spin Hall phase

The existence of the QSH state can be justified by a non-trivial Z₂-valued invariant $\nu =1$ . This Z₂ 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₂ index explicitly. This method essentially takes the advantage of IS such that the index can be related to the parity eigenvalues ${{\xi }_{2\,m}}({{\Gamma }_{i}})$ of $2\;m$ occupied energy bands at four TRIM: ${{\Gamma }_{i}}={\boldsymbol{ \Gamma }} ,{\bf X},{\bf Y}$ , and ${\bf M}$ .

For the purpose of instruction, we focus on 1/3 filling with ${{\varepsilon }_{A,B,C}}=0$ and hence only the lower band is relevant. As shown in section 3, the presence of V₂ can dynamically generate the 'spin–orbit' interaction, which is equivalent to add the following term into equation (2).

$\begin{eqnarray}&&{{\mathcal{H}}_{SO}}=\mathop{\sum }\limits_{{\bf k},\sigma }\psi _{{\bf k}\sigma }^{\dagger }{{H}_{SO}}({\bf k}){{\psi }_{{\bf k}\sigma }}\end{eqnarray} \tag{ B.1 }$

with ${{H}_{SO}}({\bf k})=$

$\begin{eqnarray}\pm {{\lambda }_{SO}}\left( \begin{array}{ccccccccccccccc} 0 & 0 & 0 \\ 0 & 0 & 2i{{w}_{{\bf k}}} \\ 0 & -2i{{w}_{{\bf k}}} & 0 \\ \end{array} \right),\end{eqnarray} \tag{ B.2 }$

where ${{w}_{{\bf k}}}\equiv \sum _{j=1}^{2}{{(-1)}^{j+1}}{\rm cos} ({\bf k}\cdot {\bf a}_{j}^{\prime })$ with ${\bf a}_{1}^{\prime }={{{\bf a}}_{1}}+{{{\bf a}}_{2}}$ , and ${\bf a}_{2}^{\prime }={{{\bf a}}_{2}}-{{{\bf a}}_{1}}$ . The $+(-)$ sign refers to spin-up (down) fermions. The eigenvalues of $H({\bf k})={{H}_{0}}({\bf k})+{{H}_{SO}}({\bf k})$ are ${{\epsilon }_{\pm }}({\bf k})=\pm \sqrt{{{b}_{{\bf k}}}+4\lambda _{SO}^{2}w_{{\bf k}}^{2}}$ and ${{\epsilon }_{0}}({\bf k})=0$ and their corresponding eigenstates can be written in the form,

$\begin{eqnarray}&&|{{u}_{n{\bf k}}}\rangle ={{G}_{n{\bf k}}}{{\left( {{q}_{A{\bf k}}},{{q}_{B{\bf k}}},{{q}_{C{\bf k}}} \right)}^{{\rm T}}},\end{eqnarray} \tag{ B.3 }$

where the expressions of the components ${{q}_{l{\bf k}}}$ ( $l=A,B,C$ ) and the normalization factor ${{G}_{n{\bf k}}}$ for each band n( = 0, ±) are given in table 1.

Table 1. The expressions for the coefficients of $|{{u}_{n{\bf k}}}\rangle$ with x_i = ${\bf k}\cdot {{{\bf a}}_{i}}$ (i = 1, 2). For ${\bf k}={\bf M}$ , the listed expressions break down with ambiguity and hence we single them out, $|{{u}_{+{\bf M}}}\rangle ={{(0,1/\sqrt{2},i/\sqrt{2})}^{T}}$ , $|{{u}_{-{\bf M}}}\rangle ={{(0,1/\sqrt{2},-i/\sqrt{2})}^{T}}$ and $|{{u}_{0{\bf M}}}\rangle ={{(1,0,0)}^{T}}$ . Note that the upper (lower) sign in the square bracket refers to spin-up (down) fermions. ${{b}_{{\bf k}}}=\sum _{i=1}^{2}{{[2t{\rm cos} ({\bf k}\cdot ,{{a}_{i}})]}^{2}}$ .

${{q}_{A{\bf k}}}$	$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \epsilon _{n}^{2}({\bf k})-4\lambda _{SO}^{2}w_{{\bf k}}^{2}$
${{q}_{B{\bf k}}}$	$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 2t[{{\epsilon }_{n}}({\bf k}){\rm cos} {{x}_{1}}\pm 2i{{\lambda }_{SO}}{{w}_{{\bf k}}}{\rm cos} {{x}_{2}}]$
q_Ck	2t[i_n(k)cos x₂∓ 2iλ_SOw_kcosx₁]
G_nk^-2	_n²k)[_n²(k)-4λ_SO²w_k²+b_k]
	+4λ_SO²w_k²(b_k+4λ_SO²w_k²)

At TRIM, we have eigenstates $|{{u}_{-{\boldsymbol{ \Gamma }} }}\rangle ={{(\frac{1}{\sqrt{2}},\frac{-1}{2},\frac{-1}{2})}^{T}}$ , $|{{u}_{-{\bf X}}}\rangle ={{(\frac{1}{\sqrt{2}},0,\frac{1}{\sqrt{2}})}^{T}}$ , $|{{u}_{-{\bf Y}}}\rangle ={{(\frac{1}{\sqrt{2}},\frac{-1}{\sqrt{2}},0)}^{T}}$ , and $|{{u}_{-{\bf M}}}\rangle ={{(0,\frac{1}{\sqrt{2}},\frac{-i}{\sqrt{2}})}^{T}}$ according to table 1, respectively. Now, let us make a gauge transformation on eigenstates such that only one coordinate ${\bf R}$ would be assigned to each unit cell (instead of ${\bf R}+{{{\bf a}}_{i}}$ for each atom within the cell), and the Bloch Hamiltonian $H({\bf k})$ now obeys $H({\bf k}+{\bf G})=H({\bf k})$ and, with the inversion operator $\hat{P}$ , $H(-{\bf k})=\hat{P}H({\bf k})\hat{P}$ , resulting in $[H({\bf k}={{\Gamma }_{i}}),\hat{P}]=0$ . Thus, the parity eigenvalue now becomes a good quantum number at TRIM. Under such transformation the eigenstates are rewritten as $|{{u}_{-{\boldsymbol{ \Gamma }} }}\rangle ={{(\frac{1}{\sqrt{2}},\frac{-1}{2},\frac{-1}{2})}^{T}}$ , $|{{u}_{-{\bf X}}}\rangle ={{(\frac{1}{\sqrt{2}},0,\frac{1}{\sqrt{2}})}^{T}}$ , $|{{u}_{-{\bf Y}}}\rangle ={{(\frac{1}{\sqrt{2}},\frac{-1}{\sqrt{2}},0)}^{T}}$ , and $|{{u}_{-{\bf M}}}\rangle ={{(0,\frac{-i}{\sqrt{2}},\frac{-1}{\sqrt{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 ${\bf R}$ is simply multiplying our eigenstate at ${{\Gamma }_{i}}$ by a factor ${{e}^{i{{\Gamma }_{i}}\cdot {\bf R}}}$ , we now see the parity eigenvalues for ${{\Gamma }_{i}}$ = $({\boldsymbol{ \Gamma }} ,{\bf X},{\bf Y},{\bf M})$ , are, respectively, P = (+,+,+,-). The Z₂ index ν is then determined by

$\begin{eqnarray}&&{{(-1)}^{\nu }}=\mathop{\prod }\limits_{i=1}^{4}{{\delta }_{i}}=-1,\end{eqnarray} \tag{ B.4 }$

where ${{\delta }_{i}}=\prod _{m=1}^{N}{{\xi }_{2m}}({{\Gamma }_{i}})$ , and N = 1 indicating either spin-up or spin-down band in our example here. Therefore, $\nu =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 2. Finally, for ${{\varepsilon }_{A}}\ne 0$ , one can follow similar procedure and the conclusion of $\nu =1$ is unchanged.