Correlation-induced insulating topological phases at charge neutrality in twisted bilayer graphene

Twisted bilayer graphene (TBG) provides a unique framework to elucidate the interplay between strong correlations and topological phenomena in two-dimensional systems. The existence of multiple electronic degrees of freedom -- charge, spin, and valley -- gives rise to a plethora of possible ordered states and instabilities. Identifying which of them are realized in the regime of strong correlations is fundamental to shed light on the nature of the superconducting and correlated insulating states observed in the TBG experiments. Here, we use unbiased, sign-problem-free quantum Monte Carlo simulations to solve an effective interacting lattice model for TBG at charge neutrality. Besides the usual cluster Hubbard-like repulsion, this model also contains an assisted hopping interaction that emerges due to the non-trivial topological properties of TBG. Such a non-local interaction fundamentally alters the phase diagram at charge neutrality, gapping the Dirac cones even for small values of the interaction. As the interaction strength increases, a sequence of different correlated insulating phases emerge, including a quantum valley Hall state with topological edge states, an intervalley-coherent insulator, and a valence bond solid. The charge-neutrality correlated insulating phases discovered here provide the sought-after reference states needed for a comprehensive understanding of the insulating states at integer fillings and the proximate superconducting states of TBG.

Twisted bilayer graphene (TBG) provides a unique framework to elucidate the interplay between strong correlations and topological phenomena in two-dimensional systems. The existence of multiple electronic degrees of freedom -charge, spin, and valley -gives rise to a plethora of possible ordered states and instabilities. Identifying which of them are realized in the regime of strong correlations is fundamental to shed light on the nature of the superconducting and correlated insulating states observed in the TBG experiments. Here, we use unbiased, sign-problem-free quantum Monte Carlo simulations to solve an effective interacting lattice model for TBG at charge neutrality. Besides the usual cluster Hubbard-like repulsion, this model also contains an assisted hopping interaction that emerges due to the non-trivial topological properties of TBG. Such a non-local interaction fundamentally alters the phase diagram at charge neutrality, gapping the Dirac cones even for small values of the interaction. As the interaction strength increases, a sequence of different correlated insulating phases emerge, including a quantum valley Hall state with topological edge states, an intervalley-coherent insulator, and a valence bond solid. The charge-neutrality correlated insulating phases discovered here provide the sought-after reference states needed for a comprehensive understanding of the insulating states at integer fillings and the proximate superconducting states of TBG.
The recent discovery of correlated insulating and superconducting phases in twisted bilayer graphene (TBG) [1][2][3] and other moiré systems [4][5][6][7] sparked a flurry of activity to elucidate and predict the electronic quantum phases realized in their phase diagrams . Because the low-energy bands of TBG have a very small bandwidth, of about 10meV at the magic twist angle, the Coulomb interaction, which is of the order of 25meV, is expected to play a fundamental role in shaping the phase diagram [1,13,20,22]. Indeed, insulating states have been reported at all commensurate fillings of the moiré superlattice [10], signaling to the importance of strong correlations. Besides correlations, topological phenomena have also been reported, including a quantum anomalous Hall (QAH) phase [52,53].
An important issue is the nature of the quantum ground state at charge neutrality, characterized in real space by 4 electrons per moiré unit cell, and in momentum space by Dirac points at the Fermi level. Experimentally, a large charge gap characteristic of an insulating state was reported in transport measurements in Ref. [10] and in STM measurements in Ref. [11], despite no obvious alignment with the underlying hBN layer. The fact that this gap is not observed in all devices has been attributed to inhomogeneity [10]. Theoretically, because the electronic states in TBG have several degrees of freedom -spin, valley, and sublattice -various possible ground states can emerge. Indeed, Hartree-Fock calculations of the continuum model at charge neutrality found various possible phases, such as orbital-magnetization density-waves, valley polarized states, and states that spontaneously break the three-fold rotational symmetry of the moiré lattice [39,43,44,54,55]. To distinguish among these different possibilities, it is desirable to employ a method that is not only unbiased, but that can also handle strong correlations.
Large-scale quantum Monte Carlo (QMC) simulations provide an optimal tool, limited only by the finite lattice sizes. Although such a limitation makes it impossible to simulate a model with thousands of carbon atoms per moiré unit cell, it is very well suited to solve lattice models on the moiré length scale. At charge neutrality, the non-interacting part of the model has only Dirac points at the Fermi level. The crucial part of the model, however, is the interacting part, which governs the system's behavior in the strong coupling regime. At first sight, based on the analogy with other strongly-correlated models, it would seem enough to consider a cluster Hubbard-like repulsion as the main interaction of the problem. Previously, some of us used QMC to simulate this model, which does not suffer from the infamous fermionic sign-problem [56,57]. The result was a variety of valence-bond insulating states, which however only onset at relatively large values of the interaction U, of the order of several times the bandwidth W. Below these large values, the system remained in the Dirac semi-metal phase.
However, microscopically, the full interaction of the lattice model can be derived from projecting the screened Coulomb arXiv:2004.12536v2 [cond-mat.str-el] 4 May 2020 repulsion on the Wannier states of TBG. The latter turn out to be quite different than in other correlated materials, as they have nodes on the sites of the moiré honeycomb superlattice and a three-peak structure that overlaps with Wannier functions centered at other sites [20][21][22]. Recent work has shown that this leads to the emergence of an additional and sizable non-local interaction, of the form of an assisted-hopping term [32,48]. This new interaction ultimately arises from the fact that, in a lattice model, the symmetries of the continuum model cannot all be implemented locally, a phenomenon dubbed Wannier obstruction [22]. Therefore, the assistedhopping interaction is not a simple perturbation, but a direct and unavoidable manifestation of the non-trivial topological properties of TBG.
In this paper, we study the impact of the assisted-hopping interaction on the ground state of TBG at charge neutrality via sign-problem-free QMC simulations. We find that such a term qualitatively changes the phase diagram, as compared to the case where only the Hubbard interaction is included. In particular, the Dirac semi-metal phase is no longer stable, but is gapped already at weak-coupling. We show that this gap is a manifestation of a quantum valley Hall (QVH) state, characterized by topological edge states. We confirm this weak-coupling result by unrestricted Hartree-Fock (HF) calculations of the same model simulated by QMC. The HF calculations, well suited for weak interactions, also show that the QVH state is a robust property of the weak-coupling regime, and is directly connected to the assisted-hopping term. As the interaction strength increases, a different type of insulating phase arises, displaying intra-valley coherence (IVC) order. This onsite IVC order breaks the spin-valley SU(4) symmetry of the interacting part of the model, resembling recently proposed ferromagnetic-like SU(4) states proposed to emerge in TBG at charge neutrality and other integer fillings [32,43]. Upon further increasing the interaction, a columnar valence bond solid (cVBS) insulator state appears, favored by the Hubbard-like interaction [56,57]. Importantly, the presence of the assisted-hopping term makes the QVH and IVC states accessible already for substantially smaller values of U/W, as compared to the case where there is only Hubbard repulsion. Therefore, the experimental observation of such quantum states in TBG at charge neutrality would provide strong evidence for the importance of non-local, topologically-driven interactions in this system.

Model, symmetry analysis and method
As shown in Fig. 1 (a), our model describes two valleys (orbitals) of spinful fermions on the honeycomb lattice that is dual to the triangular moiré superlattice. The Hamiltonian is given by H = H 0 + H , with the non-interacting tight-binding term: where c † ilσ (c ilσ ) denotes creation (annihilation) operators of electrons at site i, valley l = 1, 2 with spin σ =↑, ↓. The nearest neighbor hopping t = 1 and we use the bare bandwidth W = 6t as the energy unit in the remainder of the paper.
While this simple band dispersion displays Dirac points at charge neutrality, it does not faithfully reflect the detailed band dispersion obtained by DFT calculations and the topological features of the bands. However, the precise form of H 0 is not expected to play a dominant role in the strong-coupling regime [32], which is our focus here. For this same reason, it is fundamental to correctly capture the effective interaction of this lattice model. Here, we consider an interaction term H that contains two contributions: where U sets the overall strength of the Coulomb interaction.
The two contributions in Eq. (2) consist of the cluster charge Q ≡ j ∈ n j 3 , with n j = lσ c † jlσ c jlσ , and the cluster assisted hopping T ≡ j,σ ic † j+1,1σ c j,1σ − ic † j+1,2σ c j,2σ + h.c. . Here, the index j = 1, . . . , 6 sums over all six sites of the elemental hexagon in the honeycomb lattice. The pre-factor α controls the relative strength of the two interactions. Hereafter, we fix the electronic filling strictly at the charge neutrality point, where there are four electrons per hexagon once averaging over the lattice.
This form of H , introduced in Ref. [32], follows from projecting the screened Coulomb interaction on the Wannier states of TBG. The cluster charge term Q is analogous to the Hubbard onsite repulsion in the standard Hubbard model; the reason why it extends over the entire hexagon is because of the screening length set by the separation between the gates in a TBG device and because of the overlap between Wannier states of neighboring sites. In TBG, the Wannier wave-functions are not peaked at the honeycomb sites, but instead are extended and peaked at the centers of the three neighboring hexagons. Therefore, one single Wannier state overlaps spatially with other Wannier states on neighboring sites, leading to the cluster charging term Q . On the other hand, the origin of the assisted hopping term T is topological. As explained in Refs. [22,58], there is an obstruction to construct fully symmetric Wannier states for the isolated nearly-flat bands if one attempts to extend the symmetries of the continuum model of TBG to a lattice model. However, it is still possible to construct Wannier states by implementing the valley-related symmetry C 2 T non-locally [20,21]. Here, C 2 refers to two-fold rotations with respect to the z-axis and T , to time-reversal. As a consequence, the spatial integral of two neighboring Wannier states inside a single hexagon becomes nonzero, giving rise to the crucial T term. The coefficient α is of order unity, and depends on the details of the Wannier states [32].
An interesting feature of H is its emergent SU(4) symmetry. To illustrate this, we introduce the spinor ψ i = c i1↑ , c i1↓ , c i2↑ , c i2↓ T and rewrite the interactions as:  (2), i.e. the valley U(1) symmetry and the two independent spin SU(2) rotations for the two valleys [22,38]. Thus, strictly speaking, the SU(4) symmetry is exact for H but only approximate for H 0 . To solve the model H = H 0 + H non-perturbatively, we employ large-scale projection QMC simulations [56,57]. This QMC approach, employed in several previous studies [56,57,[59][60][61][62], provides results about the T = 0 ground state, the correlation functions (which are used to determine broken symmetries), and the electronic spectra (both single-particle and collective excitations). Despite the presence of the assisted hopping interaction, the model at charge-neutrality does not suffer from the sign-problem (see supplemental material (SM) [63]). Thus, it can be efficiently simulated by introducing an extended auxiliary bosonic field that dynamically couples to the electrons on a hexagon -in contrast to the standard Hubbard model, where the auxiliary field is local. Details about the QMC implementation, as well as comparison with results from exact diagonalization, are discussed in the SM [63].
We also complemented the unbiased QMC simulations with self-consistent HF calculations, which are well-suited for the weak-coupling regime, and can be employed even when additional terms are included in H that introduce a sign-problem for QMC. The HF approach is fully unrestricted in the sense that H 0 + H is mean-field decoupled in all channels, and free to acquire any value in site-, spin-, and valley-space . Further technical details, including the resulting coupled set of (real-space) self-consistency equations, can be found in the SM [63]. In the regime of weak interactions, we find excellent agreement between the results obtained from HF and QMC, as shown below.

Quantum valley Hall phase, intervalley-coherent insulator, and valence-bond solid
The QMC-derived phase diagram for the ground states at charge neutrality is shown in Fig. 1 (b) as a function of U/W and α. We emphasize that while U gives the overall magnitude of the total interaction term, α is proportional to the relative strength between the assisted-hopping and cluster-charge terms. We find that three types of correlated insulating phases emerge in the phase diagram: the quantum valley Hall (QVH) phase, the intervalley-coherent (IVC) phase, and the columnar valence bond solid (cVBS). The QVH phase is the ground state for small U values and is characterized by a gap in the single-particle spectrum. This gap can be extracted from the imaginary-time decay of the Green's function along a high-symmetry path of the Brillouin zone (BZ), G(k, τ) ∼ e −∆ sp (k)τ . Fig. 2(a) shows the enhancement of the single-particle gap at the K point of the BZ as a function of U for a fixed α = 0.45 (blue points). Together with Fig. 2(b), one sees the gap opens at the entire BZ at infinitesimally small U. In many honeycomb lattice models, the Dirac cone at the K point is protected by a symmetry, and the semi-metal phase is robust against weak interactions [56,57,59,64]. In TBG, however, the relevant symmetry, C 2 T , cannot be implemented locally due to the topological Wannier obstruction. This opens up the possibility of very weak interactions gapping out the Dirac cone.
In our QMC simulations, for any non-zero α that we investigated, a gap appeared even for the smallest values of U probed. This suggests a weak-coupling origin of this phase. To verify it, we performed HF calculations on the same lattice model. The results, shown by the red points in Fig. 2(a), are in very good agreement with the QMC results. We also used HF to investigate the stability of the gap against changing the phase that appears in the assisted-hopping term T [32]. This phase can be gauged away, at the expense of introducing complex hopping terms in H 0 , which introduce a sign-problem to the QMC simulations. However, they do not affect the efficiency of the HF algorithm. As discussed in the SM, our analysis confirm that the onset of the QVH phase is robust and appears regardless of the phase of T .
Importantly, we find that the gap completely disappears when α = 0, in agreement with Ref. [57]. Combined with the fact that the gap onsets for small interaction values when α 0, this suggests that the origin of the gap can be understood from a mean-field decoupling of the cross-term Q T of the interaction in Eq. 2. This cross-term can be rewritten as: (5) where l and m are valley indices and the spin index is omitted for simplicity. The terms with j = i − 1 and j = i vanish after summing over different hexagons. In the weak-coupling limit, we can do a mean-field decoupling and use c † i,l c i+1,m ∝ δ lm , due to the nearest-neighbor hopping term present in H 0 . The cross-term then becomes: Thus, the cross-term of the interaction naturally induces an imaginary hopping between next-nearest-neighbors in the weak-coupling limit. As a consequence, the mean-field Hamiltonian becomes two copies (four, if we consider the spin degeneracy) of the Haldane model [65], leading to a Chern number of ±1 for the two different valleys. For this reason, we call this state a QVH phase; it is illustrated in the corresponding inset in Fig. 1 (b). We verified that our self-consistent HF calculation generates the same pattern of imaginary NNN hopping.
One of the hallmarks of the Haldane model is the existence of gapless edge modes, despite the bulk being gapped. In the QVH phase, these edge states should be valley-polarized. To probe them, we performed QMC simulations with open boundary conditions and extracted the imaginary-time Green's functions on the edge, G edge (τ) ∼ e −∆ sp τ . As shown in Fig. 2 (c), in the regime of small U (U/W = 0.25), the Green's function on the edge decays to a constant in the long imaginary-time limit, demonstrating the existence of a gapless edge mode in the QVH phase. Note that a Chern number can be defined separately for each valley l = 1 and l = 2 (with spin degeneracy). Because the valley U(1) symmetry guarantees that these two Chern numbers must be equal, the whole system is characterized by one Chern number that takes integer values, i.e. it belongs to a Z classification [66]. Fig. 2(c) also shows that, as U increases (U/W = 2), the gapless edge mode disappears, signalling a departure from the topological QVH phase. Clearly, the bulk remains gapped, as shown in Fig. 2(a). The new insulating phase is an intervalley coherent (IVC) state, which spontaneously breaks the onsite spin-valley SU(4) symmetry. In the QMC simulations, IVC order is signalled by an enhancement of the correlation function C I (k) = 1 L 4 i, j ∈ A(B) e ik·(r i −r j ) I i I j , here, the operator I i = σ (c † i,l,σ c i,l ,σ + h.c.), l l , represents an "onsite hopping" between the two different valleys, A(B) represents the A(B) sublattice. Thus, the correlation function is a 2 × 2 matrix in sublattice space, i.e.
C  Figs. 3(a) and 3(b), we show one of diagonal compent C AA I (k). The fact that the correlation function is peaked at k = Γ implies that the IVC order is ferromagnetic-like, i.e. it does not break translational symmetry. Such an onsite coupling between opposite valleys (see the corresponding inset in the phase diagram in Fig. 1 (b)) breaks the valley U(1) symmetry, and hence the SU(4) symmetry of the model. The fact that the SU(4) symmetry-breaking pattern is ferromagnetic-like is similar to recent analytical results [32,33], which focused, however, at integer fillings away from charge neutrality.
For larger values of U/W, as shown in Fig. 3, the IVC order fades away, but the system remains insulating. The new state that emerges is the columnar valence-bond solid (cVBS) insulator, characterized by the appearance of strong nearest-neighbor bonds forming the pattern illustrated in the corresponding inset of Fig. 1 (b). The onset of cVBS order is signalled by an enhancement of the bond-bond correlation function [56,57] i,l,σ c i+δ,l,σ + h.c.) and δ denoting one of the three nearest-neighbor bond directions of the honeycomb lattice (ê 1 ,ê 2 andê 3 ). For this particular calculation, e 1 was chosen.
As shown in the lower panels of Figs. 3(a) and 3(b), we find an enhanced C B (k) at momenta K and K , demonstrating that the bond-order pattern breaks translational symmetry. However, a peak of C B (k) at these momenta does not allow us to unambiguously identify the cVBS state, as the plaquette valence-bond solid (pVBS) also displays peaks at the same momenta [56,57,64,67]. To further distinguish the two types of VBS phases, we construct the complex order parameter 3 . The Monte Carlo histogram of D K is different for the two VBS phases [64,67]: for the pVBS state, the angular distribution of D K is peaked at arg(D K ) = π 3 , π, 5π 3 , whereas for the cVBS state, it is peaked at arg(D K ) = 0, 2π 3 , 4π 3 . Our results, shown in the inset of Fig. 3(a), clearly demonstrate that the cVBS order is realized in our phase diagram.
The phase boundaries in Fig. 1(b) are obtained by scanning the correlation functions C I (Γ) and C B (K) as a function of U/W for fixed values of α. Two of these scans are shown in Fig. 3, for α = 0.4 (panel (a)) and α = 0.6 (panel (b)). It is clear that, as U/W increases, in both cases the ground state evolves from QVH to IVC to cVBS and then back to IVC. Furthermore, in the strong coupling limit U/W → ∞, the IVC order C I (k = 0) is independent of α and saturates at 0.5, consistent with our analytical calculation at the charge neutrality point [63]. The transitions between IVC to cVBS are first order, as signaled by the fact that as the system size L increases, the suppression of the IVC order becomes sharper (see for instance the region around U/W ∼ 5 and U/W ∼ 11 in panel (a)). A similar sharp drop is also featured at the QVH-IVC transition (region around U/W ∼ 2.5 in panel (a)), indicating that the QVH-IVC and IVC-cVBS transitions are all first-order. It is interesting to note that, as α increases, the values of U/W for which the IVC and cVBS phases emerge are strongly reduced.

Discussion
In this paper, we employed QMC simulations, which are exact and unbiased, to obtain the phase diagram of a lattice model of TBG at charge neutrality. Our main result is that even very small interaction values trigger a transition from the non-interacting Dirac semi-metal phase to an insulating state. Upon increasing U, the nature of the insulator changes from a non-symmetry-breaking topological QVH phase, to an onsite SU(4) symmetry-breaking IVC state, to a translational symmetry-breaking cVBS phase, and then finally back to a reentrant IVC state. This rich phase diagram is a consequence of the interplay between two different types of interaction terms: a cluster-charge repulsion Q and a non-local assisted-hopping interaction T . The former is analogous to the standard Hubbard repulsion and, as such, is expected to promote either SU(4) antiferromagnetic order or valence-bond order in the strong-coupling regime. The latter, on the other hand, arises from the topological properties of the flat bands in TBG. When combined with Q , it gives rise not only to SU(4) ferromagnetic-like order, but also to correlated insulating phases with topological properties, such as the QVH phase.
While the precise value of U/W in TBG is not known, a widely used estimate is that this ratio is of order 1 [32]. Referring to our phase diagram in Fig. 1(b), this means that certainly the QVH phase and possibly the IVC phase can be realized at charge neutrality, provided that α is not too small. While some experimental probes do report a gap at charge neutrality Ref. [10,11], additional experiments are needed to establish its ubiquity among different devices and the nature of the insulating state. The main manifestation of the QVH phase would be the appearance of gapless edge states, whereas in the case of the IVC state, it would be the emergence of a k = 0 order with onsite coupling between the two different valleys.
In a more general context beyond TBG, our work offers a promising route to realize correlation-driven topological phases. As explained above, the topological QVH insulating state appears due to the cross-term in the interaction Hamiltonian that contains both Q and T . While repulsive interactions similar to the charge-cluster one are generally expected to appear in any correlated electronic system, an interesting question is about the necessary conditions for the emergence of an interaction similar to the assisted-hopping. In our case, it arises from the projection of the standard Coulomb repulsion on Wannier states that suffer from topological obstruction. The latter, in turn, is a manifestation of the phenomenon of fragile topology [68]. Thus, interacting systems with fragile topology may offer an appealing route to search for interaction-driven topological states.

Methods
We apply the projection QMC (PQMC) method [59,62,69] to investigate the ground-state phase diagram of the model. In PQMC, one can obtain the ground state wave function |Ψ 0 from projecting a trial wave function |Ψ T along the imaginary axis |Ψ 0 = lim Θ→∞ e − Θ 2 H |Ψ T . The expectation value of an observableÔ can then be calculated as We use the eigenvector of the non-interacting H 0 as the trial wave function, with electron filling fixed at charge neutrality. To be able to perform the Monte Carlo sampling, one needs to decouple the interaction term H via the Hubbard-Stratonovich transformation, . The sum is taken over the auxiliary fields s on each hexagon, which can take four values ±1 and ±2. After tracing out the free fermionic degrees of freedom, the configurational space {s (i, τ)} with size L × L × Θ is the space in which the physical observables in Eq. (7) are computed with ensemble average. We choose the projection length Θ = 2L/t and discretize it with a step ∆τ = 0.1/t. The spatial system sizes are L = 6, 9, 12, 15. More details about the QMC simulations are give in the SM [63]. Since we are interested in the ground state properties of the system, the projection QMC is the method of choice [59,62,69].

Acknowledgements
Continuing from the discussion in the Methods section in the main text, we performed Trotter decomposition to discretize Θ into L τ slices (Θ = L τ ∆τ). Each slices ∆τ is small and the systematic error is O(∆τ 2 ). After the Trotter decomposition, we have where the non-interacting and interacting parts of the Hamiltonian are separated. To treat the interacting part, one usually employs a Hubbard Stratonovich (HS) transformation to decouple the interacting quartic fermion term to fermion bilinears coupled to auxiliary fields. We use the four components auxiliary s = ±1, ±2 as shown in Eq. (8) in the main text, and for each configuration of the auxiliary field {s }, the fermions are non-interacting and can be traced out, it is the configurational average of Monte Carlo sampling that restores the interaction effects in an unbiased manner. After tracing out the free fermions degrees of freedom, we obtain the following formula with a constant factor omitted where P is the coefficient matrix of the trial wave function |Ψ T . In the simulation, we choose the ground state wavefunction of the half-filled non-interacting system (described by H 0 ) as the trial wave function. In the above formula, the B matrix is defined as and has properties B(τ 3 , τ 1 ) = B(τ 3 , τ 2 )B(τ 2 , τ 1 ), i.e. the B matrix is an imaginary time propagator, where we have written the coefficient matrix of the interaction part as V[{s ,τ }], and K is the hopping matrix from H 0 . Every hexagon contains six sites, as shown in the figure below, so our V[{s ,τ }] matrix is a block matrix, and every block contributes a 6 × 6 exponential matrix, The Monte Carlo sampling of auxiliary fields are further performed based on the weight defined in the sum of Eq. (S2). The measurements are performed near τ = Θ/2. Single-particle observables are measured by the Green's function directly and many-body correlation functions are measured from the products of single-particle Green's function based on their corresponding form after Wick-decomposition. The equal time Green's function is calculated as with R(τ) = B(τ, 0)P, L(τ) = P † B(Θ, τ).

II. ABSENCE OF SIGN-PROBLEM
At the charge neutrality point, the model is sign-problem-free, as can be seen from the following analysis. Define W σ,l,S i as the update weight of one fixed auxiliary field at the i-th hexagon, where l = 1, 2 is a valley/orbital index and σ =↑, ↓ is a spin index. From the symmetry of the Hamiltonian, W ↑,l,S i = W ↓,l,S i . Since the model is particle-hole symmetric at charge neutrality, one can perform a particle-hole transformation (PHS) only for the valley l = 2. Then one can focus on a fixed auxiliary field, and focus only on one spin flavor, such as spin up. Eq. (8) in the main text can then be abbreviated as Applying PHS for valley 2 and using the relation, we find that Eq. (S7) becomes Because of the relations above, the total weight of the model is which is a real positive number. This implies that the QMC simulations are sign-problem-free.

III. BENCHMARK WITH EXACT DIAGONALIZATION
We employ Lanczos exact diagonalization (ED) to benchmark the PQMC results. The system contains 2 × 2 unit cells of the honeycomb lattice with periodic boundary condition (16 electrons in total). We make use of symmetries, such as the valley U(1) symmetry and the total S z conservation for each valley, to reduce the computational cost of the ED. The ground state lays in the subspace with N ↑ = 4, N ↓ = 4 in both valleys, where N ↑ (N ↓ ) is the number of electrons with spin up (down) in each valley. The dimension of the ground-state subspace is about 24 million. In the PQMC simulations, we choose the linear system size L = 2 and the projection length Θ = 100/t with Trotter slice ∆τ = 0.0005/t. We compared the ground-state expectation values of H 0 and of the double occupation as a function of U/t at α = 0.3, which are shown below. The results of both methods agree very well.

IV. HARTREE-FOCK METHOD
We solve the full Hartree-Fock Hamiltonian (H = H 0 + H H F ) self-consistently using that where κ, λ = {ilσ}, U is the unitary tranformation diagonalizing H, γ's are the eigenvectors and f (E , µ) is the Fermi-Dirac distribution of the excitation energies, E . We explicitly write the dependence on the chemical potential, µ, as we iterate this value to fulfil N −1 f (E , µ) = ν, where ν is the filling. We compute results at charge neutrality (ν = 0.5) with a total of 600 lattice sites and periodic boundary conditions. The calculations are fully unrestricted; thus we iterate all (4 × 600) 2 mean-fields, and define convergence by the condition that |∆E | < N × 10 −10 , where ∆E is the change of the excitation energies from one iteration to the next, and N is the total number of states (4×600). In Table S1 we present an example of the HF calculations, displaying results for t = 1 ,α = 0.45 and U/W = 0.5. We set the temperature T = 2.5 · 10 −5 in all computations. The values in the table are the renormalized mean-fields. It is evident that all hoppings within each hexagon are renormalized due to the interactions. The simple hopping renormalizations, however, do not open a gap in the Dirac cones. The gap is generated directly by the mean-fields c † i,1,σ c i±2,1,σ = − c † i,2,σ c i±2,2,σ = ±0.0914i, which explicitly display a spin degenerate quantum valley Hall (QVH) phase, as illustrated in Fig. 1(b) of the main text.
In Fig. S2(a) we show the single-particle gap in the QVH phase with α = 0.45 for several interaction strengths. Fig. S2(b) displays the corresponding band structures. The Dirac cone at K in the bare bands is immediately gapped out when including interactions. The renormalization initially flattens the bands with a significant gap at all high-symmetry points. As U increases the valence band gradually develops a peak at Γ while it is pushed down correspondingly at K. This behavior results in a gradual shift of the maximal gap value from Γ to K.  The topological nature of the QVH phase is manifested by edge states. To verify the emergence of said states we open the boundaries in the system and compute a self-consistent result with parameters as in Table S1 (t = 1 ,α = 0.45, U/W = 0.5, T = 2.5 · 10 −5 and N = 4 × 600). We find clear evidence of edge states as seen in Fig. S3.
Finally we present results obtained by implementing the interaction terms found in Ref. 32. As mentioned in the discussion section of the main text, we are able to solve this model within the HF approach as it does not suffer from sign problems. The assisted hopping reads The other terms, Q and H 0 , remain unchanged. To reach this expression for T , we have performed the following gauge transformation, c ilσ −→ e iθ l /2 c ilσ , i odd c ilσ −→ e −iθ l /2 c ilσ . i even (S17) The transformation introduces phases in H 0 effectively causing t to become complex. We set the phases according to Ref. [32], that is θ 1 = −θ 2 = 0.743π. The renormalized mean-fields are presented in Table S2, where we have performed the inverse gauge  transformation for direct comparison with Table S1. Input parameters are the same as those used to generate Table S1. The result is consistent with the values presented in Table S1 and clearly also features a QVH phase.  Table S1 and open boundary conditions. U/W = 0.50 c i,1,σ c i,2,σ c i±1,1,σ c i±1,2,σ c i±2,1,σ c i±2,2,σ c i+3,1,σ c i+3,2,σ c †