Tuning topological phase and quantum anomalous Hall effect by interaction in quadratic band touching systems

Interaction-driven topological phases significantly enrich the class of topological materials and thus are of great importance. Here, we study the phase diagram of interacting spinless fermions filling the two-dimensional checkerboard lattice with a quadratic band touching (QBT) point. By developing new diagnosis based on the state-of-the-art density-matrix renormalization group and exact diagonalization, we determine accurate quantum phase diagram for such a system at half-filling with three distinct phases. For weak nearest-neighboring interactions, we demonstrate the instability of the QBT towards an interaction-driven spontaneous quantum anomalous Hall (QAH) effect. For strong interactions, the system breaks the rotational symmetry realizing a nematic charge-density-wave (CDW) phase. Interestingly, for intermediate interactions we discover a symmetry-broken bond-ordered critical phase sandwiched in between the QAH and CDW phases, which splits the QBT into two Dirac points driven by interaction. Instead of the direct transition between QAH and CDW phases, our identification of an intermediate phase sheds new light on the theoretical understanding of the interaction-driven phases in QBT systems. Topological materials are attracting a great deal of attention, and proposals suggesting that topologically trivial materials with a quadratic band touching point can be driven in a topological phase by electronic interactions that break time-reversal symmetry are very appealing. However, because of the very small energy gap protecting such topological phase, this effect remained so far elusive in microscopic calculations. Tian-Sheng Zeng, Wei Zhu and Donna Sheng developed a systematic and controlled scheme that allowed them to determine the accurate phase diagram of a 2D checkboard lattice with a quadratic band touching point, finding that for weak interactions a topological phase emerges. At stronger interactions an intermediate critical phase is observed before a nematic charge-density wave phase develops. This work opens the way to a deeper understanding of other topological systems driven by electronic interactions.


INTRODUCTION
Recently topological phases of matter in the band structure with nontrivial topological Berry phase become an exciting research area of modern physics, culminating in the experimental observations of Haldane-honeycomb insulator 1 and the quantum anomalous Hall (QAH) effect in topological insulator. 2 Generally speaking, in such kind of systems with a nontrivial topological invariant Chern number, 3 exemplified by the integer quantum Hall effect in the absence of magnetic field, 4 time-reversal symmetry (TRS) breaking of band structure is typically necessary. For topologically trivial band structures with zero Chern number [5][6][7][8] , it was proposed that the strong correlation between electrons can also induce spontaneous TRS breaking at mean field level 9,10 and lead to QAH effect in Dirac semimetals. These interaction-induced topological phases are interesting, since they can significantly enrich the class of topological materials. [11][12][13] However such a mechanism is challenged by subsequent numerical simulations, [14][15][16][17] which do not support the interaction-driven QAH effect. Alternatively, it has been identified that a quadratic band touching (QBT) point has topological feature, and can be driven towards a QAH phase with TRS breaking even under arbitrary weak interactions, while strong interactions may lead to other competing phases. 18,19 The theoretical prediction of such interaction-driven topological phases, [20][21][22][23][24][25][26][27][28] stimulates extensive studies by more rigorously theoretical and numerical methods, including low energy renormalization group approach in C 4 symmetric checkerboard lattice [29][30][31][32] and bilayer graphene, 33,34 first-principle calculations in spin-dependent optical square lattice 35 and halogenated hematite nanosheets, 36 and recently the unbiased numerical exact diagonalization (ED) diagnosis in both C 4 symmetric checkerboard lattice 37 and C 6 symmetric Kagome lattice. 38 However, the stability of the QAH effect in the presence of weak interaction has not been settled. Based on the mean-field theory, the gap protecting the QAH state is exponentially small for weak interaction, making such a phase difficult to be identified. A recent density-matrix renormalization group (DMRG) study finds a semimetal phase for Kagome QBT systems with nearest-neighboring interactions, 39 while adding further long-ranged interactions (V 1~V2~V3 ) favors a robust gapped QAH state. 38 Thus, the weak interaction effect on the QBT point remains elusive, because finite-size calculations are incapable of detecting extremely small energy gap. Furthermore, the issue of the quantum phase transition between the interaction-driven topological phase and other phases is hardly touched. While ED analysis points to a first-order phase transition from the QAH phase to the charge density wave phase without intermediate phase, 37,38,40 the interesting scenario of intermediate Dirac liquid phase in such QBT systems deserves to be explored by applying controlled numerical methods.
In this work, we develop accurate numerical diagnosis for the topological phases to address these challenge issues through the state-of-the-art DMRG and ED simulations. A schematic diagram of our main results is shown in Fig. 1 for interacting spinless fermions occupying the topologically trivial checkerboard lattice with a QBT. For weak interactions V < V 1 , by applying the Hellmann-Feynman theorem, we demonstrate the emergence of interactiondriven QAH phase whose topological properties are featured by two degenerate TRS breaking ground states arising from a pair of Kramers degenerate states with opposite chiralities and integer quantized topological Hall conductances. Remarkably, for intermediate interactions V 1 < V < V 2 we establish the existence of a new gapless critical phase as a sublattice bond-ordered nematic phase, while the nematic CDW phase appears in the strongly interacting regime. Our comprehensive DMRG and ED studies can access large system sizes to establish the stability of the QAH state with weak interaction at the thermodynamic limit, and illustrate the interaction controlling of the anomalous dissipationless Hall transport properties of the QBT systems through opening and closing of the energy excitation gap.

Hamiltonian
We consider the spinless fermions in the topologically trivial checkerboard lattice model, Here a r , b r′ are the particle annihilation operators for sublattices A, B respectively, and n a r ¼ a y r a r ; n b r ¼ b y r b r the particle number operators at site r. 〈…〉 and 〈〈…〉〉 denote the nearest-neighbor and the next-nearest-neighbor pairs of sites on a checkerboard lattice, respectively. As shown in Fig. 1, the hopping amplitudes t a r;r 0 ¼ t 0 for solid lines and t a r;r 0 ¼ Àt 0 for dashed lines in sublattice A, while t b r;r 0 ¼ t 0 for solid lines and t b r;r 0 ¼ Àt 0 for dashed lines in sublattice B. The single-particle dispersion hosts a QBT at (k x , k y ) = (π, π) with Berry flux ±2π, protected by TRS and C 4 rotational symmetry. Note that by adding a nonzero opposite shift δ ≠ 0 to the hopping t aðbÞ r;r 0 ¼ t aðbÞ r;r 0 þ ðÀÞδ of both sublattices, the system H breaks C 4 symmetry down to C 2 symmetry, and the QBT splits into two Dirac points with gapless particle-hole symmetric dispersions.

QAH phase
We first address the critical issue if a robust QAH effect can be stabilized with the presence of weak and nearest interaction. We begin with studying the ground state properties up to a maximum system sizes N s = 32 based on ED. The geometry hosts both C 4 point-group symmetry and TRS, and we find an exact twofold ground state degeneracy ψ ± j i at momentum K = (0, 0) below excitation continuum for systems with finite interaction V. This pair of degenerate eigenstates are also eigenstates of C 4 rotation with eigenvalues ±i, respectively, which serves as important evidence of TRS breaking as long as the energy excitation gap remains open for large systems. In order to illustrate two TRS spontaneously breaking states with opposite chiralities, we consider the system response to the TRS breaking perturbation hopping phase e iϕ a y r b r 0 þ h:c: to the nearest-neighbor A and B sites (the phase ϕ is a tiny detecting flux per plaquette for detecting QAH order). From the Hellmann-Feynman theorem, 41,42 we can derive the TRS breaking chiral bond current J r ¼ i a y r b r 0 À b y r 0 a r between nearest-neighboring sites from the linear response of the ground state energy, as 067, implying the opposite chiralities of TRS breaking for weak interactions. To extract the topological invariants of the doublet ground states for any value ϕ, we utilize the twisted boundary conditions ψðϕ; r þ N αêα Þ = ψ (ϕ;r) exp(iθ α ) where θ α is the twisted angle in the α (x or y)direction. The system is periodic when one flux quantum θ α = 0 → 2π is inserted. Meanwhile, the many-body Chern number of the ground state wavefunction ψ ± (ϕ) is defined as 43,44 We identify the topological invariants C ± (ϕ) = ±1 for the twofold ground states ψ ± ðϕÞ j iunder an infinitesimal TRS breaking phase ϕ ( 1. Due to the adiabatic connection between the ψ ± (ϕ) and ψ ± (0) as shown in Fig. 2a, we can obtain the topological invariants for these doublet states C ± = C ± (ϕ → 0) = ±1. On the contrary, for strong interactions, we find that the expectation value of J r in the  ground state vanishes precisely, and the topological invariants C ± (ϕ) = 0, signalling a topologically trivial nematic CDW phase. Indeed, the density structure factor To access larger system sizes to establish the stability of the QAH, we exploit an unbiased DMRG approach using a cylindrical geometry up to a maximum width L y = 16 (N y = 8). By randomly choosing different initial states in DMRG simulations, we can obtain two different ground states ψ ± j i with degenerate energies E + ≃ E − and opposite chiral circulating loop currents per plaquette ψ ± h jJ r ψ ± j i, from a DMRG algorithm allowing complex wavefunctions. In Fig. 2b, we measure the current-current correlation functions ψ þ J r J r0 ψ þ between nearest-neighboring bonds r; r 0 h i and r 0 ; r 0 0 , with the distance r À r 0 j j. For different system sizes, the bond current long-range order parameter ψ þ J r ψ þ = lim ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ψ þ J r J r0 ψ þ q for weak interactions persists to a finite value at the large distance r À r 0 j j limit. Meanwhile in Fig. 2b, we also characterize the topological nature of the ground state from its topological charge pumping by inserting one U(1) charge flux quantum θ y = θ from θ = 0 to θ = 2π in the periodic y-direction of the cylinder system based on the adiabatic DMRG 45 in connection to the quantized Hall conductance. Here we partition the lattice system on the cylinder along the x-direction into two halves with equal lattice sites. The transverse transfer of the total charge from the right side to the left side in the x-direction is encoded by the expectation value QðθÞ ¼ trρ L ðθÞN L Â Ã . N L is the particle number in the left cylinder part, andρ L the reduced density matrix of the corresponding left part. Under the inserting of the flux θ y = θ in the y-direction, the change of Q(θ) indicates the transverse charge transfer from the right side to the left side in the x-direction, induced by the topological Hall conductances of the state ψ ± j i. From Fig. 2c, we obtain a nearly quantized transverse Hall conductance C ± = ΔQ = Q(2π) − Q(0) = ±1 for these two degenerate ground states ψ ± j i. Again, for strong interactions, the current-current correlation functions are vanishingly small, and the charge pumping disappears, signaling the absence of a QAH phase. Significantly, our results establish that any weak interaction would drive the system into the QAH phase.
Intermediate phase We turn to analyze the emergent intermediate phase between the QAH phase and the nematic CDW phase. Figure 3a1 depicts the evolution of the interacting many-body low energy spectrum. Near the point V = V 1 , the doublet ground states of the QAH phase undergo a level crossing with the twofold degenerate excited levels. When the interaction V increases further, there appears another level crossing near the transition around V = V 2 . For V 1 < V < V 2 , the doublet ground states host nonzero bond orders Δ r ≠ 0, which is calculated from the symmetry-breaking perturbation response Δ r = 1 2Ns ∂EðδÞ ∂δ δ¼0 = a y r a rþêx À a y r a rþêy = b y r b rþêx À b y r b rþêy (The bonds a y r a rþêx and a y r a rþêy have opposite signs and different magnitudes. Also, a y r a rþêx = Àb y r b rþêx and a y r a rþêy = Àb y r b rþêy ). For example, we have Δ r~0 .08 at V = 2.15, as shown in Fig. 3a2. The nonzero bond order indicates the bond nematic nature of the intermediate phase.
Meanwhile, we plot the evolutions of the current order parameter J r h i, sublattice bond order parameter Δ r h i, nematic CDW density structure factor S/N s and the charge-hole gap Δ c = (E 0 (N + 1) + E 0 (N − 1) − 2E 0 (N))/2 as a function of V in Fig. 3b. For weak interactions V < V 1 , J r h i has a finite expectation value, and Δ c increases with the increase of the interaction strength, manifesting the robustness of a gapped QAH phase. On the other hand, both Δ r h i and S/N s take small values consistent with the properties of a gapped QAH phase. When interaction V goes across V 1 , J r h i, Δ r h i, and S/N s experience a sudden jump, where J r h i drops down to a vanishingly small value of the order 10 −4 , and Δ c begins to decrease quickly, signaling the collapsing of a topological phase (finite size scaling results for the Δ c will be discussed below). Simultaneously, Δ r h i, and S/N s jump to a finite large value. When interaction further increases beyond the critical value V 2 , Δ r h i drops down to a vanishingly small value of the order 10 −5 , while S/N s undergoes another step jump to a larger value, where the system becomes an insulating nematic CDW phase with a large excitation gap as shown in Fig. 3a1.
As presented above, based on a given system size calculation, we see the intermediate phase between V 1 and V 2 is associated with a finite bond order Δ r h i. To inspect properties of the phase in the thermodynamic limit, we carry out a finite size scaling of bond order Δ r and nematic order δ CDW = n a r À n b r 0 = 2 Ns ∂EðmÞ ∂m . In Fig. 4, it is found that Δ r extrapolates to a finite value, but δ CDW gradually decreases down to a negligibly small value in the limit 1=N 2 y ! 0. Thus, our ED and DMRG methods provide a very strong evidence for the existence of the intermediate topologically trivial bondordered phase. Physically, the doublet states of the QAH phase host opposite eigenvalues ±1 angular momentum in the presence of C 4 rotation symmetry and TRS, while both bond-ordered and nematic CDW phase breaks C 4 symmetry down to C 2 symmetry. In ref. 37 , such an argument is used to claim a first-order phase transition between C 4 symmetric phases and C 2 symmetric phases. Our ED and DMRG study access much larger systems, which lead to the discovery of a time-reversal symmetric intermediate bondordered phase sandwiched in between QAH and CDW phases (The maximal system size in our ED calculations is 2 × 4 × 4 = 32, which goes beyond the limit in ref. 37 . The 2 × 4 × 4 system respects both C 4 symmetry and C 2 symmetry, and hosts the QBT point (π, π), which is unique and reduces finite-size effect. Our DMRG study on larger system sizes also confirms the intermediate phase and the two-step phase transitions, in agreement with our ED study of 2 × 4 × 4 system size). At the mean-field level, 18 one can show that this phase is gapless with the QBT splitting into two Dirac points (see Fig. 1). Moreover, we further provide strong numerical evidences to support the gapless nature of the intermediate phase, which point to the existence of Dirac cone structure. First, we find that the finite-size scaling of the chargehole gap Δ c in Fig. 5a gives a nearly zero value in the thermodynamic limit, which implies the gapless single-particle excitation nature. Second, the entropy dependence on the cylinder length approaches an arch structure as indicated in Fig.  5b, and it can be fitted to the universal scaling function (up to an additive constant depending on the cylinder width) 46-48 SðxÞ = c 3 log Lx π sin πx Lx with the central charge c ≈ 2, as plotted in Fig. 5c. The only deviation from the straight line fitting shown in Fig. 5c appears at larger x values due to the convergence difficulty in capturing the entanglement of a gapless system. The finite gap for finite size systems shown in Fig. 5a can be understood as the following. On the finite size systems, the available momentum points in the Brillouin zone are discrete, such that the Dirac points are not guaranteed to be exactly covered in our calculations and a finite-size gap appears. Similarly, when experienced a small energy gap, the central charge obtained from entropy behavior would deviate from the ideal value c = 2 of free Dirac fermions. Nevertheless, the finite-size scaling of the gap, and the entropy scaling approaches the Dirac liquid behavior in the thermodynamic limit, as indicated in Fig. 5a, c. These numerical results indeed support that both the charge-hole gap Δ c → 0 and the central charge c → 2 as the cylinder width L y = 2N y increases.

Phase transition
To further study phase transitions, we exploit the finite DMRG calculation on a cylindrical geometry up to a maximum width L y = 16 (N y = 8) and length L x = N x = 20. We measure five different j (δV is as small as 0.1t), the entanglement entropy S L in the middle of the cylinder, in addition to order parameters Δ r , δ CDW and J r . The first-order transition is characterized by the discontinuous behavior of these physical quantities.
As shown in Fig. 6a-c, for V < V 1 , F(V) has a large value close to 1, both S L and δ CDW exhibit featureless properties, but J r grows slightly with the increase of V, implying the robustness of the QAH phase. When V approaches a transition point V 1 , F(V) suddenly drops down to a very small value close to zero, and both order parameters Δ r , J r exhibit a sharp discontinuous jump near the transition point. Similarly, S L starts to drop at V ≃ V 1 , with a discontinuous derivative ∂S L /∂V. These results are consistent with a first order transition between the QAH and intermediate phases. 49,50 When the interaction V increases further close to V 2 , δ CDW undergoes a jump to a large saturated value close to 1, while Δ r gradually drops down to zero in the strongly interacting regime. F(V) exhibits a minimum, sharpening as the cylinder width increases from N y = 4 to N y = 8, while entanglement entropy S L also drops down to a smaller value consistent with an insulating phase at V > V 2 . These results are consistent with a first-order transition into a CDW phase at V = V 2 .

DISCUSSION AND SUMMARY
In summary, we have numerically presented a solid diagnosis of an interaction-driven spontaneous QAH phase in the checkerboard lattice with a QBT by turning on any weak interaction. Such a diagnosis relies on finite size scaling up to wide systems (L y = 20 lattice spacing) and uses detecting flux for the QAH. The QAH phase hosts twofold ground state degeneracies with opposite spontaneous TRS breaking behaviors, and a quantized Hall conductance measured by Laughlin argument of charge pumping. In particular, we demonstrate the existence of a bond-ordered phase sandwiched between the QAH phase and the nematic CDW phase, characterized by the sublattice bond order. The intermediate bond-ordered phase points to the long-sought Dirac liquid phase in the QBT systems. 18,34 We believe that this work would open a new route for the study of the possible competing intermediate phases under the interplay of interaction and frustration, and excite a more extensive investigation of the fate of the QAH phases in many other systems, such as bilayer graphene, 33,34 and C 6 symmetric Kagome lattice 20 where an intermediate gapless CDW phase was also proposed. Other future directions include a study of the interplay of nearest-neighboring and next-nearest-neighboring interactions, which may lead to rich possibility and other competing phases driven by interactions. 51 At experimental side, our work also suggests a practical way of opening and closing the QBT gap, inducing quantized dissipationless transport currents by interactions.
In the preparation of this work, we become aware of a parallel work, ref. 51 .

Exact diagonalization
We perform ED calculations on the model Eq. (1) with parameters as indicated in the corresponding text and figure captions. In the ED calculations, we study the many-body ground state of H at half-filling in a finite system of N x × N y unit cells (the total number of sites N s = 2 × N x × N y ). The energy eigenstates are labeled by a total momentum K = (K x , K y ) in units of (2π/N x , 2π/N y ) in the Brillouin zone.

Density-matrix renormalization group
For larger systems we exploit both finite and infinite DMRG on the cylindrical geometry. We keep the dimension of DMRG kept states up to 5600 to obtain accurate results (the truncation error is of the order 10 −6 ). This leads to excellent convergence for the results that we report here. The geometry of cylinders is open boundary condition in the x-direction and periodic boundary condition in the y-direction.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.