Understanding Symmetry Breaking in Twisted Bilayer Graphene from Cluster Constraints

Twisted bilayer graphene is an exciting platform for exploring correlated quantum phases, extremely tunable with respect to both the single-particle bands and the interaction profile of electrons. Here, we investigate the phase diagram of twisted bilayer graphene as described by an extended Hubbard model on the honeycomb lattice with two fermionic orbitals (valleys) per site. Besides the special extended {\it cluster interaction} $Q$, we incorporate the effect of gating through an onsite Hubbard-interaction $U$. Within Quantum Monte Carlo (QMC), we find valence-bond-solid, N\'eel-valley antiferromagnetic or charge-density wave phases. Further, we elucidate the competition of these phases by noticing that the cluster interaction induces an exotic constraint on the Hilbert space, which we dub {\it the cluster rule}, in analogy to the famous pyrochlore spin-ice rule. Formulating the perturbative Hamiltonian by projecting into the cluster-rule manifold, we perform exact diagonalization and construct the fixed-point states of the observed phases. Finally, we compute the local electron density patterns as signatures distinguishing these phases, which could be observed with scanning tunneling microscopy. Our work capitalizes on the notion of cluster constraints in the extended Hubbard model of twisted bilayer graphene, and suggests a scheme towards realization of several symmetry-breaking insulating phases in a twisted-bilayer graphene sheet.


I. INTRODUCTION
Twisted bilayer graphene (TBlG) has emerged as a versatile platform for studying competing phases in a system with strong correlations.Experiments have reported correlated insulating phases [1][2][3][4][5][6][7][8][9][10], in other words insulating phases where band theory predicts a metal, as well as nodal (unconventional) superconductivity [2,3,7,11,12].The relationship between the insulators and the superconducting phases is still unclear: either the insulator can be viewed as the parent state, which upon doping becomes superconducting, or the insulator and superconductor are competing phases, with different underlying mechanisms.
Experiments have addressed this question by studying the dependence of the phase diagram of TBlG on the electronic screening.Increasing the screening decreases the range and strength of the Coulomb interaction.In Ref. 13, the screening was tuned by varying the electron density in a metal close to the TBlG, while in Refs.8 and 14 distinct devices with different gate distances were investigated.Strikingly, the experiments observed that the insulating phases weaken or disappear, while superconductivity survives even as the screening is increased.In our work, we focus on understanding the effect of the modified screening on the correlated insulators.
TBlG sits at the intersection of two paradigmatic models of strongly correlated phases.On the one hand, the flat Chern bands are reminiscent of the Landau levels of the quantum Hall effect, an analogy that can be made precise in an idealized model of TBlG [15].In this language, the correlated insulators are generalized quantum Hall ferromagnets that may exhibit intervalley coherence (IVC) [16][17][18][19].On the other hand, the proximity of correlated insulating and unconventional superconducting phases as well as linear-in-T resistivity [20,21] suggests a connection to the phenomenology of the Hubbard model used to model cuprate superconductors.In that language, the correlated insulators are valence bond solid (VBS) or quantum valley Hall (QVH) phases [22][23][24].Aiming at studying a strong coupling theory with local constraints, in this work we consider the Hubbard model of TBlG with realistic extended interactions.In the phase diagram of this model, we identify the emergent strongly-correlated phases and discuss how to tune the interactions experimentally across the phase diagram.We emphasize the similarity of these extended interactions' ground states to the iconic spin-ice manifold in the pyrochlore lattice [25] and develop an intuitive perturbation theory for hopping terms.Finally, we provide possible STM signatures of these correlated phases.
The special extended interactions emerge due to the characteristic 'fidget spinner' form [26][27][28] of the maximally localized Wannier functions.This leads to significant overlaps between Wannier orbitals on neighbouring sites implying that longer-range Hubbard interactions need to be included.The Coulomb interaction leads to the cluster-interaction Q, where the onsite, nearest-neighbour, next-nearest-neighbour and next-to-next-nearest-neighbour interactions satisfy the ratios 3:2:1:1.This interaction can be written as a sum of perfect squares fixing the total charge per hexagon (cluster), thus the name.Besides, to mimic the screening we introduce an on-site interaction U that we tune separately [29].
We consider a model with two orbitals on each moiré honeycomb lattice site representing the two valleys of the microscopic graphene lattice.The dynamics of these FIG. 1. Phase diagram in the (Q/t, U/t) plane.The fixed-point states are sketched in the respective regions of the phase diagram.The (U, Q) pairs mark the phase boundary obtained by extrapolating the QMC data to the thermodynamic limit.The error bars indicate the direction of a scan in the plane.The uncertainty is ∆Q/t = 1.0 for horizontal and ∆U/t = 0.01 for vertical scans.Some transition points in these phase diagrams are already known from the literature: the U/t = 0 axis at Ns = 2 was studied in Ref. [22], the Q/t = 0 axes at Ns = 1 and 2 were studied in Refs.[30,31] orbitals is given by a (nearest-neighbor) hopping with strength t, which sets the energy scale.We study both the case of a single spin species N s = 1 as well as the case of two spin species N s = 2, relevant to studying phases with non-trivial pseudo-spin structure.We study the model at half-filling.The physical system always has spin, and the case N s = 2 at half-filling corresponds to studying the charge-neutrality point ν = 0 of the physical system.The model with N s = 1 is relevant for the spin-polarized phases of the physical system at ν = ±1.A spin-polarized phase at ν = −2 only has one spin species occupied and the other spin-species is at halffilling, therefore N s = 1 at half-filling is the relevant case to study.On the other hand, a spin-polarized phase at ν = +2 has one spin species fully occupied (and therefore inert) while the other spin species is at half-filling.Again, the model with N s = 1 at half-filling is relevant.We observe numerous phases in the (Q, U ) phase diagram, including the weakly-interacting Dirac semimetal phase (DSM) Kekulé VBS, Néel and charge-density wave (CDW).
Figure 1 summarizes the phase diagram we obtain within Quantum Monte Carlo (QMC) [32].The DSM phase is characterized by vanishing susceptibilities and single-particle gap.In the N s = 1 case, the three valley Néel orders τ v are present, where τ v are the Pauli matrices with v ∈ {1, 2, 3} denoting the valley degree of freedom.These order parameters form the adjoint representation of the SU (2) group and their susceptibilities are degenerate.The phase diagram is shown in Fig. 1 (a), where we depict the τ 3 order.In turn, in the Kekulé state, namely columnar VBS (cVBS) emerging in the moiré scale, bond singlets and plaquettes are formed.At the same time, in the N s = 2 setup, we do not observe a Néel state, but rather two Kekulé states, cVBS and plaquette VBS (pVBS).The phase diagram is shown in Fig. 1 (b).
In the valley Néel phase, where three degenerate orders are expected, we consider how the possible perturbations to the Hamiltonian could break the degeneracy.We find that the leading perturbation, next-to-nearest-neighbor hopping, favors the intervalley orders such as τ v with v ∈ {1, 2}.We also show that this phase has distinct features in the local electron density observable in scanning tunneling microscope (STM) experiments.
Importantly, the strong cluster interaction Q fixes the number of electrons to 6N s per each hexagon, which we call the cluster rule, and the states satisfying this constraint -the cluster-rule manifold.Since the hexagons share corners, the cluster-rule manifold cannot be decomposed into a product of simple local Hilbert spaces, similarly to the iconic pyrochlore spin ice [25], where the total magnetization per a (corner-sharing) tetrahedron is constrained to zero.This is in contrast to the case of the strong on-site repulsion U , which limits the occupation to one fermion per site and turns the Hubbard model into the Heisenberg model with simple on-site degrees of freedom.Within this cluster-rule manifold, we develop an intuitive perturbation theory and supplement the QMC results with exact diagonalization for N s = 1.
The paper is organized as follows.In Section II, we introduce the model and the relevant interactions.In Section III, we study the model within several numerical techniques, address breaking degeneracy between Néel states and show possible STM images.Finally, in Section IV we discuss the obtained results.

II. MODEL
The Wannierization of TBlG conducted in Ref. 27 leads to orbitals centered at the sites of a honeycomb lattice on the moiré scale.The sites of the honeycomb lattice represent the AB and BA stacked regions of the bilayer.Similarly to Ref. 27, we consider the fermionic model on the honeycomb lattice with two orbitals (valleys) and N s spins, with N s = 1 or 2. In the case of TBlG, the orbitals are delocalized, thus an electron on a site has its wave function density concentrated in the three blobs located in the centers of the three hexagons adjacent to the site, which correspond to the AA/BB regions of the bilayer.For simplicity, we consider a real nearest-neighbor kinetic term with hopping strength t, which yields the SU (2N s )symmetric Hamiltonian where ⟨i, j⟩ denote nearest-neighbor sites, σ enumerates the N s spin species and τ = +, − enumerates the two valleys.The operator ni = σ, τ ĉ † iστ ĉiστ measures the full charge at a site and the hexagon charge operator is given by n = (1/3) i∈ ni .The factor (1/3) stems from the delocalization of Wannierized fermionic orbitals into three peaks.
The estimates for the characteristic Coulomb energy Q = e 2 /(0.28LM ) ∼ 220 meV (here, L M is the moiré length scale) and the bandwidth ∆ ∼ 6t ∼ 2 meV given in Ref. 27 for θ ∼ 1.07 • suggest that the ratio Q/t can reach hundreds.The cluster interaction depends on the gate distance d, as shown in Fig. 2 in agreement with Ref. 8.
In the absence of the on-site interaction term U , the cluster interaction can be seen as an extended Hubbard model with the long-range repulsive contributions decaying as 27, where the interaction strength stems solely from the number of overlapping wave function density centers.In the experimental setup with the interaction screened by two parallel metallic gates placed at the distance d to the TBlG sheet, this cluster-interaction limit corresponds to the case d ≪ L M , where L M is the TBlG moiré lattice spacing [27,33].When the metallic gates are moved away from the TBlG sheet, the cluster interaction is supplemented with the screened non-local Coloumb interaction between the non-overlapping blobs.This modification of the cluster interaction violates the 3/2/1/1 extended interactions ratio, which we mimic by introducing the on-site interaction U .In Fig. 2, we show the ratio U/Q as a function of d/L M .Importantly, in an experiment with metallic gates, both positive and negative values of U/Q can be realized.In Fig. 1 (a-b), we draw the possible lines that can be induced by the dependence of Q(d)/Q(d = 0) and U (d)/Q(d) on d.The fully-screened interaction strength Q(d = 0)/t depends on the twist angle θ.Notably, these lines pass through all phases, making the variation of d a suitable way to access the observed symmetry-brkoen phases experimentally.

A. Order parameters
In this Section, we sketch the order parameters emerging in the phase diagram of the Hamiltonian Eq. (1).

Kekulé valence-bond solids
In the Kekulé VBS phase, reported previously in Ref. [22], the order is given by the operator where N b = 3L 2 is the number of bonds in the lattice, ρ i enumerates the nearest neighbors of site i, and ξ ∈ {0, . . ., 2N s − 1} labels the flavors.This order transforms as the E representation of the C 3 group and has the spatial momenta K M and K ′ M , the Dirac points in the moiré Brillouin zone.
The matrix representations M αβ of the ±K M Kekulé orders Eq. ( 2) are related by complex conjugation and have degenerate susceptibilities.In Ref. 22, it was shown that these two orders condense in two real-valued superpositions ( M+KM + M−KM )/ √ 2 and ( M+KM − M−KM )/(i √ 2), the pVBS and cVBS, respectively.The determination of which phase is realized can be done by plotting the histogram of the ⟨ M+KM ⟩ measurements within QMC in the complex plane.The characteristic distributions of the order parameter ⟨ M+KM ⟩ are shown in Fig. 1 (c-d).

Néel antiferromagnets
In the Néel AFM phase, the symmetry breaking may take place due to the condensation of the orders where ĉ † i is the collection of 2N s creation operators at the site i, l(i) = 0, 1 is the sublattice index of the site i, and t v,s is one of the (2N s ) 2 − 1 Lie algebra generators of the SU (2N s ) group.Here, in the N s = 2 case, we choose the standard representation where 0 ⩽ v, s < 4 enumerate the Pauli matrices with v + s > 0. For the N s = 1 case, t v = τ v with v > 0.

Charge-density wave
Lastly, in the CDW the order is described by condensation of

B. Quantum Monte Carlo study
The Hamiltonian Eq. ( 1) can be studied within the Quantum Monte Carlo (QMC) approach [22], which is sign-problem free at U, Q > 0. At U < 0, the approach is only sign-problem free at Q = 0, which hinders the study of the phase diagram for negative U .
We consider the equilateral L × L clusters with L = 3, 6 9 12 and 15 and vary the Trotter step t δτ = min (0.1, t/Q, t/U ) .This choice of δτ allows us to keep the Trotter errors under control even at strong interaction [34].The 4-valued Hubbard-Stratonovich field is introduced using the 4-th order O(δτ 4 ) decomposition with the errors negligible as compared to the Trotterization errors [35].
A generic order parameter has the form where α = (i α , τ α , s α ) is the composite site-flavor index.This operator is normalized as M * αβ M βα = 1.The tendency to establish a non-zero value of the order ⟨ M ⟩ ̸ = 0 can be quantified by the zero-frequency suscep-tibility where we singled-out the four-point particle-hole vertex operator Treating (γ, δ) and (α, β) as the composite in-and outindices, we view Γ (αβ),(γδ) as a matrix, whose eigenvectors are the operators transforming as irreducible representations of the translation and point-group symmetries and the eigenvalues are the susceptibilities.The irreducible representation of the eigenvector M with nonzero extrapolated susceptibility lim For all parameters studied in this work, we find that the leading eigenvector is always one of the orders defined in Section III A 1.Moreover, we reproduce the expected ((2N s ) 2 − 1)-fold and 2-fold degeneracies for the M v,s Néel and the M±KM Kekulé orders, respectively [36].
We perform QMC simulations at βt = 2L and check for selected parameters that the conclusions are temperature-insensitive at βt = 3L.The resulting QMC phase diagrams are shown in Fig. 1 (a-b).In Fig. 3 (a) we show the susceptibilities' extrapolations in the N s = 1 case.Here, the extrapolated Néel susceptibility is finite for (U/t, Q/t) = (2, 80), indicating the valley AFM phase, while the Kekulé susceptibility is finite in the cVBS phase for (U/t, Q/t) = (0, 100).In the DSM phase at (U/t, Q/t) = (0, 20), all susceptibilities extrapolate to zero.In Fig. 3 (b), where N s = 2 is shown, the Kekule susceptibilities are finite in the pVBS phase at (U/t, Q/t) = (0, 60) and in the cVBS phase at (U/t, Q/t) = (6, 60), while they extrapolate to zero in the DSM phase at (U/t, Q/t) = (0, 40).

C. Perturbation theory
The large Q/t interaction strengths of the phase transitions in Fig. 1 (a-b) suggest construction of an effective theory in the limit of the strong cluster interaction (Q ≫ t, U ), where the low-energy manifold is described by the states having exactly 6N s fermions in each hexagon.Notably, unlike the construction of the t-J model in the regime of strong U , here the resulting cluster-rule manifold can not be written as a product of Hilbert spaces of some emergent local degrees of freedom, which is similar to the pyrochlore lattice [25].The cluster-rule manifold is separated by the gap ∆ Q = Q/9 from the rest of the states in the full Hilbert space.Within this manifold, the on-site term U is diagonal, while the kinetic term t can be treated perturbatively where K, Q are the kinetic and cluster terms in Eq. ( 1), and P projects onto the orthogonal complement of the cluster-rule manifold.In this manifold, the effective Hamiltonian reads: Here, the first term, exchanges the fermions with the flavors ξ 1 , ξ 2 on the nearest-neighbor sites i, j, which is shown in Fig. 4 (a).
The second term ĥhexagon = ξ1,ξ2,ξ3 moves three fermions with the flavors ξ 1 , ξ 2 , ξ 3 over three non-adjacent bonds in a hexagon in the same direction, as shown in Fig. 4 (b).Note that each of these terms preserves the cluster rule in all hexagons.These terms annihilate a state if a move is not possible.To illustrate these moves, in Fig. 4 (c) we show a basis element satisfying the cluster rule.The figure illustrates the dynamics induced by the terms ĥbond and ĥhexagon within the cluster-rule manifold.
We treat the Hamiltonian Eq. ( 9) within exact diagonalization (ED) [37] in the N s = 1 case on the 3×3 hexagonal lattice.The cluster-rule manifold contains 5'018'802 states, which is a considerable reduction from the full Hilbert space with 9 + 9 fermions of the size ∼ 3 × 10 9 .Unfortunately, the cluster-rule manifold on the 3 × 3 lattice with N s = 2 has the size of O(10 13 ) and is out of reach.
Figure 5 (a) shows the low-energy level structure at U/t = 0.2 as a function of Q/t.We determine the phase diagram on the Q/t axis by considering the quantum numbers of the low-energy excited states, namely their spin S, momentum k and their transformation property under the point group D 3 .In the regime of small Q/t, the lowest excitations above the (featureless) ground state, are singlets S = 0 with k = K M /K ′ M , transforming as the E representation under the D 3 group, which matches the symmetry-breaking pattern of the Kekulé state.In the large-Q/t regime, we observe the magnetic tower of states signalling breaking of the SU (2) symmetry and establishing a Néel order [38].Following Refs.39 and 40, we identify the transition between the magnetic Néel phase at large Q/t and Kekulé phase at small Q/t with the crossing between the lowest-lying triplet and the symmetry-breaking singlet Kekulé state.
Negative U/t = −0.2favours doubly-occupied or free sites.There, the phase transition is pinned by the crossing between the Kekulé and the CDW (S = 0, k = 0, A 2 )   excited states, signalling the k = 0 CDW, as we show in Fig. 5 (b).With these transition criteria, in Fig. 5 (c) we show the resulting phase transition lines between the Kekulé and the CDW phases.At U > 0 the results qualitatively agree with the QMC phase diagram shown in Fig. 1 (a).The actual interaction strengths differ due to the significant finite-size effects present in L = 3 perturbation theory.Importantly, the region of U < 0 is not accessible for the exact QMC study due to a severe sign problem.

E. Fixed-point states analysis
The absence of a Néel phase for N s = 2, which is in contrast to the N s = 1 case, can be explained by constructing the fixed-point states of the orders observed within QMC and ED.The corresponding orders are depicted in Fig. 1 (a-b).A fixed-point state is obtained as a ground state of the Hamiltonian Ĥ = M , where M is one of the order parameters, i. e., the condensed order dominates the initial Hamiltonian.Finally, we project this state onto the cluster-rule manifold by removing the wave function components violating the cluster rule.
The possible Néel SU (4) state (for simplicity, we consider the σ 0 ⊗ τ 3 order) reads where the first (second) operator creates fermions in the τ = +1(−1) valley on the A(B) sublattice.In the SU (2) spinless case, the spin degree of freedom in Eq. ( 13) is omitted.
For the cVBS phase, the state reads where b enumerates the N b /3 dimerized bonds between sites (α b , β b ), while c enumerates 2Ns Ns possible ways to put N s fermionic species at the site α b and the remaining N s species at the site β b , forming the state |v α b c ⟩ ⊗ |v β b c ⟩.This state has exactly N s electrons per site and thus satisfies the cluster rule.
We compute the energies of these fixed-point states according to Hamiltonian Eq. ( 9), yielding We see that at N s = 1, the Néel and cVBS energies are equal, E Néel = E Ns=1 cVBS , which opens the prospect to realize the SU (2) valley Néel state in Fig. 1 (a) due to higher-order quantum fluctuations.On the contrary, at cVBS , and thus the SU (4) Néel state in Fig. 1 (b) is never a ground state [41].

F. Breaking the degeneracy between the Néel states
We observed that in the spinless case N s = 1, the valley Néel phase appears in the phase diagram.The three Néel states are degenerate within the given SU (2)-symmetric model Eq. ( 1).However, additional perturbations in the Hamiltonian can break this symmetry, which may lead to a splitting of the degeneracy between these Néel orders.
In Eq. ( 1), we neglected additional long-range hoppings and longer-range interactions emerging due to the delocalized nature of the Wannier orbitals [27].The next-toleading hopping obtained in Ref. [27] is the fifth-nearestneighbor hopping with the respective Hamiltonian where ⟨⟨i, j⟩⟩ denotes the fifth-nearest-neighbor pairs.Crucially, while the nearest-neighbour hopping t can always be chosen to be real via an appropriate choice of gauge, t 5 is generally complex [42].With a complex hopping, the symmetry of the model is broken down from valley SU (2) symmetry to U (1) valley charge conservation.
Since this term hops a single particle, there is no contribution in the first-order perturbation theory.The second-order perturbation theory (including only the imaginary component of the t 5 hopping since this is the piece which breaks the symmetry) yields Therefore, adding the 5-th nearest-neighbor hopping favours inter-valley Néel orders, i. e. the system is an easy-plane valley-antiferromagnet.

G. Graphene scale electron density distribution
STM topography experiments have recently revealed that the correlated insulating phases in TBlG [43] and twisted trilayer graphene [44] exhibit a Kekulé pattern on the microscopic graphene scale.This is a signature of certain types of intervalley coherent (IVC) states [45,46] and is consistent with predictions from Hartree-Fock on the continuum model [17,47,48].
Since STM has been shown to be a powerful tool for identifying the nature of symmetry breaking in TBlG, we plot in Fig. 6 for N s = 1 the potential STM images that would appear on the microscopic graphene scale of the bilayers, when the respective orders we found stabilize on the moiré scale.In particular, we show the local electron density in one of the layers for the CDW and for three different Néel orders, all of which break moiré translational symmetry.The τ 1 and τ 2 Néel orders allow for superposition of electrons from the two microscopic graphene valleys and this leads to a √ 3 × √ 3 increase of the size of the microscopic graphene unit cell.In the Fourier transform of the STM images, the graphene-scale translational symmetry breaking shows up as additional peaks besides the graphene Bragg peaks.Using the nomenclature of Ref. 43 there can be bond as well as site Kekulé order.Intrasublattice IVC leads to a site Kekulé pattern whereas intersublattice IVC leads to a bond Kekulé pattern.In addition, the τ 1 and τ 2 Néel orders both break C 3 symmetry.The degeneracy between the different Néel orders is lifted by the t 5 term which favours the intervalley Néel order (with a Kekulé pattern).

IV. DISCUSSION AND OUTLOOK
In this work, we considered an extended Hubbard model for twisted bilayer graphene with two valleys and N s spins on the honeycomb lattice.The topology of the bands of twisted bilayer graphene leads to delocalized Wannier orbitals and therefore the onsite Hubbard interaction U is subdominant.Instead, the dominant interaction term is a special cluster interaction Q within each hexagonal plaquette of the honeycomb lattice.Experimentally, the range of the Coulomb interaction and therefore the competition between U and Q can be tuned by modifying the distance between the sample and the gates.While previous works have focused on this model for U = 0 [22,27], we extend the study to non-zero U and show that the presence of this term can stabilize a new phase.We studied the phase diagram of this model in U -Q space using QMC.While the QMC is limited to U > 0 due to the sign problem, we make progress on the U < 0 side of the phase diagram by performing ED within the manifold of states satisfying a cluster rule for large Q.The different competing phases we find include a Dirac semimetal, as well as three symmetrybroken phases (CDW, Néel and VBS).In particular, the Néel phase is only stabilized for nonzero U .
We explored both N s = 1 and N s = 2 and find that the phase diagrams are different: In particular, the Néel phase only appears for N s = 1.For N s = 1, we find good qualitative agreement between the QMC and the perturbation-theory approaches in the U > 0 case, where there are solutions from both methods.However, the perturbation theory is too computationally expensive in the N s = 2 case and therefore the U < 0 part of the phase diagram remains unexplored in this case.This would be an interesting regime to study in future work, for example using density-matrix renormalization group methods.
Besides, the cluster rule of having exactly 6N s fermions per hexagon, while the hexagons share corners, appears very similar to the famous pyrochlore spin-ice rule [25], where the tetrahedra are also corner-sharing.There, this 'frustrated' local charge conservation allows to reformulate the model in terms of a frustrated gauge theory hosting a U (1) algebraic spin liquid.A similar U (1) gauge theory can possibly be built in this model with the clus-ter interaction.However, this U (1) gauge theory would be confined in two dimensions [49], in accord with us not finding any spin liquid phases.The construction of this gauge theory is beyond the scope of this work.
Finally, one of the key questions is what symmetrybreaking phases are realized in the experimental system.Many of the experiments performed on TBlG perform transport measurements and as such are unable to distinguish between the different patterns of symmetry breaking: these states all appear as otherwise featureless insulating states in transport.However, recently advances have been made in using STM in order to image TBlG [43,44].This allows one to image the form of translational symmetry breaking that is characteristic of the different phases.We computed the STM pattern for the CDW and Néel phases found in the QMC and showed that they can indeed be distinguished with such a measurement.
In our work, we studied one of the two main models employed to describe TBlG, namely a Hubbard-like model with fidget spinner Wannier orbitals descending from the symmetries of TBlG [26][27][28].This Hubbard-like model allows us to utilize the local cluster-rule constraints and access the physics of TBlG from the strong-coupling regime.This real-space description stands in contrast to the often-used momentum-space description in terms of the flat bands of the continuum model [50] whose strong-coupling ground states are generalized quantum Hall ferromagnets [16].The mechanism behind these ground states is the Stoner ferromagnetism due to the exchange interactions.On the other hand, Hubbard-like models often give rise to antiferromagnetic states such as the Néel state in our model, in a way that can be understood using Anderson superexchange in the simplest case.Which of these two descriptions is most appropriate depends crucially on the range of the interactions as well as the comparison between the exchange and superexchange energy scales.One of the great strengths of TBlG and van-der-Waals materials in general is their tunability, via gating, substrate potentials and twist angle.In our work we showed how this tunability allows access to different phases in the phase diagram of the Hubbard model studied.In future, it may be possible to tune the system between the Hubbard model regime and the continuum model regime.
FIG.1.Phase diagram in the (Q/t, U/t) plane.The fixed-point states are sketched in the respective regions of the phase diagram.The (U, Q) pairs mark the phase boundary obtained by extrapolating the QMC data to the thermodynamic limit.The error bars indicate the direction of a scan in the plane.The uncertainty is ∆Q/t = 1.0 for horizontal and ∆U/t = 0.01 for vertical scans.Some transition points in these phase diagrams are already known from the literature: the U/t = 0 axis at Ns = 2 was studied in Ref.[22], the Q/t = 0 axes at Ns = 1 and 2 were studied in Refs.[30,31].The squares identify the points of the phase diagram where we later show infinite-volume extrapolations of susceptibilities within QMC.The dashed violet lines indicate the possible paths in the phase diagrams induced by the dependence of the interaction parameters on the gate distance d shown in Fig. 2. (a) The spinless Ns = 1 case.(b) The spinful Ns = 2 case.(c-d) At Ns = 2, distributions of the complex order parameter ⟨ MK M ⟩ = Tr Ĝ MK M for Q/t, U/t = (84.4,0.24) and (96.6, 0.24), respectively.These parameters are marked as stars in panel (b).The distributions are characteristic for the (c) cVBS and (d) pVBS Kekulé orders.

FIG. 2 .
FIG. 2. (Left axis) Ratio of the on-site potential correction U to the cluster interaction Q as a function of d/LM, where d is the distance to the metallic gates from the TBlG sample.(Right axis) Dependence of the screened cluster interaction Q(d) on the gate distance d.The details of both calculations are given in Appendix B.

6 FIG. 4 .
FIG. 4. Two moves in the cluster-rule perturbation theory.(a) Bond flip exchanging fermions of flavors ξ1 and ξ2 on sites i and j.(b) The hexagon flip process that moves three fermions with the flavors ξ1, ξ2, ξ3 over three non-adjacent bonds in a hexagon in the same direction, such that the cluster rule is respected.(c) An example of a cluster-rule configuration for the Ns = 1 case.Left (right) colored semicircles indicate presence of a spin up (down).Note that in each hexagon, there are exactly 6 fermions.The arrows along the bonds indicate the possible single-particle hops between adjacent sites that are allowed by the Pauli principle.The round arrows in the centers of the hexagon indicate possible (b)-moves.Namely, there are two (b)-moves possible in the upper-right hexagon, one (b)-move possible in the lower hexagon, and no moves possible in the upper-left hexagon.

FIG. 5 .
FIG. 5. (a)Spectrum obtained within exact diagonalization of the Hamiltonian Eq. (9) at U/t = 0.2.Here, we plot the energy gap between an excitation (and mark its quantum numbers) and a featureless ground state (with quantum numbers S = 0, k = 0, A1).The lowest excitations in selected symmetry sectors are connected with lines, and their quantum numbers are given.Other excitations are shown as simple black squares.(b) The case U/t = −0.2.(c) The resulting phase diagram with phase transitions pinned from the crossings of the excited states' energies.

FIG. 6 .
FIG.6.Upper row: Local electron density ρ(r) of the CDW and the τ 1 , τ 2 and τ 3 Néel orders, which could be measured using STM.The images show the local electron density in one of the graphene layers on the microscopic graphene scale.White dots show the positions of the sites of one of the graphene layers.Bottom row: Fourier transforms ρ(q) of the respective STM images.The central peak at q = 0 is removed for visibility.The six graphene Bragg peaks are circled, the additional peaks in the Fourier transform for the τ 1 and τ 2 orders are a result of the √ 3 × √ 3 translational symmetry breaking in the presence of intervalley coherence and correspond to momenta K − K ′ .