Time-reversal symmetry breaking Abelian chiral spin liquid in Mott phases of three-component fermions on the triangular lattice

We provide numerical evidence in favor of spontaneous chiral symmetry breaking and the concomitant appearance of an Abelian chiral spin liquid for three-component fermions on the triangular lattice described by an SU(3) symmetric Hubbard model with hopping amplitude $-t$ ($t>0$) and on-site interaction $U$. This chiral phase is stabilized in the Mott phase with one particle per site in the presence of a uniform $\pi$-flux per plaquette, and in the Mott phase with two particles per site without any flux. Our approach relies on effective spin models derived in the strong-coupling limit in powers of $t/U$ for general SU$(N)$ and arbitrary uniform charge flux per plaquette, which are subsequently studied using exact diagonalizations and variational Monte Carlo simulations for $N=3$, as well as exact diagonalizations of the SU($3$) Hubbard model on small clusters. Up to third order in $t/U$, and for the time-reversal symmetric cases (flux $0$ or $\pi$), the low-energy description is given by the $J$-$K$ model with Heisenberg coupling $J$ and real ring exchange $K$. The phase diagram in the full $J$-$K$ parameter range contains, apart from three already known, magnetically long-range ordered phases, two previously unreported phases: i) a lattice nematic phase breaking the lattice rotation symmetry and ii) a spontaneous time-reversal and parity symmetry breaking Abelian chiral spin liquid. For the Hubbard model, an investigation that includes higher-order itinerancy effects supports the presence of a phase transition inside the insulating region, occurring at $(t/U)_{\rm c}\approx 0.07$ [$(U/t)_{\rm c} \approx 13$] between the three-sublattice magnetically ordered phase at small $t/U$ and this Abelian chiral spin liquid.


I. INTRODUCTION
Quantum spin liquid phases are unconventional states of matter that have gained a lot of attention in the last decades due to their fascinating properties and possible future applications in quantum devices like quantum computers [1][2][3]. From a theoretical point of view they are expected to emerge in strongly-correlated systems, for instance in Mott insulating phases. The recent progress in experiments with ultra cold atoms in optical lattices opens the exciting new possibility to simulate a broad variety of such quantum models [4]. The optical lattice allows one to adjust the lattice type as well as the interaction strength, which can be tuned sufficiently to reach the Mott phase [5][6][7]. Furthermore, a synthetic static gauge field (the analog of magnetic flux in electronic systems) can be applied [8,9]. The fundamental degrees of freedom can be adjusted by the choice of atom type. In particular fermionic alkaline earth atoms provide an SU(N ) symmetric spin degree of freedom with N ≤ 10, because of an almost perfect decoupling between nuclear spin and electronic angular momentum [10][11][12][13][14]. This new possibility to realize SU(N ) symmetric Hubbard models within the Mott phase [6,7] creates a strong motivation for further theoretical investigations to identify potential realistic hosts of quantum spin liquid phases.
These systems can quite generally be described by an SU(N ) Hubbard model with a uniform "charge" flux defined by the Hamiltonian where t denotes the hopping amplitude of the fermions with color α = 1, . . . , N . The Peierls phases φ ij are chosen such that the flux per elementary plaquette of the lattice amounts to Φ and U ≥ 0 parametrizes the repulsive interaction strength.
We use the convention t ≥ 0, which is the natural sign for fermions hopping on a lattice with a quadratic spectrum at zero momentum. In the following, we will pay special attention to timereversal symmetric models (Φ = 0 or π) since the main objective of this article is to identify chiral phases that spontaneously break time-reversal symmetry. As we shall see, in the Mott insulating phase with one particle per site, and with t ≥ 0, such a phase is realized when Φ = π, a situation to which we will refer as π-flux case. Now, on the triangular lattice, introducing a π flux is equivalent to changing the sign of the hopping integral. Accordingly, for N = 3 and in the Mott phase with two particles per site, such a phase is realized when Φ = 0 since these two phases are related by a particle-hole transformation that changes the sign of the hopping term.
In the strong-coupling limit U t, and with an average filling of one particle per site, the system is in a Mott insulating phase, and effective magnetic descriptions are given in a basis of SU(N ) spins in the fundamental representation on every site. Up to second-order in t/U the effective description is given by the SU(N ) Heisenberg model with the transposition operator P ij , which exchanges the arXiv:1802.03179v3 [cond-mat.str-el] 2 May 2020 states on sites i and j where α and β run over N different colors [15], i.e. SU(N ) spins in the fundamental representation. Up to second order, the coupling constant is given by J = 2t 2 /U . In a particular large N limit of this model SU(N ) chiral spin liquids (CSL) were first reported [16,17]. In these phases parity and timereversal symmetry are spontaneously broken, but their product is not [16], and most importantly, they exhibit topological order if the spectrum is gapped. This shows up as a 2N -fold degeneracy of the ground state on a torus [17]. By contrast, for SU (2), the Heisenberg model Eq.
(2) realizes a 120 • long-range ordered state [18][19][20][21][22]. However, subleading four-spin interactions arising in order four in t/U are known to trigger a first-order phase transition within the Mott phase to an exotic quantum disordered phase. Several investigations have led to numerical evidence in favor of a gapless spinon Fermi surface phase [23][24][25][26]. However, recent numerical studies [27][28][29][30][31][32][33] argue in favor of the presence of a chiral spin liquid or the existence of a gapped Z 2 spin liquid or even a gapless Dirac spin liquid. So the definite identification of this exotic quantum phase within the N = 2 Mott phase is not yet settled.
For SU(3) symmetric fermions on the triangular lattice, which is the system of interest in this work, it is useful to consider the expansion to third order in t/U . The effective model is then given by the J-K model defined by the Hamiltonian where the first sum runs over all nearest-neighbor sites, and the second sum over all elementary triangles of the lattice. The operator P ijk = P ij P jk is a ring exchange operator that cyclically permutes the states between the sites i, j, and k. The coupling constants are given by J = 2t 2 /U − 12 cos (Φ) t 3 /U 2 and K = −6e iΦ t 3 /U 2 , where Φ is the flux per triangular plaquette. A number of studies were performed for this J-K model already [34][35][36][37]. For purely real ring exchange, besides the conventional three-sublattice magnetically long-range ordered (3-SL LRO) phase [38][39][40], two phases have been predicted by variational Monte Carlo (VMC): a d x + id y CSL with spontaneous time-reversal symmetry breaking followed by a spin nematic phase [34]. However, the uniform π/3-flux and 2π/3-flux CSL states were not considered. A later mean-field study found that those are favored energetically for 0.41 K 6.0 and K 6.0, respectively [35]. This competition between many phases calls for investigations beyond mean-field and VMC. So far, this was only done for the J-K model with purely imaginary ring exchange K, hence explicitly broken time-reversal symmetry, and the presence of robust CSL phases was confirmed for all N from 3 to 9 by exact diagonalizations (EDs) on top of VMC [37]. The first goal of the present paper is to extend this investigation to the case of purely real ring exchange, where we find in particular a previously unnoticed phase that breaks lattice-rotational symmetry, and a spontaneous time-reversal symmetry breaking CSL.
With respect to experiments, the most relevant open question is then whether the CSL phase is also present in the SU(3) Hubbard model, for which the J-K model only provides a reliable description at very strong coupling U . Since EDs of the SU(3) Hubbard model are limited to very small system sizes (12 sites, see Sec. V B), a direct investigation of the Hubbard model is difficult. One way to overcome this difficulty is to push the expansion in t/U to higher order to extend the range of validity of the effective spin model. Indeed, not only is the effective model better because it contains more terms, but one can also test the effect of each additional order in the expansion and get a more precise idea of the range of validity of the effective model. So, the second goal of this paper is to push the expansion in t/U to higher order (we will reach order five), and to identify the physics of the Hubbard model in the range where U/t is not too large so that terms beyond nearest-neighbor exchange play a role, but still large enough to allow one to draw conclusions on the basis of the fifth-order expansion. As we shall see, ED complemented by VMC provides strong evidence in favor of a uniform π/3-flux CSL for moderate K/J, as well as for the Hubbard model in the Mott phase below (U/t) c ≈ 13.
The paper is organized as follows: At first we summarize the key results in Sec. II. In Sec. III, we review the various methods used in this manuscript -EDs (III A), variational Monte Carlo (III B), and the derivation of effective models using degenerate perturbation theory (III C). Section IV is devoted to the J-K model. We derive its phase diagram on the basis of ED and VMC, and provide evidence for a CSL. We then turn to the Hubbard model in Sec. V. We start with a detailed presentation of the fifth-order effective model in Sec. V A. We then benchmark this model in Sec. V B by comparing its properties with those of the π-flux Hubbard model on 12 sites. Finally, we use this effective model to discuss the phase diagram of the π-flux Hubbard model and provide evidence for a CSL with spontaneous symmetry breaking in Sec. V C. The work is wrapped up in Sec. VI including a proposal for the experimental realization of the CSL.

II. KEY RESULTS
Before we give a detailed description of the methods and the full presentation of the results, we summarize our main findings. Let us start with the ground-state phase diagram of the J-K model with real ring exchange K shown in Fig. 1. We use the J-K notation, as well as the angle α in the range −0.3 ≤ α/π ≤ 1, which relates to the previous parametrization by J = cos α, K = sin α [34] . At the Heisenberg point, i.e. α = 0 (K = 0, J = 1), the well known 3-SL LRO phase is present. It is stable to the addition of a non-vanishing ring exchange of either sign. For α −0.25π (K < 0, K/J −1) it eventually gets replaced by a ferromagnetic phase (FM). For K/J > 0 an increase of the ring exchange K triggers a phase transition to a spontaneously time-reversal symmetry break-  1. Comparison of the predicted phase diagram for the triangular lattice J-K model on the 12-, 21-and 27-site clusters from ED (top three rows) and the 144-site cluster from VMC (bottom row). Around the Heisenberg point (α = 0) the 3-sublattice ordered phase (3-SL LRO,green) is present. For increasing values of α up to α = π a π/3-flux chiral spin liquid (CSL,orange), a lattice nematic phase (LN,blue), a 120 • LRO phase (white) and a ferromagnetically ordered phase (FM,gray) occur. The stripe state from VMC is expected to be closely related to the LN phase from ED as discussed at the end of Sec. IV.
where the specific value depends strongly on the symmetry of the cluster under investigation. Then, an unexpected phase with broken lattice-rotational symmetry occurs. This phase has not been reported in a previous VMC study [34] and we call it a lattice nematic (LN) phase in the following [41]. In ED this LN phase extends up to α c ≈ 0.7π [K > 0, (K/J) c ≈ −1.4] for the sizes considered, i.e. one requires negative J to leave the LN phase. In VMC however, the phase transition seems to occur at a smaller value of α c ≈ 0.475π [K > 0, (K/J) c ≈ 13]. We attribute this difference to both the strong finite size effects in ED as well as the possible necessary improvements in the VMC ansatz for the LN phase. For larger α another phase is present. This phase was previously introduced as a 120-nematic or J -nematic phase in Ref. 34. However, in the context of an SU(3) symmetric model we think it is best described by a 120 • color ordered (120 • LRO) phase, as it is essentially an SU(2) Heisenberg antiferromagnet embedded into an SU(3) symmetric system. Finally, the FM occurs for The most interesting aspect of the phase diagram is the presence of the spontaneous time-reversal symmetry broken CSL phase at intermediate coupling ratios K/J. We show that this phase has strong chirality correlations within a manifold of six low-lying eigenstates from the singlet sector on a torus in ED, reflecting the topological ground-state degeneracy. The modular matrices of the apparent VMC states yield Abelian anyonic exchange statistics for four quasi-particles with topological spin ±2π/3 and chiral central charge c = 2, in agreement with predictions for the π/3-flux CSL [16].
In view of the numerous phases discussed so far in the context of the J-K model, and for the benefit of the reader, let us summarize the status of the various phases suggested so far:  (3) Hubbard model with Φ = π on the triangular lattice. In the Mott-insulating phase for small t/U the 3-SL LRO and π/3-CSL are found by ED and VMC. The metalinsulator transition is estimated in Sec. V B as (t/U ) mi c ≈ 1/8.5 ≈ 0.12 (ED on 12 sites for the Hubbard model). The large uncertainty on this value and the nature of the Mott phase in that area is indicated by the grey area.
• Ferromagnetic phase: present in all investigations that looked at this parameter range.
Second, we find that the physics of the J-K model at small and intermediate coupling ratios, including the CSL, are present in the Mott phase with one particle per site of the SU(3) π-flux Hubbard model. To this end, we derive the fifthorder effective spin model for the SU(N ) Hubbard model in the strong-coupling regime t/U → 0 for general flux. For Φ = π corresponding to positive ratios K/J and N = 3, this spin model is qualitatively converged up to ratios t/U including the phase transition point (t/U ) c ≈ 0.075 from ED and (t/U ) c ≈ 0.067 from VMC between the 3-SL LRO and the spontaneous time-reversal symmetry broken CSL phase. The 12 site ED results for the Hubbard model suggest a metalinsulator transition at (t/U ) mi c ≈ 1/8.5 ≈ 0.12. This implies the existence of the spontaneous time-reversal symmetry broken CSL in the Mott phase of the SU(3) Hubbard model with Φ = π on the triangular lattice as illustrated in Fig. 2.
Finally, we note that the CSL is also present in the Mott phase with two particles per site of the SU(3) Hubbard model without a flux.

Spin models
For the ED of spin models we employed a basis within irreducible representations (irreps) of the SU(3) group, which was introduced by two of the present authors [42]. Details are given in Appendix A. The method takes advantage of the full SU(N ) symmetry, and thus reduces the size of the Hilbert space under study much more effectively for N > 2 than the usual approach of U (1) color conservation together with lattice symmetries. A drawback is that the approach is more contrived to obtain common quantum numbers like the lattice linear momentum or the angular momentum, which are useful for identifing the appearing states [43]. To this end, we mainly used the following observables. For magnetically ordered phases, the most prominent observable is the spin-spin correlator S(i, j) between sites i and j. It is given by the expectation value of the transposition operator To find the ordering momentum, we studied the spin structure factor where N s is the number of sites and the last term gives the auto-correlator. In addition, we also used the dimer correlator, which is the connected correlator of two bonds. In terms of transposition operators, it reads The identification of the CSL phase is supported by the real space connected scalar chirality correlator C( i, j, k , l, m, n ) for two individual triangles with sites i, j, k and l, m, n oriented in the same way. The scalar chirality for a single triangle is defined as Then, the connected scalar chirality correlator is given by With the overall chirality or chirality signal we refer to the sum of all connected chirality correlators where i, j, k is a fixed reference triangle and l, m, n are all triangles which share no site with the reference triangle and are labeled in the same orientation. For a CSL, which features long-range chiral correlations, the chirality signal should be strong and the real space connected chirality correlator should be positive and uniform across the lattice. The calculations were performed on the 12-site, 21-site, and 27-site clusters shown in Fig. 19.

Hubbard Model
The full Hilbert space of the SU(3) Hubbard model contains states with multiple occupied sites and we cannot employ the same basis as before. Instead, we use an ordinary U(1) quantum number basis, where we label sectors by the numbers of fermions of each of the three colors a, b, c as (N a , N b , N c ).
We also attempt to pin down the metal-insulator transition point by computing the partice-hole charge gap. For a general SU(N ) Hubbard model this can be done by calculating the lowest energies of three different Hilbert space sectors, namely the sector with the same number of particles per color at 1/N filling, the sector with one additional fermion, and the sector with one fermion less for one color compared to the nominal filling. Then, the particle-hole charge gap is where the numbers in parenthesis label the deviation from the Mott filling. In order to obtain an estimate for the metalinsulator transition, one has to extrapolate the large-U linear part of the particle-hole charge gap. The point where the extrapolated line crosses zero, so where the charge gap closes, is an estimate for the location of the metal-insulator transition.

B. Variational Monte Carlo
The VMC approach is based on projected 1/3 filled noninteracting parton wave functions. We start from a nearestneighbor tight-binding model where α denotes the color, and f † i,α (f i,α ) are creation (annihilation) operators of an α-fermion (parton) at site i. For each color the one-fermion states are filled up to one-third filling. Then, a Gutzwiller projection [44][45][46] excludes the charge fluctuations. The energies, symmetries, and overlaps of these variational states can be evaluated using a Monte Carlo sampling. The hopping amplitudes and the on-site fields of the tight-binding model serve as variational parameters. Similar studies of the J-K model for real K were carried out by Bieri et al. [34]. They only considered time-reversal symmetry breaking through complex pairing terms (p + ip and d + id states in their paper) and did not consider any state with complex hoppings and non-time-reversal invariant flux configurations. The latter is a natural way to create CSL variational states [16,36,37]. In the following, we provide a list of the states we considered. We carried out these simulations on two clusters of 6 × 6 = 36 and 12 × 12 = 144 sites, both compatible with all the listed scenarios. A detailed description of our results is given in Sec. IV B, where we also compare the energies of the most competitive states to the ones proposed in Ref. 34.
CSL states-To create spin liquid variational states, we set all the h i,α on-site terms to zero, and the t α i,j hopping amplitudes to unity. The phases of the hopping amplitudes are such that the flux on every triangular plaquette is the same. We considered the cases with mπ/6-flux (m ∈ Z) per plaquette. We also made calculations on states with opposite fluxes on up and down triangles with fluxes ±mπ/3.
Plaquette ordered states-By increasing the magnitude of the hopping amplitudes around distinct plaquettes in Eq. (11), we can create variational states with lower bond energies around these plaquettes. In particular, we focused on the possibility of strong bonds around triangular plaquettes corresponding to the formation of SU(3) singlets. We considered all the possible plaquette coverings on the 36-site cluster, with a uniform zero-or π-flux per plaquette, changing the magnitude of hoppings around the plaquettes from 0.5 to 2 (the magnitude of the other hopping amplitudes were set to 1). We examined the all-up plaquette covering [shown in Fig. 3(a)] in more detail for both the 36-and the 144-site system. On top of changing the magnitude of the bonds around the triangular plaquettes, we considered various flux configurations with either zero-or π-flux per plaquette for each of the three different types of triangles [marked with φ 1 , φ 2 , and φ 3 in Fig. 3(a)]. We also considered the combination of the all-up plaquette order with a uniform mπ/6-flux (m ∈ Z) per plaquette.
Stripe order-Similarly to the plaquette ordering, by strengthening the hoppings along chains in Eq. (11), we generate variational states of weakly coupled SU(3) chains. We considered the cases of straight and zigzag chains [c.f. Fig. 3(b) and (c)] with 0 ≤ t inter-chain /t intra-chain ≤ 2 on the 36-and 144-site clusters. Note, that the zigzag chains are compatible with the 36-and 144-site clusters, but not with the 21-and 27-site clusters used for ED. On top of varying t inter-chain /t intra-chain , we considered uniform fluxconfigurations of 0-, π/3-, π/2-, 2π/3-, and π-flux, as well as opposite fluxes on up and down triangles with the same values. The straight stripe π-flux state is related to the LN phase from ED in Sec. IV A.
Color order-If the chemical potentials are different for different colors, the resulting variational states will have SU(N ) symmetry breaking color order. We considered states with a three-sublattice color order as shown in Fig. 3(a). On each sublattice the on-site term for only one color is non-zero, i.e. h i,α = h if i is in sublattice α. The magnitude of the onsite field is varied between 0 < h < 2. We also combined this color order with the all-up plaquette order and the straight stripe order of the previous points.
Gutzwiller projected states can also be used to access different topological sectors by introducing twisted boundary conditions in Eq. (11) before Gutzwiller projection [37]. Then, if the number of non-zero eigenvalues of the overlap matrix is independent of how many different twisted boundary conditions are considered, it gives the number of linearly independent states. The phases φ1, φ2, and φ3 were either chosen to be uniform with values mπ/6 (m ∈ Z), or set to a combination of 0 or π. For the 3-SL LRO (shown by the color of the vertices) the on-site chemical potential was changed from 0 to 2. For the (b) stripe order, and (c) zigzag stripe order, the strengths of the hopping on the purple bonds are set to 1, and the strengths on the thin black bonds between the stripes are tuned between 0 and 2. In both stripe cases we considered either uniform flux configurations or opposite fluxes on up and down triangles with φ = {0, π/3, π/2, 2π/3, π}.

C. Derivation of effective models
The effective model for the Mott phase of a general SU(N ) Hubbard model with an arbitrary flux for an average filling of one particle per site in the strong-coupling limit is stated and analyzed in Sec. V A. Here, we briefly motivate the underlying calculations, which are based on a linked-cluster expansion (LCE) and degenerate perturbation theory. A similar approach was successfully applied to the SU(2) Hubbard model on the triangular lattice in Ref. 25.
A LCE is a technique to extend results on small finite clusters to larger clusters or to the thermodynamic limit. We use a very appealing LCE approach along the lines of a white-graph expansion [47], which we employ to simplify the subtraction process, as well as to include complex phases. This is explained in detail in Appendix B.
On every linked cluster, we use degenerate perturbation theory [25,[48][49][50][51] about the strong-coupling limit t → 0 in Eq. (1), and determine effective descriptions in the subspace of states with exactly one fermion per site. A link is created by the perturbative hopping of a fermion between two sites. We find that in order k a linked cluster only yields a nonvanishing contribution if the number of links that are part of a loop plus twice the number of links that are not part of a loop is smaller or equal to the order. In order 5 on the triangular lattice, six linked-clusters give a non-zero contribution, namely the dimer, the trimer, the triangle, a triangle with one additional site, a four-site loop, and a five-site loop. For each cluster one derives the associated reduced effective Hamiltonian, which can be written in a compact form (12) where the coupling constants depend on t/U and Φ. Here and in the following, all effective Hamiltonians and effective couplings are given in units of U . Every linked cluster yields a [4,4,4] 12 sites, [6,6,0] 21 sites, [7,7,7] 21 sites, [11,10,0] 27 sites, [9,9,9] 27 sites, [14,13, Ground-state energies from ED on clusters with Ns = {12, 21, 27} sites for the J-K model with the coupling constants J = cos α and K = sin α. Different point shapes indicate different symmetry sectors. The singlet sector [Ns/3, Ns/3, Ns/3] yields the ground states of the 3-SL LRO, CSL and LN phase below α/π ≈ 0.7. For 0.7 α/π 0.85 the 120 • color order state, which contains only two colors from the sector [Ns/2, Ns/2, 0] or [(Ns + 1)/2, (Ns − 1)/2, 0], is present. In the FM phase for α/π 0.85, where the energy per site is identical for all clusters, only a single color is present and therefore the ground states lie in the symmetry sectors [Ns, 0, 0]. Wherever the phase boundaries are not the same for different lattice sites Ns, the color of the boundary indicates the corresponding cluster size. different set of exchanges that are embedded on the full system in the second step. Some of these exchanges are unique to the specific cluster, whereas other interactions emerge from contributions of a variety of linked-clusters. We perform these calculations for all interactions in the effective model in order 5 for the thermodynamic limit (Sec. V A) and in order 4 for the 12-site cluster with PBCs.
To improve and assess the quality of convergence of the coupling constants, we employ Padé extrapolants [52]. These are rational functions and can be denoted with [m, n] where the degree of the numerator is m and that of the denominator is n. At certain parameters, unphysical divergences emerge, which have to be excluded.

IV. J-K MODEL
In this section, we determine the phase diagram of the J-K spin model for arbitrary real values of the J and K parameters by ED and VMC. The results of both methods are compared at the end of the section.

A. Exact diagonalization
Using ED, we studied the J-K model on the 12-, 21-, and 27-site triangular clusters with periodic boundary conditions (PBCs). We start with a discussion of the ground-state energy per site plotted as a function of the parameter α in Fig. 4. We show the energies for three different SU(3) irreps: the singlet, the ferromagnet and the irrep which corresponds to the lowest effective SU(2) sector (for reasons which become clear below). We discuss the evidence leading to the labeling of the different phases in the following. Note that the CSL seems to be squeezed at the interface between two quite extended phases, similar to scenarios in SU(2) quantum magnets, where spontaneous time-reversal symmetry breaking CSL have been observed [53][54][55].
The appearance of the CSL when increasing K > 0 (J > 0), or α > 0, can be understood with the low-energy spectrum of the 21-site cluster from ED in Fig. 5, where the marker size illustrates the overall chirality signal for the lowest three singlet energy levels (whose degeneracies are not resolved here). For small values of K/J the first excited state above the singlet ground state is in the adjoint irrep [8,7,6], which corresponds to the tower of states expected for the 3-SL LRO phase [56]. With an increasing ratio K/J two singlets cross the adjoint irrep around (K/J) c ≈ 0.31 (α c ≈ 0.096π) and become the lowest excitations. In the same parameter range the total chirality of the ground state, as well as that of the two levels from the singlet sector above, increases, which is an indication for a chiral phase. More importantly, the degeneracies of the three lowest singlet levels have been determined to be 1-4-1 corresponding to a total of six states, which may form a six-fold degenerate ground state in the thermodynamic limit. Such a degeneracy is associated with spontaneous chiral symmetry breaking in the SU(3)-symmetric π/3-CSL on a torus [17]. However, these states are not very well separated from higher excited states on the 21-site cluster. The signature is more pronounced than on the 12-site cluster though. On the 27-site cluster the splitting between the six low-energy states and the states above is comparable to that on the 21-site cluster.
In the upper panel of Fig. 6 the magnetic structure factor at the ordering momentum of the 3-SL LRO phase, the K  [8,7,6] singlet states 1-3 [7,7,7] other singlet states [7,7,7] FIG. 5. Spectrum of the J-K model on the 21-site cluster from ED. The marker size corresponds to the overall chirality signal and is plotted for the lowest three states (the degeneracy of these states is 1, 4, and 1). The three (six when we count the degeneracies) lowlying singlets with strong chiral signal indicate the presence of a CSL phase in an extended parameter space at intermediate values of K/J. relator for the first three singlets, as expected for a CSL phase. The third excited state shows no pronounced chiral order anymore. The signature of the CSL phase varies in its extension in parameter space for different cluster sizes. It is larger for the 21-site cluster than for the 12-and 27-site clusters. This is due to the nature of the adjacent phase at larger values of K/J (α), which seems to have its peak in the magnetic structure factor close to the X point. While the 12 sites cluster has this particular point, the 21 and 27 sites cluster do not, which may lead to an increased parameter space for the CSL phase on those clusters. The phase next to the CSL is a previously unreported spatial symmetry breaking phase which is characterized by strong bonds along one lattice direction and weak bonds along the other lattice directions. It therefore breaks the rotational symmetry of the lattice, leaving a strong signature in the dimerdimer correlations, as can be seen in Fig. 8. Given the limitations in reachable cluster sizes using ED, it remains an open question whether there is any sort of magnetic ordering in this LN phase or not.
Increasing K/J (α) further leads to an SU(2)-like behavior with one color less. We confirm this by finding an SU (2) tower of states at large K/J where the ground state is in the irrep [N s /2, N s /2] (or [(N s − 1)/2 + 1, (N s − 1)/2] if the number of lattice sites N s is odd) of SU(3) as shown in Fig. 9. According to Eq. (13), the ground-state energies in this region can be compared to ED results on the nearest-neighbor spin-1/2 triangular lattice [57], and we find a perfect match.

B. Variational Monte Carlo
As mentioned above, we found three relevant phases in our VMC calculations: a three-sublattice color-ordered state with 0-flux per plaquette, the π/3-flux per plaquette CSL state, and a π-flux striped phase. We plot the energies of these phases together with the states proposed in Ref. [34] in Fig. 10. Due to the different focus of their paper, Bieri et al. use a nomenclature more suitable for spin-one systems. However, we will reinterpret their findings in the SU(3) language.
In the coupling regime −0.248π < α < 0.064π (J > 0, −0.99 < K/J < 0.20), we find that the energy of the 120 • AFM order proposed by Bieri et al. [34] and the energy of the 3-SL LRO we propose are almost exactly matched. In fact, for special values of the parameters used by Bieri et al. their variational states are equivalent to our 3-SL LRO variational states. Namely, for sin η = 2/3 in Eq. (15) of their paper, the three d vectors become mutually orthogonal to each other, and their mean-field ansatz can be transformed to ours with an SU(3) rotation. We find that the 3-SL LRO is the more appropriate identification of this phase, as the 120 • order is just a special case, and in fact not all members of the ground-state manifold have dipole ordering due to the SU(3) symmetry of the model, which can mix the dipolar and quadrupolar moments of the spins 1.   For 0.064π < α < 0.195π (K > 0, 0.20 < K/J < 0.70) we find that a π/3-flux CSL phase is selected. This phase was not considered by Bieri et al. [34], but was predicted by the mean-field study of Lai [35]. We find 6 linearly independent states with similar energies by considering states with ±π/3 flux per triangle (3 states each) and using twisted boundary conditions (compare Sec. III B). This fits the ground-state degeneracy of the CSL breaking time-reversal symmetry spontaneously on a torus with genus g = 1. We studied further topological properties by calculating the modular matrices and find a qualitative agreement for the predicted chiral central charge and a quantitative agreement for the anyonic exchange statistics. Details are given in Appendix D. Beyond the π/3-flux CSL phase, between 0.195π < α < 0.475π, a striped state with stronger bonds along straight chains with t inter-chain /t intra-chain ≈ 0.2 and a uniform π-flux has the lowest energy, superseeding both the d+id phase proposed by Bieri et al. and the 2π/3-flux CSL proposed by Lai [35].
For 0.475π < α < 0.85π (K > 0, K/J > 12.7, and K/J < −0.51), the 120-nematic or J -nematic phase proposed by Bieri et al. is clearly dominant; its energy is not matched by any of the states we considered. In their construction the on-site terms before the projection select states with directors in an 'umbrella' configuration (see Eq. (16) in Ref. 34), interpolating between a ferroquadrupolar and a 120 • quadrupolar order. However, due to the SU(3) symmetry of the model there are many other degenerate states, containing not only nematic states, but other states obtained by global SU(3) rotations. Therefore, we find that a 120 • LRO state is more suitable in this context. In the SU(3) language, a 120 • fully color ordered state is only built out of two colors on every site, while the third color is missing in the system. Interestingly, in this subspace the J-K Hamiltonian simplifies since the nearest-neighbor exchange and the ring exchange terms are not independent. Namely, for three sites around a triangle P ij + P jk + P ki = P ijk + P −1 ijk + 1 (13) in the two-color subspace, just as in the spin-1/2 case. As a result, at the special point K = −J/2, the Hamiltonian for the two-color states is a constant, giving the same energy for any state in this subspace. This macroscopic degeneracy can also be seen in ED results, and the transition between the 120 • LRO and ferromagnetic SU (3)  The phase diagrams for the J-K model by ED and VMC, depicted in Fig. 1, are in good qualitative agreement. For small ring exchanges K around α = 0 we find the 3-SL LRO phase, which is expected in a slightly smaller area for K > 0 in VMC than in ED. This is most likely linked to the π/3-flux CSL at larger ring exchanges K, which is particularly well described as a variational state. The CSL states from VMC and ED show a direct correspondence in terms of symmetries. Details are given in Appendix E. In the VMC the next competitive state is the striped state, which is linked to the LN phase in ED, since both break rotational symmetry. We therefore assume that these phases correspond to the same ground state in the thermodynamic limit, and stick with the term LN in the following. The extension of the 120 • LRO phase is strongly dependent on the method, whereas the one of the FM is quite similar from ED and VMC.

V. HUBBARD MODEL
We found that the J-K model hosts a spontaneous timereversal symmetry breaking CSL for intermediate values of Here we turn to the question whether this CSL is also present in the Mott phase of the corresponding Hubbard model (1). In the limit of small values of t/U the Mott phase is well described by the leading nearest-neighbor Heisenberg model stabilizing the 3-SL LRO. With increasing t/U , the next subleading interaction, the third-order K-term, triggers a phase transition to the CSL. The ED critical value is (K/J) c ≈ 0.31 (α c ≈ 0.096π), corresponding to (t/U ) c ≈ 0.064 [(U/t) c ≈ 15.6] taking the bare third-order series. In the following we show that higherorder contributions do not prevent the occurrence of the CSL phase, and that the critical value remains similar.

A. Effective description
In this section we state and discuss the effective model for the Mott phase of a general SU(N ) Hubbard model at a commensurate 1/N filling for arbitrary uniform fluxes Φ for the thermodynamic limit. The model of the 12-site cluster with PBCs, which we use for the direct comparison between the effective description and the Hubbard model in Sec. V B, is given in Appendix F.
In fifth-order perturbation theory the effective model in the thermodynamic limit contains 13 different types of interactions involving permutations on up to five sites. The effective Hamiltonian reads where the pictogram underneath every sum illustrates which sites on the full lattice are addressed. This has to be understood as follows. The relative angles between the bonds of a graph are fixed and characterize the graph (e.g. the graph of L 3sp d can not be transformed into the one of L 3sp s ), however every possibility of rotation has to be included (e.g. the graph of L 3sp d can be rotated around the axis defined by i and j by π). Then, every distinct set of sites contributes to the Hamiltonian. The prefactors starting in second and third order are We note that the complex phases of ring exchanges in subleading orders of the effective interaction are not identical to the flux through the apparent plaquettes in the Hubbard model. For example, the complex phase of a cyclic permutation on a triangle only equals the flux through a triangle in order 3. The representation of the effective model given above is not unique for real exchange constants if N < 5. This is due to the fact that consecutive permutations acting on more sites than spin colors in the model can be rewritten in terms of sums of various other permutation operators. The best known example occurs for the three-site permutation acting on SU(2) spins, which can be expressed by two-site permutations plus a constant [compare Eq. (13)]. The situation becomes much richer if one considers permutations between more spins. For instance, in the case of SU(3) spins, the four-site permutation can be rewritten in terms of two-spin, three-spin and various four-spin interactions  14) is the more natural one in terms of perturbation theory. If and how a systematic reduction of higher-order interactions to only already included exchanges is possible for orders larger than 5 remains an open question.
In the following, we discuss the most interesting case Φ = π, where the convergence of the series works particularly well. Additional results for other fluxes are presented in Appendix G. The dependence of the largest coupling constants of every order as bare series in t/U is illustrated in Fig. 11. The most important subleading corrections to the nearest-neighbor Heisenberg term come from ring exchanges around triangles, squares, and 5-site trapezoids. The shaded background keeps track of where the CSL is observed in only VMC (light yellow) and in ED and VMC (darker yellow) (compare Sec. V C). The Padé-extrapolants for the dominant nearest-neighbor coupling J and the subleading three-site ring exchange K are shown in the insets of Fig. 11 and Fig. 12, respectively. The relevant ratio K/J is plotted in Fig. 12, where we take the ratio of the extrapolations of J and K. In respect to unphysical divergences, we found a composition of [3,2]-Padé extrapolations for 0 , J, and K, [2,1]-Padé extrapolations for L 2sp s , L 2sp d , L 3sp s , L 3sp d , and L 4sp r , and all other couplings as bare series to work best for flux Φ = π.

B. Comparison
In order to confirm the validity of the effective description, we compare the spectrum of the SU(3) Hubbard model (1) with the spectra of the J-K model (4) and the O(4) (order 4) effective spin model (F2) on the 12-site cluster (compare Appendix F) from ED. The results are shown in Fig. 13. While for the spin models the numbers in the square brackets label a certain irrep of the SU(3) group, the three numbers in the round brackets for the Hubbard model label the number of particles of a certain color, and are thus U(1) quantum numbers. For U ≈ 30t the spectra of both effective models are in reasonable agreement with the spectrum of the Hubbard model. The first excitation in the spin models is in the adjoint representation [N s /3 + 1, N s /3, N s /3 − 1] followed by three singlet levels. For decreasing couplings U , the higher lying singlets start to cross each other in the spin models, as the corresponding excited states do in the Hubbard model. In the effective models at U ≈ 15 − 20t, the first excited singlet crosses the low-energy state of the adjoint representation, which for even smaller values of U is also crossed by two more singlets. One can observe a similar behavior in the spectrum of the Hubbard model, even though the order in which the crossings occur is not exactly the same. In general, the level crossings of the O(4) model seem to match those of the Hubbard model slightly better than those of the J-K model, as expected.
To further check the effective description, we compare the differences between the ground-state energies of the spin models with those of the Hubbard model as can be seen in Fig. 14. The errors of the ground-state energies decay with one order higher than the corresponding order of the effective description, as expected for a valid perturbative approach. Overall, there is a qualitative agreement between the effec- tive description and the Hubbard model in the strong-coupling Mott regime. An important aspect of the problem is that the effective description breaks down at the metal-insulator transition. Therefore, the physics found in the effective spin model is only valid in the Mott phase and it is important to estimate the metal-insulator transition point. To this end, we compute the particle-hole charge gap for the π-flux SU(3) Hubbard model on the 12-site cluster. The results are shown in Fig. 15. They indicate that the metal-insulator transition is located at (U/t) mi c ≈ 8.5. This value is further discussed at the end of the next Subsection, which mainly focuses on the presence of the spontaneous time-reversal symmetry broken CSL in the fifth-order effective model within the Mott phase.

C. Chiral spin liquid
Next, we analyze the effective spin model in Eq. (14) for the specific case N = 3 with π-flux on the triangular lattice. For small values of t/U the couplings up to third order are dominant, hence the J-K model is well converged. The phase transition between the 3-SL LRO and the CSL occurs at (K/J) c ≈ 0.31, which translates to (t/U ) c ≈ 0.09 [(U/t) c ≈ 11] in bare fifth-order. Interestingly, this critical value changes only slightly to (t/U ) c ≈ 0.1 [(U/t) c ≈ 10] when applying Padé extrapolations to the couplings J and K. For both couplings several Padé extrapolants give essentially the same result in this t/U -regime so that the extrapolations work well for the most important interactions of the effective model (see also the inset in Fig. 11). The small difference between the critical ratios from bare series and extrapolation within the J-K model indicates that even the bare series is almost converged up to these t/U values for the couplings J and K.
Since the impact of all the other smaller terms is not obvious a priori, we explicitly study them with the set of Padé   extrapolants for the coupling constants discussed in Sec. V A using ED and VMC. The low-energy spectra of the fourth-and fifth-order effective model for flux Φ = π on the 21-site cluster from ED are illustrated in Fig. 16, where again the point size for the lowest three singlets illustrates the chirality signal.
For large values of U/t the tower of states characteristic of the 3-SL LRO phase appears in the spectrum. This is in agreement with the large structure factor and its extensive scaling at the K point, given in the upper panel of  from the symmetries discussed in Appendix E. Therefore, the CSL is most plausibly present here as well. In the same parameter range U/t ≈ 13 the chirality signal of the ground state increases, as can be seen in the bottom panel of Fig. 17, whereas the structure factor at the K point decreases. This behavior agrees perfectly with the phase transition from the 3-SL LRO phase to the CSL. Within the area of the potential CSL the manifold of six low-lying states is not very well separated from the excited states, but the indications are stronger on the 21-site cluster than on the 12-site one (not shown). Since ED including all exchanges is a lot more costly than for the J-K model, we did not study the 27-site cluster. However, with VMC other singlet states [7,7,7], O(4) singlet states 1-3 [7,7,7], O (5) other singlet states [7,7,7], O(5)  The CSL states behave very similar to the lowest eigenstates from ED regarding energy spectra and symmetries. All these findings strongly point to the realization of a π/3-flux CSL phase with spontaneous breaking of time-reversal symmetry.
The energy spectra of the fourth-and fifth-order model from ED in Fig. 16 behave fairly similar for U/t 10. For large ratios U/t the eigenenergies approach each other, as expected. With increasing perturbation (decreasing U/t) the differences increase. Nevertheless, the same manifolds of six low-lying states emerge. Also on the level of observables the signature of the CSL is present in both models as shown in Fig. 17. So, even though the effective model is not quantitatively converged in the relevant area 11 U/t 13, the signature of the CSL occurs in the third-, fourth-, and fifthorder model, which implies, that it is a definite feature of the effective description of the Hubbard model in this U/t regime. For coupling ratios U/t 10 the behavior of the eigenstates is rather different in the fourth-and fifth-order model, so no statements about this parameter range can be made, apart from the fact that the Mott phase breaks down eventually.
The question is whether the potential CSL phase is within the Mott phase of the SU(3) Hubbard model with flux Φ = π. In the ED calculations on 12 sites, the charge gap (given in Fig. 15) shows a linear behavior above U ≈ 10|t|, which together with the effective model (compare Appendix F) provides the estimate (U/t) mi c ≈ 8.5. Can we interpret this result in first-order perturbation theory around the limit of strong coupling at filling 1/3? Assuming linear behavior for the charge gap the transition point yields ∆ charge ≈ U − 8.5t, and the first order term must arise from the kinetics of charge excitations. These are given by doublon-hole pairs. In the case of the SU(2) model, both the hole and the doubly occupied site are featureless objects, in fact SU(2) singlets, and they behave similarly. In the case of the SU(3) Hubbard model, the doubly occupied site forms the three-dimensional anti-symmetrical irreducible representation, and the motion of the doubly occupied site is more complicated than the motion of the hole. The energy −8.5t originates from the hoppings of the hole and the doubly occupied sites, the contribution from the hole is E 0 (−1) − E 0 (0) ≈ −0.02U − 3.30t, and the contribution of the doubly occupied site is E 0 (+1) − E 0 (0) ≈ 1.01U − 5.10t [c.f. Eq. (10)]. We can see, that the kinetic energy of the doubly occupied site is larger than that of the hole.
The metal-insulator transition from ED on the 12-site cluster at (U/t) mi c ≈ 8.5 [(t/U ) mi c ≈ 0.12] occurs for weaker coupling strengths U than the phase transition towards the CSL. For the J-K model the estimated transition point between 3-SL LRO and CSL is located at slightly larger values (U/t) c ≈ 10, and therefore lies within the crudely estimated extension of the Mott phase. However, the fifth-order effective model includes a larger variety of quantum fluctuations and is therefore more reliable with the CSL below (U/t) 13. We thus expect that the Mott phase of the SU(3) π-flux Hubbard model on the triangular lattice realizes, besides the 3-SL LRO phase, a spontaneous time-reversal symmetry broken CSL phase before the Mott insulating phase breaks down. The full phase diagram for the Hubbard model is illustrated in Fig. 2.

VI. CONCLUSION AND OUTLOOK
To summarize, we have provided very strong numerical evidence that the stabilisation of a chiral phase can be achieved with cold atoms in a simple and realistic setting: three flavours (e.g. 87 Sr or 173 Yb) and a triangular optical lattice. To achieve this, we have first studied the SU(3) J-K model on the triangular lattice and shown that it has a CSL phase that spontaneously breaks time-reversal symmetry for a positive but moderate value of three-site exchange term K. In the search for realistic models with a CSL phase that could be implemented experimentally, this is an important step forward because the two interactions that appear in this model, the nearest-neighbor exchange and the three-site permutation around triangles, are the most relevant ones in the Mott phases of the model when coming from strong coupling. So in principle it is accessible to cold atom experiments on SU(3) fermions by just monitoring the ratio of the repulsion between fermions U to their hopping amplitude.
To make this proposal more concrete, we have discussed in great details the connection between the Hubbard model and the J-K model with positive K. For the physically most relevant sign of the hopping, −t with t > 0 to get a band with a quadratic spectrum around its bottom at zero momentum, a positive three-site term K can be reached in two situations: either in the Mott phase with one particle per site if one introduces a π flux per plaquette to effectively change the sign of the hopping, or in the Mott phase with two particles per site without any flux since it is equivalent to the Mott phase with one hole per site starting from the full system, and the sign of the hopping changes under a particle-hole transformation.
It is also important to check if, upon reducing the repulsion U , the three-site term gets large enough to induce the transition from the three-sublattice color order into the CSL before the Mott transition into a metallic phase takes place. To get an estimate of this Mott transition, we have investigated the Hubbard model directly on small clusters, and indeed the three-site interaction reaches values large enough to be in the CSL phase before the system turns metallic.
Finally, to lend further support to this proposal, we have pushed the strong-coupling expansion for the Mott phase to higher order to check the effect of residual interactions, and the conclusion is that a CSL phase appears to survive the inclusion of these terms.
Altogether, the proposal that there is a CSL phase in the SU(3) Hubbard model on the triangular lattice between the three-sublattice color ordered phase and the metallic phase is we believe fairly solid, and we hope that the present paper will encourage the experimental investigation of this model with cold fermionic atoms.
In that respect, we would like to briefly discuss the temperature effects. Since the chiral phase reported here spontaneously breaks time-reversal symmetry, it is expected to give rise to an Ising transition at finite temperature. It would be nice to have an estimate of this Ising temperature to see if it is accessible for the lowest entropies per particle that can be reached with state-of-the-art cooling protocoles. Unfortunately our methods are limited to zero temperature, and this has to be left for future investigation.
Finally, as a byproduct, we have also come up with a new version of the full phase diagram of the J-K model with real K, including all signs of J and K, and we have shown that, just above the CSL, it contains a large lattice nematic phase that had been overlooked so far. Whether such a phase is also present in the Hubbard model before the Mott transition remains to be seen.
The multiplicity L(c) is the number of ways in which a cluster can be embedded on the full lattice of interest and determines the impact of each cluster. The weight of a cluster W P (c) to the property P is given by the inclusion-exclusion principle. This means that only reduced contributions appearing on the considered cluster but not on its subclusters are taken into account One calculates the physical property P (c) on cluster c and subtracts the weights of all subclusters s. Evidently, all processes on a disconnected cluster are already included in the sum of its pieces and the weight vanishes. Therefore, only linked clusters contribute and the weight of a cluster contains only properties that arise from all bonds of the cluster. This fact can be used to calculate the weights directly with a whitegraph expansion [47], which makes the subtraction in Eq. (B2) unnecessary. To this end, every bond on a linked-cluster is labeled with a different exchange constant during the calculation. For instance, on the triangle we take three exchange constants h 1 , h 2 , and h 3 connecting different sites. The per-turbation is written as The subtraction of the effective Hamiltonian derived in this way is then achieved by taking only terms that include every exchange constant at least once and hence emerge from perturbations that link the whole cluster. This procedure can be extended to include complex phase factors by splitting up a hopping process on a link into the hopping from left to right and from right to left. For the triangle one can choose Appendix C: Eigenvalues of the Quadratic Casimir Operator for SU(2) and SU (3) In Fig. 9 it seems that every eigenvalue of the SU(2) quadratic Casimir operator is shifted to larger values by some offset, but reappears in the SU(3) observable. Since this is not obvious, we determine the origin of this effect in the following. To this end, we first show that the lowest eigenvalue of the SU(2)-like irreps inside SU(3) (the irreps with only two colors) depends only on the number of lattice sites N s . Second, we prove that the eigenvalues have the same slope as their SU(2) counterparts. We start from the SU(3) case, and consider an irrep [n a , n b , n c ], which can also be labeled by only two numbers In the language of Young tableaus p and q give the difference in the numbers of columns between the first and second row, and the second and third row, respectively. The eigenvalues of the quadratic Casimir operator in the SU(3) case are then given by For effective SU(2) irreps inside SU(3) it is n c = 0, which yields We can thus write the eigenvalues of the SU(3) quadratic Casimir operator for two-color states as (n a ) = n a + n 2 a − n a · N s + N 2 which only depends on N s and n a . We then compute the difference of the eigenvalues for increasing n a From Eq. (C6) it follows that the eigenvalues grow with n a . The smallest eigenvalue of the effective SU(2) irreps inside SU(3) is therefore given by inserting the lowest allowed value n * a , while n c remains 0. Since the number of columns in Young tableaus must not increase going from top to bottom, this value is determined by The lowest eigenvalue of effective SU(2) irreps inside SU (3) is thus given by (n * a ) and increase by p + 2 where p starts from 0 (or 1 for odd N s ) and increases by 2 for constant N s . Starting from SU(2), the quadratic Casimir operator only depends on p since the number of rows of the corresponding Young tableaus is 2. It can be written as which follows by setting S = p/2 into the usual eigenvalue expression C SU(2) 2 (S) = S(S + 1). For even N s the lowest eigenvalue is 0, while for odd N s it is 3/4. If p is increased, the eigenvalues of the SU(2) quadratic Casimir operator changes by Thus the slope is the same as for the SU(3) quadratic Casimir operator. So, every eigenvalue of the SU(2) quadratic Casimir operator also appears in SU (3), shifted by C SU(3) 2 (n * a ). E.g. for 21 sites, n * a = 11 and C SU(3) 2 (11) = 48, which perfectly agrees with Fig. 9. Mathematically, one may summarize the situation as follows. Let G ∈ {SU (2), SU (3)} be a Lie-Group and σ (G) denote the spectrum of the Quadratic Casimir operator of the group G, then: ∀s ∈ σ (SU (2)) ∃t ∈ σ (SU (3)) such that t − s = In Sec. IV, we argue that the six-fold degenerate ground state corresponds to an Abelian CSL. Here, we focus on the three variational states with π/3-flux and compare the numerical results from VMC with the predictions for the π/3-flux CSL.
We use the variational state overlap method [58][59][60] to confirm the topological properties of the anyonic excitations in the CSL phase. Similar VMC constructions were shown to produce topological fractional quantum Hall states in SU (2) systems [55,[61][62][63]. The information on the topological spins and mutual statistics of the anyonic quasi-particles can be determined by calculating the matrix elements of the generators of the modular transformations in the ground-state manifold on the torus [64]. Modular transformations are the mappings of the torus onto itself, mapping each vertex of the lattice to another one. All such transformations can be given as a composition of powers of two generators, S and T (Dehn twist), and they form the SL(2, Z) group. A general modular transformation can be represented by an integer valued 2×2 matrix with unit determinant, which describes how the vectors ω 1 and ω 2 defining the torus periodicity change under the transformation. The two generators are given by On the square lattice the S transformation is equivalent to a π/2 rotation around a site, but this is not true in general. In case of the triangular lattice S is not a symmetry transformation either, but T −1 · S is equivalent to a π/3 rotation, while T S gives a rotation by 2π/3. Calculating the matrix elements of S and T over the ground-state manifold in the basis of minimum entropy states (MES) [64] gives access to the topological spins and the exchange statistics of the anyonic quasi-particles where |Ξ a and |Ξ b are the MES [60]. The S and T transformations are not symmetries of the triangular lattice and they do not preserve the connectivity of the sites, therefore their matrix elements only give the desired topological quantities up to prefactors exponentially small in the system size [60]. In Eq. (D2), S ab gives the phase the ath quasi-particle acquires when going around the bth quasi-particle. In the MES basis the matrix of T is diagonal and the entries of θ a give the topological spin of each quasi-particle. One of the topological angles is always 1 corresponding to the trivial topological sector. The topological or chiral central charge c is the difference of the central charges for the left and right moving edge excitations [65]. However, the states obtained from the overlap calculations are not necessarily in the MES basis and diagonalizing the matrix of the T transformation does not immediately give the MES basis if several quasi-particles have the same topological spin. Therefore, some other criterion needs to be used to identify the MES basis. Before presenting the results from our VMC calculations, we determine the expectations for the π/3-flux CSL. On top of the trivial sector, the two other sectors correspond to two anyons with topological spins −2π/3, which we denote as a andā, emphasizing that they are the antiparticles of each other. The fusion rules are The Abelian nature is manifested in the fact that fusing two quasi-particles always results in only one quasi-particle. Based on the fusion rules and the topological spins, we can deduce the elements of the S-matrix based on the Verlinde formula [66] where N c ab are the coefficients in the fusion rules, d c are the quantum dimensions of each quasi-particle (1 for all quasiparticles in an Abelian theory), and D = d 2 c is the total quantum dimension. The topological spins in our case are θ 1 = 1 and θ a = θā = e −i2π/3 . The topological central charge c satisfies e ic2π/8 = α θ α d 2 α /D 2 [67,68], which is consistent with c = −2. The topological spin −2π/3 can be understood at the mean-field level. In order to exchange two fermions, they need to be moved around a triangular plaquette and thus the overall phase is a combination of the π phase of exchanging two fermions plus the extra π/3 phase picked up from the gauge flux. Also, for the π/3-flux case the fermions go around counter-clockwise in the bulk, thus the skipping modes on the edge travel clockwise, which agrees with c = −2. Based on this and Eq. (D4), the modular matri-ces T and S have the form and since it is an Abelian spin liquid all elements of S have the same magnitude.
Next, we present our VMC results for the matrix elements of the S and T modular transformations and show that they agree with the predictions for the Abelian CSL phase. The modular matrices are based on 10 10 and 10 9 independent measurements in 100 batches for the 6 × 6 and 9 × 9 clusters, respectively. The MES basis is identified with the case where, on top of the T -matrix being diagonal, all elements of the Smatrix have the same magnitude. In order to fix the overall phase of the minimally entangled states, one can use the fact that the row and column of the S matrix corresponding to the trivial sector should be 1/ √ 3. This yields  for the 6 × 6and 9 × 9-site system, respectively.
In the left panel of Fig. 20, we show how the real part of the prefactors Reµ (S) and Reµ (T ) scale with the system size.
Our results fulfill an L 2 scaling, and we find an estimate for the chiral central charge c ≈ −2.34, which is consistent with the expected c = −2. In the right panel of Fig. 20, we present the imaginary part of the prefactors Imµ (S) and Imµ (T ) corrected by the phase factors ±π/2 · L 2 , which does not change the extrapolated value. We note that the values of Imµ (S) and Imµ (T ) can be extracted only mod 2π. So the values for the different system sizes could be shifted by an integer multiple of 2π independently of each other. This can change the extrapolated value at L = 0 by an integer multiple of 2π/5, which would result in a shift of an integer multiple of 4.8 in the extrapolated chiral central charge c. So, this is clearly not the source of the discrepancy between the extrapolated and theoretical values.
We used the results from independent batches of calculations (compare Sec. III B) to estimate the statistical errors of the simulation. These are collected in Tab. I. However, the main source of error is systematic, and stems from the finite system sizes as well as from the fact that the projected states are not perfect topological CSL states. Carrying out calculations on larger systems could help to further verify our results, but based on the calculations on the 36-and 81-site clusters we find that for the 12 × 12 system the exponential prefactor would be O(10 −5 ). Since the error of Monte Carlo measurements scales as the inverse square root of the number of samples, it would require a minimum of 10 12 samples to see anything beyond the numerical error. This goes beyond the scope of our current VMC method.
The eigenvalues for the six chiral states on the 12-and 21site clusters are given in Tab. II and III, respectively. The results from the ±π/3 CSL variational states from VMC and from the low-lying energy eigenstates at K/J = 0.35 (α ≈ 0.11π) from ED show a perfect match.
The spectrum of the 12-site cluster is discussed in Sec. V B. On this particular system the ground-state manifold of six states is intertwined with another state. However, if the states are analyzed with respect to their symmetry values, the apparent six low-lying CSL states can be identified.
The symmetry properties of the chiral states for clusters larger than 21 sites were only studied by VMC. For systems where the vectors defining the torus lie in the r 1 and r 2 directions all six chiral states have a wave vector zero. The smallest example is the 36-site cluster for which we provide the symmetry properties in Tab. IV. For this specific finite size, only one pair of states is degenerate, the other states are not.
Appendix F: Effective model on 12-site cluster with PBCs On a cluster consisting of a finite number of sites the effective model is different from the model in the thermodynamic limit. In the case of PBCs this is due to the fact that a finite number of fermionic hoppings in one direction leads back to the starting site. For the 12-site torus this becomes relevant in order 4 in t/U , where the four-site plaquette can be embedded surrounding the cluster via the PBCs. This leads to an additional effective interaction around the 12-site torus, but also to modifications of other coupling constants compared to the infinite lattice. The easiest approach to derive the effective Hamiltonian for the most interesting case Φ = π is to consider the model without phases and then perform the transformation t → −t, which is identical to Φ = 0 → Φ = π. For general fluxes Φ the calculation is done by choosing the gauge for the 12-site cluster depicted in Fig. 21 where the depicted graphs under every sum indicate in which orientations the exchanges lie on the cluster. For the three site ring exchange on a triangle K all triangles have to be taken. The four-site ring exchange around the torus is only stated. The explicit shapes differ and lead to either a phase factor or not; In total every site is part of 18 four-site ring exchanges around the PBCs, wherein 10 contribute without a phase factor and 8 with a phase factor. In total these exchanges double due to the hermitian conjugated exchanges. For the real exchange constants (J, L 2sp s , L 2sp d , L 3sp d, vert , L 4sp r, pbc,Φ=0 ) the turning direction is arbitrary, since the Cosine is symmetric and the coupling constants are identical under the transformation Φ → −Φ. For the exchanges contributing with a non-trivial phase factor the turning direction has to be chosen accordingly. For 2-site interactions such a turning direction cannot be defined, which is why such couplings contribute with purely symmetric phase factors.
The effective nearest-neighbor exchange constants are given in Eq. (F1). All other couplings and the constant contributions are  (F3) Note that one could also derive an effective model for the 12-site cluster in which neither phases around the torus nor distinct exchange constants for topologically equivalent interactions (like for J hor and J dia ) occur. This can be achieved by fulfilling the Φ-flux condition on every triangle without assigning specific phases to specific bonds. However, this does not matter since the eigenenergies do not depend on the gauge. For the J-K model at Φ = π/2, hence for purely imaginary ring exchange, a π/3-flux CSL was discovered for SU(3)symmetric spins in Ref. 37. As for the spontaneous timereversal symmetry breaking CSL found for the same model but at Φ = π, discussed in the main part of the manuscript, the relevant question is whether this CSL is a feature of the Hubbard model, and therefore reachable in experiments with artificial gauge fields.
We first study the fifth-order effective model in Eq. (14) for Φ = π/2 and find that the same set of Padé extrapolations as described for Φ = π in Sec. V A works best. The convergence behavior of all coupling constants is shown in Fig. 22. A number of interesting and subtle features of the effective model become clear. The three-site ring exchange K is purely imaginary only in order 3. The fourth-order term partly arises from fluctuations around two triangles leading to a flux of 2Φ, therefore a real part is present in higher orders. Similarly, the imaginary part of the fourth-order contribution to the four-site ring exchange vanishes and it effectively becomes an order 5 term. As a consequence the model is dominated by a real nearest-neigbor and an imaginary three-site ring exchange. For exactly this subset of interactions the CSL phase occurs for Im(K/J) 0.3 [37]. In Fig. 23 the ratios of Padé extrapolants for the imaginary part of −K and the real amplitude J (green larger diamonds) and the direct Padé extrapolation of the ratio −ImK/J (black smaller diamonds) are shown. We see that in the regime of the CSL not only the extrapolations, but also the bare fifth-order (solid line) series are well converged.
Second, the numerical study of the full fifth-order effective model yields the signature of the CSL phase by three low-lying chiral states. Using ED on the 21-site cluster we find (t/U ) ≈ 0.04. This is plausible since the VMC captures CSL phases more naturally than long-range ordered phases. The apparent areas in Fig. 22 and Fig. 23 are shaded in yellow.
This CSL is stable in the J-K model for fluxes π/2 ≤ Φ ≤ π. The fifth-order effective model in this range partly suffers from spurious poles in the Padé extrapolants of several coupling constants. More precisely, for L 2sp d around Φ ≈ 0.6π and for L 3sp s around Φ ≈ 0.55π and Φ ≈ 0.75π divergences occur, depending on t/U . Here one has to use the bare series, which generally decreases the quality of convergence.
Furthermore, similar to Φ = π/2, cancellation effects of different orders at specific values of the flux Φ lead to additional subtleties in the convergence behavior. For most parts of the full Φ phase diagram we find a good convergence by comparing the ED and VMC results in different orders. It is only in the area around Φ ≈ 3π/4 that the results between fourth-and fifth-order change qualitatively. A possible explanation is the strong suppression of a lower-order term compared to a higher-order one for certain values of Φ and t/U , which can lead to a decreased quality of convergence. Here the real part of the leading fourth-order ring exchange vanishes in order 5 due to an alternating behavior by contrast to a large monotonic five-site ring-exchange. The interplay of these two effects arises only in order 5.
In summary, the π/3-flux CSL phase at Φ = π/2 is extended in the J-K model for larger fluxes Φ up to the spontaneous time-reversal symmetry broken case at Φ = π. At Φ = π/2 the CSL is a robust feature within the effective fifth-order model describing the SU(3) Hubbard model in the strong-coupling regime. If the metal-insulator transition occurs at values similar to the ones at Φ = π, the π/3flux CSL is also a feature of the SU(3) Hubbard model at Φ = π/2 on the triangular lattice. Whether this phase is directly connected to the spontaneous time-reversal symmetry broken CSL present at Φ = π remains an open question.