Kagome chiral spin liquid in transition metal dichalcogenide moir\'e bilayers

At $n=3/4$ filling of the moir\'e flat band, transition metal dichalcogenide moir\'e bilayers will develop kagome charge order. We derive an effective spin model for the resulting localized spins and find that its further neighbor spin interactions can be much less suppressed than the corresponding electron hopping strength. Using density matrix renormalization group simulations, we study its phase diagram and, for realistic model parameters relevant for WSe$_2$/WS$_2$, we show that this material can realize the exotic chiral spin liquid phase and the highly debated kagome spin liquid. Our work thus demonstrates that the frustration and strong interactions present in TMD heterobilayers provide an exciting platform to study spin liquid physics.

In this Letter, we show that the generalized Mott-Wigner states at n = ±3/4 filling, reported for WSe 2 /WS 2 bilayers [12,13], can realize both a chiral spin liquid [26,27] and the kagome spin liquid (KSL) [28][29][30][31][32][33]. At this particular filling, electrons are localized on an effective kagome lattice, which is known for its high degree of geometrical frustration. Here, we demonstrate how realistic model parameters lead to an effective spin model on this kagome lattice and investigate the model using extensive state-of-the-art density matrix renormalization group (DMRG) simulations [34,35]. The tunability of TMD bilayers -changing twist angle, gate tuning, material and dielectric environment choice, pressure, and so forth -thus allows for a systematic pursuit of spin liquid phases [36][37][38][39].
Model.-The moiré pattern of TMD heterobilayers is formed due to the lattice mismatch between the two layers, where the effective moiré length can be tuned by adjusting the twist angle. The interlayer band alignment ensures that the first conduction or valence flat band is completely localized in one of the layers. Based on our earlier work [22], we describe the resulting flat bands by a spin-orbit coupled extended Hubbard model on the triangular lattice, where jk denotes nearest-neighbor sites. We include the hopping matrix elements t jk up to third-nearest neighbor, where φ jk represent their phases induced by spin-orbit coupling. When the nearest-neighbor repulsion is sufficiently large, a charge density wave is stabilized at commensurate fillings. In particular, at n = ±3/4 filling, the charge order forms a kagome lattice, as shown in Fig. 1(c) [40,41]. In the Supplemental Material (SM) [42][43][44], we show, using a simple mean field theory, that the charges are almost completely localized on the kagome lattice when V /t 1 ≥ 5. In order to study the spin degree of freedom of these localized charges, we derive an effective spin model on the kagome lattice, starting from the extended Hubbard model of Eq. (1) [45,46]. In our strong coupling expansion, we keep all terms of second and third order in the hoppings, and up to fourth-order contributions proportional to the nearest-neighbor hopping t 1 . We employ the method introduced in Ref. [47] and the derivation and coefficients of the model are provided in detail in the SM [42].
The resulting spin model is given by  (2) for U/t1 = 75 and V /t1 = 10.5 -corresponding to the interactions for ε ≈ 9.5 in WSe2/WS2 at θ = 0 • -as a function of t2/t1 and t3/t1. The appearing phases are the chiral spin liquid (CSL), the kagome spin liquid (KSL) connected to the ground state of the nearest-neighbor Heisenberg Hamiltonian, a valence bond crystal (VBC), and cuboc1 * , q = (0, 0) * and √ 3 × √ 3/ferromagnetically ordered phases. The √ 3 × √ 3 and ferromagnet are energetically degenerate. The red dot indicates the t2/t1 and t3/t1 values for WSe2/WS2 at θ = 0 • . Arrows show how this point would shift qualitatively when tuning twist angle, pressure or changing the effective mass by a different material choice. (b) Phase diagram for t2/t1 ≈ 0.145 and t3/t1 ≈ 0.08 -corresponding to WSe2/WS2 at θ = 0 • -as a function of U/t1 and V /U . Here, the cuboc2 * phase emerges as an additional magnetically ordered phase. The ratio of V /U can be tuned by the gate distance while U/t1 changes with different dielectric environment. Crucially, the system can be tuned into the CSL by mereley changing the gate distance across almost the entire range of interaction strength. The hatched area denotes the region in which J1 is negative (ferromagnetic) and the red dot indicates the interaction values of panel (a). The data underlying the phase diagrams has been obtained on an infinite YC8 cylinder. (c) Extended Hubbard model at 3/4 filling with all charges localized on a kagome lattice. The unit cell of the charge ordering is indicated by the blue-lined box. When a spin-↑ particle hops from a site j to a site k along the direction of an arrow, it picks up a phase of φ jk . (d) Interactions of the resulting spin model on the kagome lattice.
where theφ ij are linear combinations of the φ jk phases from Eq. (1), and we neglected very small four-spin terms. The sum over ij runs over neighbors as illustrated in Fig. 1(d). The spin model of Eq. (2) contains XXZ and Dzyaloshinskii-Moriya (DM) terms caused by the nonzero phases φ jk in the extended Hubbard model. These phases are constrained by symmetry [22] and translate into the phases for the spin model as follows: φ 1 = 4π/3,φ 2 = 0, andφ 3 = 2π/3 for nearest, nextnearest and next-next-nearest neighbor couplings, re- spectively. This combination allows for a local three sublattice gauge transformation (a spin rotation in the x-y plane) [42] that brings the model into SU (2)-invariant form, hence, the model still exhibits a hidden SU (2) symmetry [48][49][50]. As a result, the structure of the phase diagram of our kagome spin model (2) coincides exactly with that of an SU (2)-invariant J 1 -J 2 -J 3 -J 3 model on the kagome lattice, despite the presence of the DM interactions. The phase diagram of this model for J 3 = 0 has been studied previously with DMRG [51]. We remark here already that the numerical results we report are consistent with this previous work in the range of parameters studied in Ref. [51]. Note, however, that the spin patterns in the magnetically ordered phases are changed by the gauge transformation relative to the phases of the SU (2)-invariant model.
Before presenting the numerical results, let us analyze how the spin model coefficients emerge from realistic material parameters. In Fig. 2, we show the contributions of the different orders of the expansion to J 1 , J 2 and J 3 for WSe 2 /WS 2 . It is evident that the third and fourth orders are extremely important to capture the correct physics. Being ferromagnetic for J 1 , the third and fourth orders suppress J 1 turning it even negative for larger twist angles (smaller interactions). In the case of J 2 and J 3 , on the other hand, the spin interactions are boosted by the FIG. 3. CSL region as a function of t2/t1 and t3/t1 for various combinations of U and V from DMRG on an infinite YC8 cylinder as detected by the absolute value of the chiral order parameter Si · (Sj × S k ) averaged over all small nearest-neighbor triangles on the kagome lattice. (a) U/t1 = 75 and V /U = 0.14 corresponding to the phase diagram in Fig. 1(a). (b) For decreased U/t1 = 70 with V /U = 0.14 unchanged, the CSL region moves slightly to the lower right, but narrows. (c) Opposite effect when increasing to U/t1 = 70, still at V /U = 0.14. (d) Changing V /U = 0.13 has a similar effect as decreasing U , (e) larger V /U = 0.15 behaves comparable to increased U .
higher orders. In this way, J 1 , J 2 and J 3 can be of the same order of magnitude permitting the rich phase diagrams tunable with experimental parameters we report below, despite t 2 and t 3 being an order of magnitude smaller than t 1 . Note that this is not a sign of a breakdown of our expansion, but rather comes from the fact that the higher orders include virtual processes that do not involve intermediate states with a double occupancy and whose contribution is therefore not suppressed by factors of 1/U .
Numerical results.-To obtain the ground state phase diagram of our spin model, we perform DMRG simulations on an infinite cylinder of YC8 geometry [29,42]. We map out the phase diagram for fixed U/t 1 = 75 and V /U = 0.14 for varying t 2 /t 1 and t 3 /t 1 , and for fixed t 2 /t 1 ≈ 0.145 and t 3 /t 1 ≈ 0.08 for varying U/t 1 and V /U , shown in Fig. 1(a) and (b), respectively. The interaction strengths of part (a) correspond to the estimates for WSe 2 /WS 2 at θ = 0 • twist angle with dielectric constant ε ≈ 9.5. The same holds for the t ratios of part (b). The derivation of these model parameters from ab initio calculations is detailed in our SM [42]. For fixed interactions in Fig. 1(a), we find two spin liquid phases, namely the CSL [51][52][53][54][55] and the KSL, which is connected to the nearest-neighbor Heisenberg point. For small t 2 and t 3 , we observe a phase in which a ferromagnet and a √ 3× √ 3 state [56][57][58] in the x-y plane are the degenerate ground states due to the gauge transformation [42]. In an SU (2)invariant version of the model, these would be the two degenerate √ 3 × √ 3 states with opposite vector chirality. Next to it, we find a valence bond crystal (VBC) with spontaneous bond order. Above the diagonal, the phase diagram is dominated by the cuboc1 * state, the gauge transformed version of the cuboc1 state, a state with finite scalar chirality [59]. On the bottom right, the ground state is the q = (0, 0) * state, the gauge transformed version of the coplanar q = (0, 0) order [56,58,60]. The relation between the magnetic orders, the gauge transformation and the states of the SU (2)-invariant Hamiltonian is further discussed in the SM [42].
The phase diagram in Fig. 1(b), with the hopping values fixed at our estimates for WSe 2 /WS 2 at θ = 0 • similarly exhibits a finite CSL region in the center. For stronger interactions, the KSL takes over. For smaller U and V , we find the cuboc2 * state, the gauge transformed version of the cuboc2 magnetic order [59], and again a region with degenerate ferromagnetic and √ 3 × √ 3 ground states. Most of these two phases coincide with the area in which the nearest-neighbor spin interaction turns ferromagnetic, in agreement with classical phase diagrams [59].
Chiral spin liquid.-To identify the CSL, we primarily use the chiral order parameter (OP) S i · (S j × S k ) where i, j and k denote the sites around a small triangle in the kagome lattice formed out of nearest-neighbor bonds. The chiral OP for the various values of U and V is depicted in Fig. 3 and clearly indicates the region of the CSL. We observe that the CSL region widens or narrows and shifts with changing interactions. We note that both the cuboc1 * phase as well as the q = (0, 0) * phase can attain a nonzero chirality on the small triangles, but in a staggered pattern such that it averages out over the unit cell.
Since the chiral OP alone is not an unambiguous signature of the CSL, we also compute the momentumresolved entanglement spectrum and the spin Hall conductivity from flux insertion. In the entanglement spectrum, a momentum k y around the cylinder and S z eigenvalue of the corresponding Schmidt state can be assigned to each level. The chiral SU (2) 1 Wess-Zumino-Witten (WZW) conformal field theory describing the edge of the CSL then predicts a certain multiplet structure in each S z sector [61][62][63], which we confirm in Fig 4(a). We show the response of the system when threading spin flux trough the cylinder in Fig. 4 so that a spin up picks up a phase of e iφext when going around the circumference. After φ ext = 2π flux insertion, the expectation value of the spin in the left half of the system increases by S z = 1/2 which implies a quantized spin Hall conductivity of σ spin xy = 1/2 as expected for the Kalmeyer-Laughlin CSL [26].
Kagome spin liquid.-The presumed ground state of the kagome spin model with only nearest-neighbor Heisenberg coupling is also a spin liquid, whose nature remains under debate [28][29][30][31][32][33][64][65][66][67][68][69]. In our t 2 -t 3 phase diagram in Fig. 1(a), we find a small strip of the KSL below the CSL. However, the separation between the KSL the q = (0, 0) * state is subtle to detect. The spins in the q = (0, 0) * can partly point out of the x-y plane which happens in the region of the phase diagram that we ascribe to the q = (0, 0) * phase. The part that we identify as the KSL has S z i = 0. The latter region could also be a weakly ordered q = (0, 0) * state with spins lying in the x-y plane. However, at negative t 3 /t 1 ≈ −0.032 and t 2 /t 1 ≈ 0.04, J 2 and J 3 almost vanish and we obtain a nearly only nearest-neighbor spin model. Since we find no signs of a quantum phase transition between this point and the S z i = 0 region in question, we assign it to the KSL phase and take the line at which a finite S z i develops as the phase boundary. The details of this reasoning are given in the SM [42]. We emphasize that it is not within the scope of this work to give further insight into the nature of the KSL phase, but that we identify a phase with paramagnetic features which is distinct from the q = (0, 0) * phase and adiabatically connected to the ground state of the nearest-neighbor only model. By this identification of phases, the entire upper right region in the V -U phase diagram of Fig. 1(b) falls into the KSL phase as well. FIG. 5. Absolute values of J ratios from our DFT estimates for WSe2/WS2 as a function of twist angle for different values of ε. For ε = 7 and 11, the ratio J2/J1 changes sign due to J1 becoming negative indicated by the blue (positive) and red (negative) shading. At ε = 15, J1 is negative over the entire twist angle range. J2 and J3 are positive everywhere while J 3 is always negative. We chose the absolute values here for better presentation clarity on a log scale. We also include the ratio of V /t1 (red line) and V /t1 = 5 (red dashed line). Above this value, almost the entire particle density is localized on the kagome lattice ensuring the validity of our spin model description.
Experimental realization and detection.-The red dots in our phase diagrams in Fig. 1 mark our estimate for the hopping and interaction values at ε ≈ 9.5 for the first flat valence band in aligned WSe 2 /WS 2 [22,42]. We thus predict that aligned WSe 2 /WS 2 falls just onto the transition line between the CSL and the KSL, suggesting a real material manifestation of these exotic spin states. As in any TMD heterobilayer, the interaction strengths U/t i and V /t i are tunable through engineering the dielectric environment. The ratio V /U can be changed by adjusting the screening length, which can be modified by the distance between the conducting gates and the bilayer. The influence of these two tuning knobs is demonstrated in Fig. 1(b) which shows that the system can be driven deeper into one of the spin liquid phases. In addition to the dielectric environment and gate distance, there are several other tuning parameters. The choice of TMD material influences the effective modelmost notably, compounds with Mo have a larger particle effective mass than We-based TMDs, which then leads to flatter bands and larger effective interactions Similarly, applying uniaxial pressure onto the bilayer increases the interlayer moiré potential, which strengthens interactions as well. We found that these two factors also lead to a slight increase in the t 3 /t 2 ratio. On the other hand, the interaction strength can be reduced by increasing the twist angle.
The values of the resulting spin interactions in WSe 2 /WS 2 as a function of twist angle for different values of ε are shown in Fig. 5. Generally, the magnitudes of the coefficients are distributed as expected with |J 1 | > |J 2 | > |J 3 | > |J 3 |. As discussed above, J 1 turns negative and we obtain a ferromagnetic model for larger ε and/or twist angle while J 2 and J 3 always stay posi-tive and J 3 negative. The relevant energy scale for the spin physics we consider is given by the nearest-neighbor exchange constant J 1 which is rather small due the large length scale in moiré systems. For the value of U/t 1 = 75 in Fig. 1(a), our estimates lead to J 1 ≈ 0.03 meV corresponding to ≈ 350 mK which severely challenges experimental detection of the exotic spin phases. It has been proposed that magnetic order can be diagnosed by the splitting of exciton resonances [70]. Further promising techniques include magnetic resonance force microscopy (MRFM) [71], spin-polarized scanning tunnel microscopy (STM) [72], and nitrogen vacancy (NV) centers [73]. The detection of spin liquids beyond the absence of magnetic order is even more challenging. One possible approach is to use the optical access to the spin degree of freedom in TMDs due to spin-valley locking [74,75] which may allow for the dynamical detection of the time-reversal symmetry breaking or quantized spin Hall conductivity of the CSL. Recently, magneto-optical Faraday rotation was proposed to detect the CSL in the triangular lattice Hubbard model [37,76].
Conclusion.-We have demonstrated that a variety of magnetic phases can be realized in an effective spin model on the kagome lattice which describes TMD bilayers at a filling of 3/4 holes or electrons. In particular, the chiral spin liquid as well as the kagome spin liquid can emerge for experimentally realistic parameters, in addition to several magnetically ordered phases. Moreover, the tunability of TMD moiré systems allows for a systematic search of elusive spin liquid physics, and as such it opens up a promising new direction in the search of highly entangled quantum matter. Apart from our approach, two additional proposals for a kagome charge arrangement in twisted TMD bilayers have recently been put forward whose spin physics has yet to be investigated [77,78].
The data and code used to create the reported results are available at [79].

MODEL PARAMETERS FOR WSe2/WS2
We derive the model for the flat bands based on the continuum model with spin-orbit coupling of Refs. [1,2]. The continuum model of a valence flat band in aligned WSe 2 /WS 2 is described by the following Hamiltonian, where m * = 0.36m e is the effective hole mass, Q is the momentum relative to the K/K point (depending on the spin), and V (r) is the moiré potential. The latter is characterized by the moiré reciprocal lattice vectors ) with a M = 7.98 nm, and V 1 = V moiré e iψ with (V moiré , ψ) = (7.7 meV, 106 • ) determined using ab initio density functional theory calculations [2].
The dependence of the hopping parameters on m * , V moiré and the twist angle is shown in Fig. S1. Increasing effective mass and moiré potential can lead to a higher |t 3 /t 2 | ratio (orange line) which generally favors the CSL (see Figs. 1 and 3 in the main text). The Wannierization provides us also with a shape of the localized Wannier orbitals on the moiré triangular lattice. For our aligned WSe 2 /WS 2 , the Wannier orbital is a Gaussian centered at the W/Se region of the moiré unit cell with width σ = 1.39 nm.
FIG. S1. Behavior of the hopping parameters as a function of effective mass, moiré potential and twist angle. Based on the Wannier orbital size, the onsite and nearest-neighbor repulsion can be calculated using a screened Coulomb potential where d is the screening length. For an infinite screening length d = ∞, the resulting onsite Hubbard and nearest neighbor repulsion are expressed analytically in terms of the Wannier orbital width σ, For ε ≈ 9.5 and d = ∞, the values for aligned WSe 2 /WS 2 are U/t 1 = 75 and V /t 1 = 10.5.

MEAN FIELD THEORY FOR KAGOME CHARGE ORDER
The simplest description of the kagome charge order amounts to a mean field theory for spinless electrons with nearest neighbor hopping on a triangular lattice and nearest-neighbor repulsion V ij n i n j . The mean-field decoupling means replacing n i n j → − n i n j + n i n j + n i n j . The resulting mean-field Hamiltonian is solved self-consistently for fixed particle filling n = 3/4. The T = 0 expectation value for the occupation difference between the kagome sites and the empty site is shown in Fig. S2. We also find within our mean field theory that for V /t > 0.4 a full gap in the spectrum appears, with the charge excitation gap at large V scaling as ∆ ≈ 2V − t.
At V /t = 5 the charge on the occupied sites exceeds n A = 0.98. We choose this as threshold for charge localization, and therefore as a limit on the applicability of our strong coupling expansion of the effective spin model.

EFFECTIVE SPIN MODEL
We expand the extended Hubbard Hamiltonian from Eq. (1) in the main text in the ratio of hoppings to interactions according to Ref. [3]. We keep all terms of second and third order in the hoppings, and all fourth order terms ∝ t 4 1 . We write the XY part of the effective spin model in a more compact notation with S + and S − operators instead of S x and S y . The expressions can be straightforwardly transformed into XX/Y Y and Dzyaloshinskii-Moriya (DM) interactions. To second order, we obtain the usual antiferromagnetic expressions: In the half-filled Hubbard model (one particle per site) without magnetic field, all third order terms vanish since they would break particle-hole symmetry [4]. In addition, they would generate three-spin interactions that would break time-reversal symmetry. In our case, however, there is no particle-hole symmetry and we have empty sites on the triangular lattice whose involvement can create two-spin terms at third order. They are given by Here, the denotes the empty site at the center of a hexagon which forms a t 2 1 t 2 triangle with two nearest neighbors. The • is the empty triangular lattice site inside a hexagon of the kagome lattice (between third nearest neighbors) and the • the filled site between third nearest neighbors on a line of sites. See illustrations in Fig. S3. Note that the terms in the second to last line therefore do only act on third nearest neighbors across a hexagon and the ones in the last line only on third nearest neighbors along a line of sites. The fourth order terms ∝ t 4 1 read Here, i • j denotes summation over third nearest neighbor pairs with empty triangle sites in between and i • j over pairs with a filled site in between, as before. The fourth order terms mostly generate two-site spin operators as well, but also some four-spin terms. These are present on all bond pairs that are connected by a nearest neighbor bond, indicated by jk − lm. Note that they are only generated due to the presence of the nearest-neighbor interaction V and vanish for V = 0. Even for finite V , they give a very small contribution and we neglect them in all the numerical computations of this work. The full spin Hamiltonian we consider is given by where the prime in H 4 denotes H 4 after dropping the four-spin terms.

INFORMATION ON DMRG SIMULATIONS
DMRG calculations are performed on infinitely long cylinders with YC8 and YC12 geometry [5] using the TeNPy package [6]. In the YC geometry, one set of bonds of the kagome lattice is oriented along the circumference as depicted for the YC12 cylinder in Fig. S4. We use S z conservation and keep a matrix product state bond dimension of up to 3200 leading to truncation errors ∼ 10 −5 for the YC8 states. All the phase diagrams in the main text are based on data of YC8 cylinders while the figures characterizing the states in Fig. 4 of the main text and Figs. S6 to S14 in this Supplemental Material come from data on a YC12 cylinder.

GAUGE TRANSFORMATION AND MAGNETIC ORDER OF THE DIFFERENT STATES
As mentioned in the main text, the model with XXZ and DM interactions we study can be transformed into an SU (2) invariant Hamiltonian by a gauge transformation, as illustrated in Fig. S5, which rotates the spins on different sublattices in the xy plane. This transformation is related to the invariance of the spectrum of the underlying Hubbard model under the introduction of 2π flux through a triangular plaquette that has been pointed out in the literature [7,8]. The classical regular orders that can emerge in an O(3) invariant model on the kagome lattice have been classified [9] and it is indeed the gauge transformed quantum versions of those which we identify in our phase diagram. For all phases, we display the spin structure factor (SF) in the extended Brillouin zone with µ = x, y, z. Since we conserve S z and our infinite cylinder is quasi-one-dimensional, we cannot break the remaining U (1) symmetry of the model and S xx (k) = S yy (k). Although we computed the phase diagram on a YC8 cylinder due to computational feasibility, the states we show here were computed on a YC12 cylinder since this geometry contains all the relevant high symmetry points in the Brillouin zone.

Ferromagnet and
Let us start with the simplest order which is the ferromagnet. Its SF and real space correlations are shown in Fig. S6. Since we conserve S z and work in the m z = 0 sector, the order has to develop in the xy plane which clearly manifests itself in a peak at k = 0 in S xx/yy (k) and the uniformly positive xx/yy correlations in real space.
The next order we consider is the coplanar √ 3 × √ 3 state whose SF and real space correlations are depicted in Fig. S7. This state comes in two versions of different vector chirality with (S i × S j ) ·ẑ positive or negative when i and j are adjacent sites in going around a triangle in counterclockwise direction (A and B in Fig. S8). As indicated in Fig. S8, the gauge transformation permutes these states together with the ferromagnet (C). In the SU (2) invariant model, the region for small t 2 and t 3 in Fig. 1(a) of the main text would have the two different √ 3 × √ 3 states A and B as its degenerate ground states. However, since our spins are transformed, the ground states are B and C and the DMRG spontaneously coverges to a √ 3 × √ 3 or an in-plane ferromagnet in that region. With the appropriate initialization, we can reach both states at t 2 = t 3 = 0 as demonstrated in Figs. S6 and S7.  The q = (0, 0) state of the SU (2) invariant model is a coplanar order of spins pointing outward in a 120 • pattern on all upward or downward pointing triangles of the kagome lattice in the orientation it is drawn in Fig. S8. However, the ordering plane does not necessarily have to be the xy plane since any SU (2) rotated version describes the same order. We observe states whose spins are rotated out of the xy plane as demonstrated in the right panel of Fig. S9. Remember that the gauge transformation does not influence the spin value in z direction. We call the gauge transformed version q = (0, 0) * with an asterisk state to indicate the relation to the q = (0, 0) of the SU (2) invariant model. If we denote the gauge transformation by an operator G and an SU (2) rotation by R and we find a magnetically ordered ground state |ψ in our model, then any state |ψ = GRG −1 |ψ is also a ground state. We note that a q = (0, 0) * state with nonzero S z generically acquires a finite chirality since the gauge transformation rotates some of the spins out of the ordering plane. In particular, upward and downward pointing triangles develop opposite chirality which is why the average shown in Fig. 3 of the main text remains zero. In z direction, the peaks of the q = (0, 0) * state stay at the M points of the extended Brillouin zone as in the untransformed q = (0, 0) state, shown in Fig. S9. In the xx/yy SF, however, they are shifted to the K points of the first Brillouin zone by the gauge transformation. We compared the SF found by DMRG to a classical one where we provide the rotation angles of the ordering plane out of the xy plane and find good qualitative agreement.

Cuboc states
The last spin rotation symmetry broken phase in the phase diagram of Fig. 1(a) of the main text is the cuboc1 * phase. In the SU (2) invariant model, the cuboc1 is a state with a 12-site unit cell in which the 12 spins point towards the corners of an eponymous cuboctahedron. Since the gauge transformation has a 9-site unit cell with √ 3 × √ 3 structure, the unit cell of the cuboc1 * generally contains 144 sites and it is not very instructive to illustrate the exact spin orientations. As in the case of the q = (0, 0) * state, we compare the structure factors of Fig. S10 with a classical SF. We therefore start from the cuboc1 pattern as classified in Ref. [9], rotate it and perform the gauge transformation. The rotation angles are given in Fig. S10.
In Fig. S11, we also provide the structure factors and S z for the cuboc2 * state that appears in the phase diagram of Fig. 1(b) in the main text. The same reasoning as for the cuboc1 * applies and we again compare to the SF of a classical rotated and transformed cuboc2 state. We note that the chirality of the small triangles in the kagome lattice of both cuboc * states averages to zero over the unit cell which is why only the CSL displays a finite value in Fig. 3 of the main text. of the cuboc1 * state at U/t1 = 75, V /U = 0.14, t2/t1 = 0.06, t3/t1 = 0.08. The lower two panels show the structure factor of the classical state given in Ref. [9] plus rotations and gauge transformation as in Fig. S9. The rotation angles are φ = 74π 90 and θ = π 2 .

Valence bond crystal (VBC)
The valence bond crystal is a state beyond classical order that is characterized by singlet formation over certain bonds. It does not break any spin rotation symmetries, but the translation symmetry of the lattice and was also reported in the J 1 -J 2 -J 3 kagome model study of Ref. [10]. We show the structure factors and bond spin correlations in Fig. S12. No clear peaks are visible in the SF indicating the absence of spin rotation symmetry breaking. In the bond spin correlations, the pattern of strengths is consistent with the one found in Ref. [10]. Note that the xx correlations are smaller by a factor of 1/2 compared to the zz correlations which is again a consequence of the gauge transformation. The S x i operator of a site at which the spin operators are rotated by 2π/3 expressed in the SU (2) invariantS i operators reads Any pair of sites ij along a nearest neighbor bond has a relative 2π/3 gauge rotation between them so that the generic xx expectation value for such pairs is Since S y iS x j has to be zero in an S z conserving state without spin current and taking into account that the state is SU (2) invariant in terms of theS variables, it follows that

CSL and KSL
We show the structure factors and bond spin correlations of the two spin liquid phases in Figs. S13 and S14. As expected for magnetically disordered phases, there are no sharp features in the SF and no pattern in the bond correlations. The SF of the KSL shows slight bumps at the same values as the q = (0, 0) * state due to the proximity of the two phases in parameter space.
DISTINGUISHING KSL AND q = (0, 0) * As explained in the previous section, the q = (0, 0) * order can come in different orientations in our phase diagram of Fig. 1(a) of the main text. On the YC8 cylinder, we observe two orientations as identified by the S z expectation value on the individual sites. These two are illustrated in the right two panels of Fig. S15. Orientation A (central panel) at t 2 /t 1 = 0.15 shows S z values of a, 0 and −a at sites 1, 2 and 3 of a small triangle whereas orientation B (right panel) at t 2 /t 1 = 0.2 has ( S z 1 , S z 2 , S z 3 ) = (−b, +2b, −b). The KSL state at t 2 /t 1 = 0.1 (left panel) has S z i = 0 everywhere. The transitions between these regions are clearly visible in Fig. S16, where we show the average | S z i |, the average absolute value of chirality over different triangles and the entanglement entropy However, this could still be a state with q = (0, 0) * order and all spins lying in the xy plane. In order to get some more insight into this question, we investigate the coefficients of our spin model for t 3 < 0. For t 3 /t 1 = −0.03246 and t 2 /t 1 = 0.04, both J 2 /J 1 and J 3 /J 1 are smaller than 10 −4 and J 3 < 7 × 10 −4 so that the model at this parameter point can be regarded as the nearest