Accessing the topological Mott insulator in cold atom quantum simulators with realistic Rydberg dressing

The interplay between many-body interactions and the kinetic energy gives rise to rich phase diagrams hosting, among others, interaction-induced topological phases. These phases are characterized by both a local order parameter and a global topological invariant, and can exhibit exotic ground states such as self-trapped polarons and interaction-induced edge states. In this work, we investigate a realistic scenario for the quantum simulation of such systems using cold Rydberg-dressed atoms in optical lattices. We consider spinless fermions on a checkerboard lattice, interacting via the tunable-range effective potential induced by the Rydberg dressing. We perform a detailed analysis of the phase diagram at half- and incommensurate fillings, in the mean-field approximation. We furthermore study the stability of the phases with respect to temperature within the mean-field approximation and with respect to quantum fluctuations using the density matrix renormalization group method. Finally, we propose an implementation protocol, and in particular identify attainable regimes of experimental parameters in which the topological properties of the model become accessible. Our work thereby opens a realistic pathway to the outstanding experimental observation of this predicted phase in state-of-the-art cold atom quantum simulators.

The interplay between many-body interactions and the kinetic energy gives rise to rich phase diagrams hosting, among others, interaction-induced topological phases. These phases are characterized by both a local order parameter and a global topological invariant, and can exhibit exotic ground states such as self-trapped polarons and interaction-induced edge states. In this work, we investigate a realistic scenario for the quantum simulation of such systems using cold Rydberg-dressed atoms in optical lattices. We consider spinless fermions on a checkerboard lattice, interacting via the tunable-range effective potential induced by the Rydberg dressing. We perform a detailed analysis of the phase diagram at half-and incommensurate fillings, in the mean-field approximation. We furthermore study the stability of the phases with respect to temperature within the mean-field approximation and with respect to quantum fluctuations using the density matrix renormalization group method. Finally, we propose an implementation protocol, and in particular identify attainable regimes of experimental parameters in which the topological properties of the model become accessible. Our work thereby opens a realistic pathway to the outstanding experimental observation of this predicted phase in state-of-the-art cold atom quantum simulators.

I. INTRODUCTION
Quantum simulators offer a powerful avenue for the study of many-body physics. These quantum systems mimic the dynamics of complex quantum matter in a highly controllable environment. They are in fact ideal candidates to solve many-body problems whose computational cost on classical computers scales exponentially with the system size. Theoretically proposed in the 80s [1], they are nowadays a reality and can be realized in various physical systems such as photonics, superconducting qubits, and cold ions or neutral atoms [2][3][4]. Here, we focus on cold atomic simulators based on atoms excited to Rydberg states [5,6], which offer rich opportunities for quantum information processing, owing to their long-lived nature and strong long-range interactions leading to the paradigmatic Rydberg blockade effect [7]. Furthermore, individually controlled Rydberg atoms in optical tweezers [8] have emerged as a powerful platform for quantum computation [9][10][11][12][13][14], and for the simulation of quantum spin models [15][16][17][18][19][20], as highlighted by recent observations of 2D spin liquid phases [21,22].
While the strong interactions in Rydberg arrays are typically well captured by spin Hamiltonians, in which kinetic terms accounting for the itinerant nature of the particles can effectively be neglected, one of the challenges in the field is to achieve comparable kinetic and interaction energy scales in order to observe the interplay of interaction and motional effects. This can be achieved in a variety of platforms that feature long-range interactions, such as dipolar quantum gases [23][24][25][26][27][28], or polar molecules [29,30], and in the presence of an optical lattice this allows one to simulate extended Hubbard Hamil-tonians with both non-local interactions and tunneling terms [31][32][33][34][35][36]. Rydberg dressing [37][38][39] has emerged as a powerful alternative in this context. In this approach, instead of exciting the atoms resonantly to a highly excited Rydberg state, in which the energy scale of the strong dipole-dipole interactions dominates over the itinerant dynamics, the atomic gas in the ground state is coupled off-resonantly to the Rydberg state, thereby admixing a reduced amount of Rydberg character to the electronic ground state. Compared to other techniques, Rydberg dressing offers the possibility to tune the strength and shape of interactions, which can be highly adjusted by a proper choice of the atomic and laser parameters of the underlying dressing protocol. Such degree of control has allowed to generate Bell pairs in optical tweezers [40], to engineer long-range [41][42][43] or even distanceselective [44] interactions in Ising Hamiltonians, and to realize extended Fermi-Hubbard Hamiltonians [45] with interaction strengths and kinetic terms of the same order of magnitude. The latter has led to the observation of quench dynamics of a Fermi gas with long-range interactions [45], paving the way for the simulation of other novel phases of quantum matter resulting from the interplay between non-local interactions and the kinetic energy. In particular, this Rydberg dressing toolbox is perfectly suited for the simulation of interaction-induced topological insulators [46], which requires a high control over the ratio of interactions in the presence of a finite tunneling term.
Topological insulators constitute a new paradigm of quantum matter [47,48]: characterized by a global topological invariant, they escape the standard classification of phases of matter and are very robust against local per-turbations such as disorder or interactions. While these phases have been realized in quantum simulators [49][50][51][52][53][54][55], they generally require the engineering of an external gauge field [56,57]. Alternatively, topological insulators can also arise solely from interactions through a symmetry breaking mechanism. In a seminal work [58], it was shown that such an interaction-induced topological insulator, also called topological Mott insulator, can arise for fermions on a hexagonal lattice, with sufficiently strong inter-site interactions. In particular, next-nearest neighbor interactions can give rise to a ground state which breaks the time-reversal symmetry and is characterized by a non-zero topological invariant, the Chern number. Subsequent studies also found topological Mott insulators in other lattice geometries [59][60][61][62][63][64]. Interactioninduced topological phases are quite different from externally induced topological phases [46]. One of the most striking differences is the ground-state degeneracy. In the case of externally induced topological phases, the ground state is non-degenerate, whereas the ground state of the topological Mott insulator is two-fold degenerate, with each of its two sectors being characterized by oppositevalued Chern numbers. These two degenerate ground states with opposite Chern numbers can give rise to interesting effects around half filling such as the appearance of self-trapped polarons or interaction-induced topologically protected edge states, discussed in a previous work by some of us [65].
In this work, we address the timely question of whether the TMI phase can be accessed in quantum simulators based on dressed Rydberg atoms in an optical lattice, under realistic experimental conditions. To this end, we go beyond previous models [62][63][64][65] relying on the simplified assumption of only nearest and next-nearest neighbors interactions, and for the first time properly account for the long-range nature of the Rydberg potential up to fourth order neighbors. Furthermore, we examine the sensitivity of the TMI phase with regard to finite temperature. Our extensive numerical analysis combines mean-field and density-matrix-renormalization group techniques, and is complemented by a thorough discussion of an experimental implementation proposal. Thereby, our study clearly establishes this phase in a robust parameter window, and furthermore provides a clear and experimentally feasible route towards the quantum simulation of the considered topological Mott insulator phase.
The article is organized as follows. In Section II, we review the phase diagram of the model featuring interactions up to next-nearest neighbors, and we introduce the different order parameters characterizing the charge orders and the quantum anomalous Hall (QAH) phase. In Section III, we present a scheme based on dressed Rydberg atoms for the quantum simulation of the model. We review the ingredients required, crucially observing that all such elements have been demonstrated in stateof-the-art setups. We then perform an in-depth study of the phase diagram in Section IV. We discuss the im-pact of longer-range interactions, present in the Rydberg dressing scheme, on the interaction-induced QAH phase. In particular, we show that these can stabilize the QAH phase. We then study how this interaction profile affects the phases at incommensurate fillings around half filling. We additionally probe the stability of the phases at finite temperature. Furthermore, we confirm the stability of the phases with respect to quantum fluctuations with the help of a density matrix renormalization group analysis. Finally, in Section V, we discuss possible parameter regimes, accessible in state-of-the-art experiments, where the QAH phase can realistically be observed.

A. Model
The emergence of a TMI phase has been extensively studied [58][59][60][61][62][63][64][65][66] in lattice systems of spinless fermions described by the extended Fermi-Hubbard Hamiltonian, (1) The first term of the equation describes spinless fermions hopping on a two-dimensional lattice, withĉ † i (ĉ i ) being the fermionic creation (annihilation) operator at lattice site i. The second term represents repulsive interactions, V ij > 0, between fermions on different lattice sites, with local particle number operatorsn i ≡ĉ † iĉ i . In the original proposal, Raghu et al. [58] considered the honeycomb lattice at half filling, for which the noninteracting band structure obtained from the hopping matrix t ij is topologically trivial and exhibits a linear band touching, i.e., Dirac cones. The authors showed that, in the mean-field approximation, the repulsive interactions open a topological gap, leading therefore to an interaction-induced topological phase that they termed Topological Mott Insulator. Subsequent exact diagonalization and DMRG studies of Dirac semimetals, including the semimetallic model of the initial proposal, showed that, beyond the mean-field approximation, interactions favor trivial charge orders with lower energy than the TMI phase [67][68][69][70][71][72][73].
In parallel, perturbative analyses in several models, for which the non-interacting fermionic band ofĤ EFH exhibits instead a quadratic band touching (QBT) also suggested the appearance of a TMI phase [59,[74][75][76]. In the perturbative limit, the TMI phase of such QBT systems was shown to be more stable that in Dirac semimetals, which are more robust with respect to instabilities driven by small symmetry-preserving interactions [59]. More recently, researchers have confirmed, using nonperturbative numerical methods such as DMRG or exact diagonalization, the existence of the TMI phase in many of these QBT systems both for weak and intermediate values of the interactions. This is for example the Hamiltonian comprises a NN interactionṼ1 and hopping amplitude t, and a NNN interactionṼ2 and hopping amplitudes x /t = 0.5 and M = 2, at half-filling and at zero temperature. The Chern number ν is 0 for trivial Mott insulating phases and ν = ±1 in the TMI phase, indicating a ground-state double degeneracy; brighter blue/green corresponds to larger stripe/site-nematic charge order, darker red to larger current loop order. (c-e) Instances of ordered states with spontaneously broken symmetry: (c) site-nematic order, imbalanced density between A and B (two configurations); (d) stripe order, alignment of the particles along S1 or S2, or along their orthogonal axes (four configurations); (e) QAH order: average half-filling and current loops with chirality ij = +1 on square plaquettes of nearest-neighbors, see Eq. (7); for ij = −1, the loops has opposite chirality.
In this work, we focus on this latter case, that is, we consider a checkerboard lattice with a Hamiltonian whereĤ 0 is the non-interacting Hamiltonian [62][63][64][65] and reads [77]: Here, µ is the chemical potential and fixes the particle number in the grand-canonical ensemble at temperature T , t is the nearest-neighbors (NN) hopping amplitude, and J α η is the next-nearest-neighbors (NNN) hopping amplitude, which depends on the sub-lattice α ∈ (A, B) and hopping direction η ∈ (x, y) [see Fig. 1(a)]. The non-interacting band structure exhibits a quadratic band touching for the choice of NNN hopping J A x = J B y = 0.5t and J A y = J B x = −0.5t, corresponding to a π-flux through the unit cell of both sub-lattices. For more general designs, especially in the case of homogeneous and isotropic NNN hopping, the dispersion is linear, as shown in Appendix B. For the interaction, we consider a general Hamiltonian with repulsive interactions of the densitydensity type, which readŝ Here the second sum is performed over the m-th order neighbors ij m of the checkerboard lattice, e.g., ij 1 corresponds to NN terms. The isotropic repulsive interaction between m-th neighbors is then parametrized by the potentialṼ m > 0. As will be discussed in Sec. III, in this work we consider the Hamiltonian in Eq. (4) with interactions up to M = 4, which faithfully describes the repulsive interactions experienced by dressed Rydberg atoms in an optical checkerboard lattice. To provide background, we begin our analysis by first reviewing some known results [62][63][64][65] for M = 2.

B. Half-filling interacting phases
The TMI phase is captured already at the mean-field level. By means of a standard Hartree-Fock decoupling, the repulsive density-density interactions of amplitudeṼ 1 andṼ 2 are approximated aŝ with ξ ij ≡ ĉ † iĉ j andn i ≡ n i , leading to the Hartree-Fock Hamiltonian The Hartree-Fock values ξ ij andn i are found by solving iteratively the resulting self-consistent quadratic Hamiltonian, as described in Appendix A. Figure 1(b) shows the half-filling phase diagram ofĤ HF zero temperature [63,65]. In the limit of vanishing hopping t → 0, the phase diagram hosts two insulating phases which spontaneously break the lattice translational symmetry, as can be seen in Fig. 1(c)-(d). The state resulting from the symmetry breaking is determined by the competition be-tweenṼ 1 andṼ 2 . Consequence of the repulsive densitydensity interaction is an energy cost ofṼ 1 on pairs of particles occupying nearest-neighboring sites, and ofṼ 2 for next-nearest-neighbors. For dominantṼ 1 , low-energy states are characterized by a minimal number of nearestneighboring pairs, conjoined with a maximal density imbalance ρ n ≡n A −n B between the two sub-lattices, thereby giving rise to the so-called site-nematic order. By the same argument, for dominantṼ 2 , the energy penalty of next-nearest-neighbor pairs favors states with stripe density order, characterized by a finite value of the density imbalance ρ s ≡n S1 −n S2 between, e.g., the stripes S 1 and S 2 in Fig. 1(d). As shown in Fig. 1(b), the transition between these two charge-ordered phases happens along the lineṼ 2 =Ṽ 1 /2, when interactions dominate over the tunnelling amplitude. However, when the kinetic energy becomes comparable to the interactions, quantum fluctuations lead to frustration between the two competing charge orders close to the phase transition. In this scenario of charge homogeneity (translational symmetry), the ground state can still be insulating due to the appearance of a current loop order across nearest neighbors which spontaneously breaks time-reversal symmetry [see Fig. 1(e)]. The local order parameter is defined as the staggered sum of currents in a closed loop of nearestneighbors bonds, where ij = +1 if the bond i → j follows the red arrow convention of Fig. 1(e), and ij = −1 otherwise. This phase is known as topological Mott insulator or interaction-induced quantum anomalous Hall (QAH) phase, as each of its two symmetry-breaking ground states with opposite current chiralities is characterized by a global topological invariant, the Chern number [78], Here u 0 k is the lowest single-particle Hartree-Fock band of the Hamiltonian in Eq. (6), which includes the effect of interactions at the mean-field level. The integral is performed over the first Brillouin zone (BZ) of the checkerboard lattice, assuming translational invariance of the two-site unit cell. The Chern number is quantized to integer values in systems with a band gap and is related to the Hall conductivity by 79]. For the topological Mott insulator it assumes one of the two non-trivial values ν = ±1 corresponding to the two sectors of the spontaneous symmetry breaking. For the two other insulating phases present in the phase diagram, the topological invariant takes the value ν = 0, indicating that these phases are topologically trivial. In the next sections, we show that this QAH phase remains present also in the scenario of the realistic long-range interaction potential that describes the interaction between laserdressed Rydberg atom pairs.

III. QUANTUM SIMULATION USING RYDBERG ATOMS
Numerical analyses aiming at unveiling the presence of a topological phase in quantum models are typically carried out in the thermodynamic limit in the meanfield approximation or, when including interactions, using exact or quasi-exact methods but considering systems of limited size. In general, it is computationally hard to study the ground-state properties of interacting two-dimensional systems in the thermodynamic limit and observe its phenomena, such as a spontaneous symmetry breaking and the emergence of a quantized Chern number. When direct observation in, e.g., quantum materials is not practicable, quantum simulation offers an alternative way to reveal theoretically predicted physical properties. In this context, ultracold gases trapped in optical lattices represent a pre-eminent platform for the quantum simulation of interacting Hubbard models [80,81] such as the one given by Eq. (2). The platform enjoys a high level of experimental tunability, allowing for the control on tunnelling and on-site interaction [82]. Furthermore, various detection methods are available for state inspection, from time-of-flight measurements to quantum gas microscopy [83], or magnification [84] techniques.
In this Section we discuss how cold Rydberg gases in suitable lattice geometries represent an ideal platform for the engineering of the interacting Hamiltonians that give rise to the discussed TMI phase. In particular, we will show how the phase can be realized in a checkerboard lattice with a π-flux and with the required long-range interaction terms. Remarkably, the demonstration of all essential elements of the Hamiltonian (2) has been reported in currently available experimental setups for the parameters of our concern.
A. Free Hamiltonian The TMI has been numerically identified in various lattices with a quadratic band touching, including the kagome [60,64] and the checkerboard lattice [62][63][64]. In this work we are interested in the latter case, where the checkerboard is obtained from a square lattice in which a sub-lattice-dependent π-flux on NNN plaquettes is introduced [Eq. (3)]. Regarding the lattice, the design of a wide variety of optical lattice geometries, including the square, has been demonstrated experimentally [85,86] and can be realized by properly adjusting the interference pattern of the standing laser beams. The injection of an artificial flux on NNN plaquettes generates the checkerboard lattice with a quadratic band touching; this can be resolved via band mapping techniques, already used to certify the presence of Dirac points in a free Fermi gas on a tunable honeycomb lattice [86]. The flux insertion has been demonstrated experimentally in cold gases quantum simulators [56,57,87] and requires control over the magnitude, sign and complex phase of the hopping amplitude t. The dynamics of cold gases in lattices, a tight-binding system, occurs via hopping between nearest-neighbors and, marginally, next-nearestneighbors. Coherent control of the hopping amplitude can be attained with several methods. Periodic perturbations of the optical lattice (Floquet techniques) [88] make it possible to reduce, suppress and eventually change the sign of the tunnelling amplitude [89]. Combining a strong lattice tilting with assisted tunnelling allows to exert selective control on hopping terms. The tilting inhibits the tunnelling by introducing inter-site energy barriers larger than the hopping amplitude. Then, the hopping can be activated again in a selective manner and with control on the hopping amplitude, using lattice amplitude modulation [90], or Raman-assisted tunnelling [87]. In particular, Raman-assisted tunnelling allows one to engineer hopping terms with complex amplitudes e i2πφ jkĉ † kĉ j , which can result in finite effective magnetic fluxes on closed paths [56]. The engineering of artificial fluxes was a crucial step for the experimental simulation of static Abelian gauge fields [49,91]. The method can be readily adapted to the π-flux case φ = 1/2 discussed in this paper, which induces the quadratic band touching present in the checkerboard lattice. It is worth stressing that the π-flux does not break explicitly the time-reversal symmetry, as opposed to generic finite fluxes φ = 1/2. As discussed, the symmetry breaking in TMIs occurs by effect of the interactions.

B. Interacting Hamiltonian
Let us next discuss how a density-density inter-particle interaction as in Eq. (4) can be engineered in cold gases experiments. We consider a scheme based on effective interactions between Rydberg-dressed atoms. Rydberg states are electronically excited atomic states with a large principal quantum number [8,[92][93][94]. The laser-driven (single-or two-photon) transition that couples the electronic ground-state or low-lying (meta-)stable state |g and an excited Rydberg state |r [94], in rotating-wave approximation and in the co-rotating frame, is described by the single-particle HamiltonianĤ c = (Ω |r g| + H.c.) + ∆ |r r|, with effective Rabi frequency Ω and detuning ∆. In this work, we consider the repulsive twobody van der Waals interactions experienced between two atoms in the same Rydberg state |r , described by the van der Waals potential U vdW (r) = C 6 /r 6 , where r is the inter-atomic distance and C 6 depends on the Rydberg state [93,94]. The van der Waals interaction between Rydberg atoms is long-ranged and strong at short distances; as an example, U vdW (r 1 ) = 90 MHz for the |28P Rydberg state of 6 Li at r 1 = 752 nm (in this case attractive) [45]. A characteristic effect of the strong Ry- Amplitude of the interactions between Rydbergdressed atom pairs, in a checkerboard lattice. The continuous curve shows the effective interaction potential V (r) from Eq. (9) including only distance-dependent terms to fourth order in Ω/∆ -see Sec. III B. At large inter-particle distances, V (r) decays as β 4 r −6 , consistently with the repulsive van der Waals interaction in the doubly-excited component of the twoatom dressed state (dashed line). Within the blockade radius, r < rc, V (r) converges to an energy plateau. The colored slabs show the ranges of V1 and V2 for which a topological QAH phase emerges, as presented in Figs. 5 and 10. The inset shows an excerpt of the lattice, illustrating the succession of inter-site distances ri: e.g., nearest-neighboring atoms sit on the cyan circle of radius r1. The proximity of r3 and r4 and the comparable magnitude of V3 and V4 require the inclusion of both terms in an analysis beyond V2. Here, r1/rc = 0.67. dberg potential at short distances is the dipole blockade: the laser excitation to the Rydberg state of multiple atoms within a certain exclusion volume is inhibited, as the strong interaction shifts the energy level of a state with multiple Rydberg atoms by more than the line width [93]. The Rydberg blockade has been observed in numerous experiments (see e.g. [9,[95][96][97][98][99][100][101]) and lies at the heart of Rydberg-based analog quantum simulation, e.g., of quantum spin models [8]. Also, the implementation of entangling gates based on the blockade has been demonstrated [9,10], and represents the basis for potential applications in quantum computations [102], under rapid development in recent years [11][12][13][14]103].
The strong interaction within the blockade radius can also be advantageously used in a Rydberg-dressing scheme [61] for the quantum simulation of extended Hubbard-type models. In the limit of far off-resonant laser coupling, i.e., for small values of the parameter β = Ω/∆ 1, the transition between |g and |r are energetically suppressed. The laser induces a weak hybridization of the electronic ground state with the Rydberg state; the strong van der Waals interaction occurring in the marginal Rydberg component of the admixture results in a finite and attenuated spatial-dependent soft-core potential [93]. One can understand the twobody interaction and derive the resulting effective poten-tial surface by looking at a two-atom system.
The ground-state energy of a laser-dressed two-atom system can be obtained as a power series of the perturbation parameter β using, e.g., van Vleck's perturbation theory [61,104]. The spatial-dependent correction to the unperturbed electronic ground state energy up to fourth order in Ω/∆ reads (see Appendix C for the derivation): Here, we have not included the interaction-independent single-particle AC Stark shift 2( We can conveniently fix E p = 2Ω 4 /∆ 3 as an energy scale for the effective interaction and define in the usual way the critical length, or blockade radius, at the full width at half maximum of V (r). Also, we express the discrete inter-site distances in the lattice in units of the lattice spacing r 1 , 2} and i is an index labelling neighbors radii (see Fig. 2). As a result, the inter-site densitydensity effective interaction reads (11) Figure 2 shows the plot of V (r) as well as the discrete V i values. At large distances, the effective potential decays as this shows the suppression of the bare van der Waals interaction by the small prefactor β 4 , that is the probability to find the dressed two-atom system in a doubly excited Rydberg state. At short interatomic distances, the strong enhancement of the van der Waals interaction is accompanied by a vanishing population of the Rydberg-Rydberg component of the admixture, resulting in a soft-core potential, The energy plateau picture does not hold at very short distances, where the overlap of the electronic wave functions becomes more relevant and the van der Waals interaction ceases to correctly describe the interparticle interaction [92]. It should be emphasized at this point that the effective interaction potential V i is a particular case of the interaction Hamiltonian (4), i.e., V i is a constrained parametrization of the more generalṼ m , and depends on a range of controllable independent laser and atomic parameters: the Rabi frequency Ω, the detuning ∆, the lattice spacing a latt , and the van der Waals interaction coefficient C 6 . Below, in Sec. V, we show how these parameters can be adequately tuned to adjust V i and access QAH states in a quantum simulation. Here, we note that fixing V 1 and V 2 , or any other pair of V i , uniquely determines E p and r 1 /r c , and thereby determines the value of the remaining V i . We also remark that the ratio between two interaction amplitudes V i /V j , with i < j, is a monotonically increasing function of r 1 /r c , with a lower bound set by V i /V j → (d i /d j ) 6 in the limit r 1 r c . In particular, V 2 /V 1 has a lower bound of 1/8. The other limit r 1 r c corresponds to an unphysical regime where Eq. (11) is no longer valid, as all long-distance V i would be comparable in magnitude to V 1 .
Since we are interested in studying how the physical properties of model (2) with M = 2 change when we include a finite number of sub-leading long-distance interaction terms, we limit our investigation to a regime in which a truncation of the effective Rydberg potential to V 4 represents a meaningful approximation of the entire effective potential, including the tail. To this aim, we chose to set the condition of V 1 always being at least an order of magnitude larger than the largest discarded interaction term, i.e., V 1 > 10V 5 , corresponding to r 1 /r c ≥ 0.51.
Having comprehensively introduced the model Hamiltonian, we can now proceed to presenting the results of our numerical study of this model.

IV. PHASE DIAGRAM WITH RYDBERG INTERACTIONS
In this Section, we present an extended numerical analysis of the Hamiltonian in Eq. (2) in the presence of long-range interactions beyond next-nearest neighbors, motivated by the long-range character of the effective Rydberg potential, Eq. (11). After showing the effect of adding an arbitraryṼ 3 interaction in the ground-state phase diagram at half filling, we focus on the particular shape of interactions given by the effective Rydberg potential. For the latter, we study the presence of the QAH phase in the phase diagram with the mean-field Hartree-Fock method, and we also discuss the effects of incommensurate fillings on finite-sized systems. Then, in the prospect of a quantum simulation, we examine the robustness of the QAH phase against thermal fluctuations with the finite-temperature Hartree-Fock method. Furthermore, we analyze the stability of the phase beyond the Hartree-Fock ansatz using the DMRG method at zero temperature, which accurately describes the ground states of gapped two-dimensional systems in cylinder geometries with finite widths [105].
A. Hartree-Fock phase diagram

Half filling
To inspect the ground-state phase diagram, we perform a Hartree-Fock study in a large unit cell containing eight sites, illustrated in Figure 3(a), which can host long-range correlators and capture charge orders with a large spatial periodicity.
First, we survey the phase diagram for various, uncon-strainedṼ 1 . As seen in Sec. II B, forṼ 1 (Ṽ 2 ) much larger than any other energy scale in the Hamiltonian, the system is in a gapped site-nematic (stripe) phase. For dom-inantṼ 3 interactions, the density distribution presents two types of charge orders, depicted in Figs. 3(b)-(c). A first observation to make is that these two orders favoured byṼ 3 are incompatible with the density orders generated byṼ 1 andṼ 2 . The consequence of this is an enhanced competition between charge orders asṼ 3 becomes larger. Figure 4 shows the size of the QAH region in theṼ 1 −Ṽ 2 plane for different values ofṼ 3 . One can observe that a fi-niteṼ 3 augments the area of the QAH phase in parameter space. The QAH phase benefits, indeed, from the frustration between competing charge orders:Ṽ 3 supports a different charge order than that of the site-nematic or the stripe phases, ultimately favoring the topologically  Fig. 1: site-nematic in green, stripe in blue and topological QAH phase in red. The QAH phase has a larger order parameter at large V1/t and V2/t. The diagram is bounded below by the largest ratio V1/V2 = 8 attainable within the Rydberg-dressed potential and above by the truncation condition V1 > 10 V5. ordered phase.
Let us now come to the Hamiltonian describing dressed Rydberg atoms, where we emphasize that V m is constrained by Eq. (11) and by our truncation condition V 1 > 10V 5 . Notice that the Hamiltonian includes a finite V 3 term which promotes the stabilization of the QAH phase, as discussed above, but also a finite V 4 term, which favors the sitenematic order generated by V 1 . Notwithstanding, since V 4 is a subleading term, we expect the appearance of the QAH forĤ R also. This is indeed what we observe in the phase diagram ofĤ R , shown in Fig. 5. The sitenematic phase prevails in a large part of the phase diagram, owing to the predominance of V 1 over the other interactions. A QAH phase emerges as V 2 /V 1 increases, and it can approximately be located in the window of V 1 /t ∈ {2, 6} and V 2 /V 1 ∈ {0.9, 0.5}. This latter corresponds to r 1 /r c ∈ {0.51, 0.74}, as can be easily verified using Eq. (11). This interval of r 1 is indicated in Fig. 2 by a cyan slab; the orange slab shows the corresponding range of r 2 /r c . Note that r 2 /r 1 is determined by the lattice geometry; consequently, V 2 does not span the orange slab independently from V 1 . These slabs illustrate, for the checkerboard model and in the presence of dressed van der Waals interactions, where the QAH is to be found on the soft-core potential curve. A firstneighbors distance r 1 too close to the critical distance r c leads the system into a deep site-nematic phase, because all ratios V 1 /V i increase for increasing r 1 /r c ; this deter- mines the right limit of the cyan band, r 1 = 0.74. On the opposite end, the slab is limited by the criterion of V 1 being an order of magnitude larger than V 5 , which we impose to work with a potential truncated to V 4 . The current loop order parameter ξ QAH takes larger values at larger V 1 /t and V 2 /t, as indicated by the darker red color. For V 1 /t ≤ 2 (not shown) we find no presence of either orders, as the system enters a metallic phase. The behavior of ξ QAH in both limits is congruent with what we observed for theṼ 1 −Ṽ 2 model, in Fig. 1, and reaffirms the emergence of the current loop order from an interplay between kinetic energy and interactions. As compared to theṼ 1 −Ṽ 2 model (6), however, we can appreciate a considerably larger QAH region with the Rydberg dressing, as effect of the frustration introduced by the competition between multiple charge orders.

Incommensurate fillings
The interaction-induced QAH phase presents several differences in contrast to a non-interacting fermionic Chern insulator. On the one hand, it exhibits a twofold degeneracy of the ground state at half filling, corresponding to the two sectors of the spontaneous time-reversal symmetry breaking. On the other hand, the rigid band picture around half-filling breaks down due to the presence of correlations, and localized states can appear inside the topological gap. These properties lead to exotic solutions at incommensurate fillings, such as self-trapped polarons or domain walls interpolating between the two sectors of the spontaneous symmetry breaking [65].
We find these solutions also in the presence of the effective Rydberg potential, as shown in Fig. 6, with the unrestricted Hartree-Fock method described in Appendix A. The quantity δ counts the number of particles added to the half-filled state. In the case δ = 1, Fig. 6(a) shows that the added particle does not populate the conducting band but instead occupies a midgap localized state induced by interactions, a self-trapped polaron. In this solution, the local current loop order ξ QAH changes its sign inside the polaron region, which can be understood as a collapsed domain wall. As δ increases, the number of mid-gap states and the polaron size increases [see Fig. 6(b)]. Eventually, we observe the formation of a ringshaped domain wall separating an inner and outer region with opposite current chiralities [see Fig. 6(c)] which correspond to opposite Chern numbers inside and outside the ring.

Finite temperature analysis
We have seen above that, at zero temperature, the QAH phase appears within the Rydberg potential for a wide range of interactions. Let us now study the stability of the phase with respect to temperature by means of the finite-temperature Hartree-Fock method (see Appendix A). Here, the occupations of the Hartree-Fock single-particle states with energies E i are given by the Fermi-Dirac distribution, Moreover, we study the typical temperatures needed in order to resolve the spatial structures around half filling, shown in Fig. 6.
FIG. 7. Gap in the energy band structure at finite temperatures, for effective Rydberg interactions with V2/t = 3, corresponding to the yellow cut in the inset (clip of the phase diagram in Fig. 5). At V1/t 5 we observe the transition between two insulating phases with a finite gap, from QAH to site-nematic. No evident effects to the zero-temperature gap Egap(T = 0) (not shown) are observed for kBT /t ≤ 0.2. At higher temperatures, the gap begins to close, affecting first and mostly the QAH phase. Homogeneous phase at half filling. Figure 7 presents the energy band gap E gap along a cut in the phase diagram, indicated by the yellow line in the inset. The gap refers to the Hartree-Fock single-particle band structure, at half-filling and for different temperatures. While the traditional notion of topology is typically defined at zero temperature, one can still use the notion of a topological invariant at finite temperature for the density matrix [106][107][108][109][110], provided that the thermal energy scale k B T is lower than this insulating gap E gap .
As a first remark, we note that the zero-temperature gap E gap (T = 0), not shown, is indistinguishable from the gap at k B T /t = 0.2. The sharp discontinuity at V 1 5t pins the phase transition between the QAH and the site-nematic phase. This jump in the value of the gap, in agreement with the first-order nature of the transition, can be understood from the fact that this quantity is correlated with the value of the order parameter of the respective phase: when approaching the transition from the QAH side, both the current loop order as well as the gap are enhanced, whereas when approaching it from the site-nematic phase both the charge order and the gap vanish. With regards to the QAH phase, the gap E gap is of the order of the hopping rate t in most of the QAH region, taking the maximum value of about 4t around V 1 = 5t, V 2 = 3t. From the zero temperature gap analysis, one would estimate that the topological QAH phase is robust for temperatures up to a few t/k B . However, given the interacting nature of the Hartree-Fock band structure, a finite temperature calculation of the gap is required in order to establish the critical temperature of the QAH phase. As shown in Fig. 7, the gap decreases non-linearly with increasing temperature, affecting most rapidly states with smaller gap at T = 0. Ultimately, we can roughly estimate a critical temperature for the appearance of the QAH phase of about T c = t/k B , well below the temperature T 4t/k B suggested by the gap structure at T = 0. Using t = 1.7 kHz, a value on the scale of current experimental realizations [45], we obtain a critical temperature of T c = 82 nK around V 1 = 5t and V 2 = 3t.
It is also worth to visualize the effect of the temperature without restrictions to a point or a line of the phase diagram. To this end, in Figure 8, we show the whole phase diagram ofĤ R at three different finite temperatures. Up to a temperature of k B T /t = 0.2 no appreciable alteration to the phase diagram is observed. The QAH phase emerges from the competition between kinetic energy and the frustrated charge order driven by the interaction. As such, it results to be most fragile against thermal fluctuations. As the temperature increases, the valence-conduction gap of all insulating phases progressively reduces, with a major impact on the QAH phase. The QAH gap closes first in the region of lower V 1 and V 2 , where the zero-temperature current loop is smaller, leading to a gapless semi-metallic phase.
Defects around half filling. Along with a closing gap, at rising temperatures, the mid-gap states progressively disappear, as they mix with the lower bulk band. As an example of such behavior, we study the effect of finite temperature for the case δ = 5, which at zero temperature corresponds to a ring-shaped domain wall [ Fig. 9(a)]. In Fig. 9(b) we show the results for a finite temperature k B T /t = 0.05; while there is no appreciable difference of the order parameter in real space compared to the zero temperature case, the conducting band starts to be populated. For an even higher temperature k B T /t = 0.2 the spatially homogeneous QAH phase without mid-gap states is recovered, as seen in Fig. 9(c).  However, notice that in this homogeneous solution the excess particles are distributed in the upper band, destroying the gap insulating nature of the phase.

B. DMRG phase diagram
In order to corroborate the stability of the QAH phase beyond the mean-field approach used in the previous section, in the case of Rydberg interactions, we perform a density-matrix-renormalization group (DMRG) study in its matrix-product-state (MPS) formulation [105,111,112]. We consider a cylindrical geometry of the checkerboard lattice with an infinite size along the longitudinal direction (iDMRG). Due to the one-dimensional nature of the DMRG algorithm, the cylinder is mapped to a onedimensional chain in a snake-like folding along the radial direction, at the cost of introducing effective long-range couplings. The latter limits us to cylinder widths up to L y = 6 unit cells (12 physical sites). By using a maximum bond dimension χ max = 3000, we get truncation errors of the infinite MPS of the order 10 −5 at most.
As shown in Fig. 10(a), the DMRG calculation confirms the site-nematic to QAH phase transition when varying V 2 /t for a fixed V 1 /t = 4, which determines the value of V 3 /t and V 4 /t according to the dressed Rydberg potential of Eq. (11). One observes a shift of the QAH boundary compared to the Hartree-Fock phase diagram of Fig. 5, which is expected since mean-field methods are known to be less accurate in the vicinity of a quantum phase transition. Specifically, at V 1 /t = 4 with DMRG the QAH phase appears at V 2 /t 1.5 and disappears for V 2 /t > 2.5. In contrast, with a Hartree-Fock ansatz the QAH phase is not present until V 2 /t > 2 and disappears for V 2 /t > 3.5. The shift points in a direction advantageous from an experimental point of view; the realization of coherent systems with strong long-range interactions is difficult, therefore the possible appearance of the QAH phase at lower V 2 , already smaller than V 1 , is a favorable sign. We also establish the existence of the QAH phase for a wide range of V 1 /t ratios, as shown in Fig. 10(b). As hinted by the previous Hartree-Fock calculations, the Rydberg potential stabilizes the phase in a larger window of interaction strengths compared to the simplifiedṼ 1 -Ṽ 2 model [63].

V. EXPERIMENTAL PARAMETERS ANALYSIS
In this section we study the relevant experimental parameter regimes to simulate the model of Eq. (14) with Rydberg dressed atoms [5,[37][38][39] in an optical lattice. Notice that this dressing technique has been already widely used in several recent experiments [40][41][42][43][44][45], which have succeeded in engineering tunable long-range interactions in two-dimensional systems. Of particular interest for the simulation of the TMI phase is the experiment of Ref. [45], which has allowed for the observation of a long coherence time in a two-dimensional Fermi lattice gas, in the presence of tunnelling and inducing strong non-local interactions [45]. Such system paves the way towards the quantum simulation of other Fermi-Hubbard Hamiltonians with long-range interactions, including models with topological properties.
Motivated by this prospect, we investigated and identified parametric regions for which a similar experimental system would well approximate the model of (14) and allow one to reach the interaction-induced topological phase. To begin the analysis, we consider the coherence time as a crucial figure of merit in the context of manybody quantum simulations. In a Rydberg-dressed cold gas, coherence is affected by spontaneous decay of atoms from the Rydberg state, which limits the time scale up to which the effective Hamiltonian (14) faithfully describes the system; in a simple single-particle picture, the time scale is given by the effective Rydberg decay rate where Γ 0 is the bare decay rate of the Rydberg state. Thus, Γ eff sets a lower limit for both the hopping and the interaction rates, namely that t, V i > Γ eff . We now show that the experimental conditions allowing one to engineer and observe the topological phase can be attained by a suitable choice of atomic and laser system parameters. The choice of the atomic species and the target Rydberg state determines the van der Waals interaction strength C 6 , the bare decay rate Γ 0 , the lattice spacing a latt at which the species can be trapped, and the order of magnitude of the tunneling amplitude t [113]. The remaining free parameters are the detuning ∆ and the Rabi frequency Ω of the dressing laser fields.
We envision the quantum simulation being based on a light atomic species in order to favour the itinerant character of the TMI phase. An example would be the fermionic isotope 6 Li of Lithium. This species can be trapped in optical lattices with a lattice spacing a latt = 752 nm [45], and is accompanied by a fast tunneling which we assume to be t 3 kHz in the following. A repulsive and isotropic van der Waals interaction can be found between atoms in |nS Rydberg states, for which the series of C 6 values can be calculated using, e.g., the library from Ref. [114]. Here we consider a Rydberg coefficient C 6 100 MHz a 6 latt , which can be attained for principal quantum numbers roughly above n = 33. For these states, the radiative lifetime is estimated to be approximately 1/Γ 0 30 µs [115]. We fix the parameters mentioned so far and leave out Ω and ∆ as independent variables. Figure 11(a) shows the figure of merit t/Γ eff , i.e., the ratio between the tunneling rate and the effective Rydberg decay rate, as a function of the free parameters Ω and ∆, as given by Eq. (16). An important remark here is that, for each point of this Fig. 11(a), the values of the interactions V i are determined by the rest of the fixed experimental parameters discussed above. In particular, we only show t/Γ eff in the region comprised between the equipotential lines V 1 = t and V 1 = 10t given by the condition Notice that, for this range of interaction strengths, relevant for the simulation of the TMI phase, the hopping rate is always larger than the decay rate, i.e., t/Γ eff > 1. We now discuss a specific region of ∆ and Ω that allows one to reach the TMI phase. In particular, we consider the yellow line at a fixed ∆ = 6 MHz in Fig. 11(a), which corresponds to fixing the ratio V 1 /V 2 = 0.62, as can be inferred from Eqs. (10)- (11). The interaction parameters spanned by this yellow line are mapped into the phase diagram of Rydberg dressed atoms in Fig. 11(b). The current loop order parameter ξ QAH along this cut is also shown in Fig. 11(c): within this parametric span one can access the topological phase.
As an example, we can choose to pin V 1 = 4t, which corresponds to V 2 = 2.5t, V 3 = 0.63t, V 4 = 0.34t; this specific set of V i is illustrated by the circles shown in Fig. 2. For this example, we obtain a large current loop order ξ QAH 0.08, a moderate decay rate leading to t/Γ eff 4, and the ground state is relatively robust to thermal fluctuations with a critical temperature of k B T 50 nK, obtained from the finite-temperature gap analysis of Fig. 7.
As a conclusive note, we remark that the observation of collective phenomena in the quantum simulation set-up considered here hinges on the validity of Eq. (16) to describe the effective decay rate for Rydberg-dressed atoms. Previous experiments have observed a rather large decay rate, scaling with the number of particles. This phenomenon is modelled as a blackbody-driven collective resonant decay from the dressed state, an avalanche mechanism triggered by the first individual atom decay that drives the entire system out of its simulation task. The suppression of this avalanche mechanism will be crucial in order to scale the system size and was to a good extent achieved by the authors of [45], thereby opening the door to quantum simulating the extended Fermi-Hubbard model and topological Mott insulating phase discussed in our present work.

VI. CONCLUSIONS AND OUTLOOK
In this work, we investigated a topological Mott insulating model and addressed the question regarding its quantum simulation with state-of-the art experimental methods. Precisely, we considered a Fermi-Hubbard Hamiltonian on a checkerboard lattice with inter-site density-density interactions, a model which is known in the literature to host a quantum anomalous Hall phase in the V 1 − V 2 case, i.e., when featuring interactions up to next-nearest neighbors [62][63][64][65]. Here, we studied the impact of longer-range interactions on the V 1 − V 2 model and observed that the interaction between thirdneighbors stabilize and enlarge the QAH region in the ground-state phase diagram. This result eases the requirements for an experimental realization, given that realistic long-range interaction potential profiles generally comprise a non-vanishing coupling beyond secondneighbors. We focused on the Hamiltonian modelling a lattice gas of Rydberg-dressed Fermi atoms, exhibiting a long-range two-body effective interaction potential, and we studied the physical properties of this model. The choice of the effective Rydberg potential is motivated by the technological progress reached in the field of optically trapped cold Rydberg atoms [8], as highlighted by recent experiments [45].
Interestingly, we observed that the Rydberg-dressed model Hamiltonian hosts a larger QAH phase in the ground-state phase diagram, as compared to the V 1 − V 2 model. This constitutes an encouraging result for the purpose of quantum simulations, since realistic longrange interactions generically present finite beyondsecond-neighbors coupling. We substantiated our analysis addressing real-system effects that can arise under ordinary conditions in a laboratory: we studied the fate of self-trapped polarons and interaction-induced domain walls at incommensurate fillings, we analyzed the stability of the QAH phase with respect to temperature with a mean-field approach and with respect to quantum fluctuations with DMRG. Furthermore, we discussed realistic ranges of the experimental parameters in a cold gases setup, which allow to access a TMI state of matter.
This work provides a clear route towards the experimental realization of an interaction-induced topological phase in cold-atom quantum simulators. In this context, it is important to better understand how to bring the quantum state into the interaction-induced topological phase. The task requires finding paths in the phase diagram to cross from, e.g., an initial state in the charge ordered phase, which can be prepared in experiments, to the QAH phase via a continuous phase transition. One route could be an adiabatic state preparation through a ramping protocol [116,117]. Another interesting question concerns the detection of the QAH phase. For this, detection schemes developed for non-interacting topological phases in cold atom quantum simulators exist [50,[118][119][120][121][122][123] which require to be generalized and adapted to interacting systems. Finally, as the system also exhibits topological defects, it would be interesting to study the dynamics of the formation of defects [124][125][126] when changing the speed of a ramping protocol for state preparation, and by means of that characterize the topological nature of the interaction-induced QAH phase.
Code availability. The codes for the unrestricted and restricted Hartree-Fock studies of this work are openly available in the repository [127].
which one constrains the Hartree-Fock values to be periodic within a certain unit cell, resulting in a numerical simplification of the algorithm.
Starting from an initial guess for the value of these Hartree-Fock parameters, the Hamiltonian matrix h ij can be diagonalized by means of a Bogoliubov transformation U defined by (f †  1 ,f 1 , . . . ,f † d ,f d ) T = U (ĉ † 1 ,ĉ 1 , . . . ,ĉ † d ,ĉ d ) T , such that the Hamiltonian takes a diagonal form:Ĥ Notice that, in the mean-field approximation, it is feasible to consider large system sizes d because the numerical complexity scales polynomially in d (diagonalization of a d × d matrix h ij ) in contrast with the exponential scaling of the general many-body case. Finally, in the equilibrium state at temperature T , the occupation number of the Bogoliubov modesf i is given by the Fermi distribution where µ is the chemical potential that is used to fix the total particle number trough the condition At half filling and T = 0 one gets f † if i = 1 for the N/2 lower energy states, and 0 for the other ones.
Therefore, any expectation value in the original fermionic basis can be computed using Eq. (A4) together with the Bogoliubov transformation U . In particular, one can compute the new values of the Hartree-Fock parametersn i and ξ ij , which give a new Hamiltonian matrix h ij that can be diagonalized again following the procedure described above. This process is iterated until convergence of the Hartree-Fock parameters is achieved. In order to avoid metastable solutions, one needs to compare the free energies F HF of different converged solutions. The free energy is defined as which, at T = 0, is simply given by the expectation value ofĤ HF . Importantly, in order to converge to solutions that break Hamiltonian symmetries, it is crucial that the initial Hartree-Fock parameters already break them. For instance, the QAH phase requires initial complex values of ξ ij to break time-reversal symmetry, that is ξ ij should be initialized to numbers with a finite imaginary part, and charge ordered phases require spatially inhomogeneous distributions ofn i .

Restricted Hartree-Fock with eight-site unit cell
The checkerboard lattice is a bipartite two-dimensional Bravais lattice that can be uniquely defined with a twosite unit-cell coordinate, accounting for sub-lattices A and B, and two unit vectors. The Fourier-transformed free HamiltonianĤ 0 of Eq. (3) has thus a twodimensional matrix form in k-space, accounting for hopping and interactions between the two sites A and B of the unit cell.
The inclusion of density-density interactions inĤ int of Eq. (4) breaks the block diagonal structure ofĤ 0 in k-space, as interactions represent scattering processes which couple modes with different ks. However, when working with the Hartree-Fock HamiltonianĤ HF described in the previous section, one can artificially impose a certain spatial periodicity of the Hartree-Fock parameters, and recover a block diagonal structure ofĤ HF in k-space.
This is known as the restricted Hartree-Fock method, and is typically used at particle fillings commensurate with the lattice size, where one expects that interactions preserve a certain translational symmetry. In this method, it is important to do a proper choice of the unit cell size. A too small cell size may lead to constrictions: in k-space, two-body energy terms acting on a distance larger than the extent of the cell become identical to existing shorter-range terms and renormalize them. Also, charge orders with a periodicity on larger length scales can not be captured. In particular, in order to resolve the charge density distribution associated to dominantṼ 3 interactions, shown in Fig. 3 of the main text, a four-site square cell is not sufficient. Therefore, in this work, we consider the eight-site cell depicted in Fig. 12. In what follows, we use this unit cell to expressĤ HF in k-space and also to derive the Hartree-Fock parameters.
The Fourier transform of the real-space annihilation operator is defined aŝ where N uc is the total number of unit cells, in this case corresponding to the number of sites, and r mn is a linear combination of the Bravais vectors [e 1 , e 2 ] spanning over the lattice sites: In model (14), we include inter-site interactions up to fourth-neighbors, each one yielding a hopping term of the same neighboring order in Hartree-Fock approximation. The eight-sites lattice cell in Figure 12 accommodates all hopping terms distinctly. The geometry of the cell correctly allows for a complete tiling of the lattice via [a 1 , a 2 ], the two unit-cell Bravais vectors: In the unit-cell framework, one can map the spatial coordinates r mn onto a combined coordinate of the cell index mn and an additional index α that identify the cell sites, as indicated by the enumeration in Fig. 12: r mn → R αmn = ma 1 + na 2 + r α , with m, n ∈ Z.
(A9) Ultimately, we can rewrite the Fourier transform as: a. Bare hopping t, first-neighbors Let us now expressly illustrate, as an example, the mathematical steps yielding the nearest-neighbors hopping in k -space; the procedure is identical for other elements in the Hamiltonian. Expanding the kinetic energy termĤ t = −t <i,j> ĉ † iĉ j + H.c. , we have: N uc p,q∈FBZ <α,β> mn ĉ † α (p)ĉ β (q)e i(q−p)·R1mn e −iq·r αβ + H.c. , where <, > symbolizes the restriction over all firstneighboring pairs and r βα = R βmn − R αmn is the associated vector, also illustrated by the green segments in Fig. 12; here, R 1mn points at cell mn, whose origin coordinate we arbitrarily fixed to the bottom-left site α = 1. The sum over mn gives a Dirac delta function in q − p, which reduces the k -dependency to q only: H t = −t q∈FBZ <α,β> ĉ † α (q)ĉ β (q)e −iq·r αβ + H.c. .
(A12) Ultimately, dropping the explicit q-dependency of the operators to ease readability and following the labelling in Fig. 12 In real-space coordinates, the Hartree-Fock parameters ξ ij could be in principle different from one another. Assuming translational invariance over the unit cells, the number of first-neighboring pairs reduces to the sixteen elements shown in Figure 12. The expectation values are calculated as follows: c. First-neighbors interaction Lastly, using Eq. (A1), we can derive the k -space form of the first-neighbors density-density interaction, in the mean-field approximation. The operator in real-space coordinate reads: (A15) The mean-field approximation introduces constant and diagonal on-site terms and renormalizes the bare hopping. After Fourier transformation, we obtain: As discussed in Section II, numerical evidences [58][59][60][61][62][63][64][65][66][67][68][69][70][71][72][73][74][75][76] suggest that a topological Mott insulating phase emerges in the presence of a quadratic band touching. The free Fermi-Hubbard model on a checkerboard lattice with only nearest-neighbors hopping does not possess this property in the band structure. In fact, conduction and valence bands touch but present a linear dispersion; this can be changed by introducing a bipartite second-neighbors hopping on the two sub-lattices A and B, J A x /t = J B y /t = +0.5 and J A y /t = J B x /t = −0.5, as shown in Fig. 1 in the main text. This specific design of the hopping introduces a so-called π-flux on square second-neighbors plaquettes of both sub-lattices, which gives the quadratic band touching. We can see this analytically by studying the model HamiltonianĤ 0 [Eq. (3)] on a two-site unit cell. In this case, one can show that the dispersion relation reads E kx,ky = −[S x cos(2k x ) + S y cos(2k y )] (B1) ± [D x cos(2k x ) + D y cos(2k y )] 2 + 16t 2 cos(k x ) 2 cos(k y ) 2 where S x = J A x + J B x , D x = J A x − J B x and likewise for S y , D y , and where k x,y = k · e 1,2 . For J A µ = J B µ = 0, S µ = D µ = 0 [see Fig. 13(a)] the system reverts to the bare Fermi-Hubbard gas with a linear dispersion, around the points where conduction and valence bands touch: E kx,ky = ±4t cos(k x ) cos(k y ). (B2) The particular case J A x = J B y = 0.5t = −J A y = −J B x considered in the article, gives S µ = 0, D µ = ±t: E kx,ky = (B3) = ±t [cos(2k x ) − cos(2k y )] 2 + 16 cos(k x ) 2 cos(k y ) 2 .
In this case, the bands touch at a zero energy point in k x , k y = ±π/2. The quadratic behaviour, shown in Fig. 13(b), is easily derived, for instance fixing k y = π/2: