Topological chiral spin liquids and competing states in triangular lattice SU($N$) Mott insulators

SU($N$) Mott insulators have been proposed and/or realized in solid-state materials and with ultracold atoms on optical lattices. We study the two-dimensional SU($N$) antiferromagnets on the triangular lattice. Starting from an SU($N$) Heisenberg model with the fundamental representation on each site in the large-$N$ limit, we perform a self-consistent calculation and find a variety of ground states including the valence cluster states, stripe ordered states with a doubled unit-cell and topological chiral spin liquids. The system favors a cluster or ordered ground state when the number of flavors $N$ is less than 6. It is shown that, increasing the number of flavors enhances quantum fluctuations and eventually transfer the clusterized ground states into a topological chiral spin liquids. This chiral spin liquid ground state has an equivalent for the square lattice SU($N$) magnets. We further identify the corresponding lowest competing states that represent another distinct type of chiral spin liquid states. We conclude with a discussion about the relevant systems and the experimental probes.


I. INTRODUCTION
SU(N ) Mott insulators are representative examples of quantum systems with a large local Hilbert space where quantum fluctuations can be strongly enhanced and exotic quantum phases could be stabilized. All the SU(N ) spin operators are present in the effective model for the Mott insulators and would be able to shuffle all the spin states rather actively in the local Hilbert space. The system can be quite delocalized within the local Hilbert space, that is to enhance quantum fluctuations and induce exotic quantum phases. This aspect is fundamentally different from the SU(2) Mott insulators with large-S local moments that also has a large local Hilbert space. For the large-S SU(2) Mott insulators, the pairwise Heisenberg interaction is quite ineffective to delocalize the spin states in the large-S Hilbert space, and thus quantum fluctuations are strongly suppressed. Thus, it is the conventional wisdom not to search for exotic quantum phases among the large-S SU(2) Mott insulators, but among the spin-1/2 quantum magnets with a strong frustration. In contrast, the emergence of the SU(N ) Mott insulators brings a new searching direction for exotic quantum phases. SU(N ) Mott insulators are not a theoretical fantasy, but exist in nature. It has been shown that the ultracold alkaline-earth atoms (AEA) on optical lattices can simulate quantum many-body physics with an SU(N ) symmetry without any fine-tuning [1]. The nuclear spin of fermionic AEA can be as large as I = 9/2 for 87 Sr, while the outer shell electrons give a total spin S = 0 and makes the hyperfine coupling inactive. This observation effectively extends the realization of SU(N ) magnets up to * gangchen@hku.hk N = 2I + 1. Early efforts by Congjun Wu explored total spin-3/2 alkaline fermions on optical lattices where the SU(4) symmetry can be achieved with fine tuning [2][3][4][5]. Quantum Monte Carlo simulations were introduced later to study magnetic properties of the SU(2N ) Hubbard model [6,7]. There is also some effort in searching for an emergent SU(N ) symmetry in real materials particularly for the SU(4) case. The two-orbital Kugel-Khomskii model can become SU(4) symmetric after some fine tuning [8,9]. Experimental and numerical evidence also suggests that Ba 3 CuSb 2 O 9 could be a prominent candidate on a decorated honeycomb lattice [10], though the Cu-Sb dumbbell is quenched rather than an active degree of freedom. The SU(4) Heisenberg model has further been proposed for the spin-orbit-coupled Mott insulator α-ZrCl 3 where the degree of freedom is the spinorbit-entangled J = 3/2 local moment [11] on a honeycomb lattice [12]. More recently, it has been proposed that the Mott insulating and superconducting behaviors in twisted bilayer graphenes can be captured by a twoorbital Hubbard model with an emergent SU(4) symmetry on a moiré triangular lattice [13]. Other work even suggested that double moiré layers, built from transition metal dichalcogenides or graphenes, separated from one another by a thin insulating layer would be a natural platform to realize Hubbard models on the triangular lattice with SU(4) or SU(8) symmetries [14,15]. Owing to enhanced quantum fluctuations for SU(N ) Mott insulators, the pioneering theoretical works [16,17] by Hermele et al, have obtained the topological chiral spin liquid (CSL) ground states with intrinsic topological orders even for an unfrustrated square lattice when N ≥ 5 [16][17][18]. Depending on the atom occupation numbers, the system can support both Abelian and non-Abelian statistics for the anyonic excitations. Since unfrustrated lattices such as square and honeycomb lattices already bring exotic and interesting physics [19][20][21], the frustrated lattices could further harbor nontrivial quantum phases, for example the triangular lattice with SU(3) and SU(4) spins [22][23][24]. In this work, we focus on the SU(N ) Heisenberg model on a triangular lattice where each lattice site comprises the fundamental representation of the SU(N ) group. For AEA, this corresponds to one atom per lattice site and is known to be most stable against the three-body loss. Because the 1/N filling is kept throughout, the large-N limit differs fundamentally from the large-N extension of the spin-1/2 SU(2) models where the 1/2 is kept and two sites can form a SU(2) spin singlet [25]. Here, as N sites are needed to form a singlet, the valence cluster solid (VCS) state is generically disfavored in the large-N limit. Instead, by our large-N calculation, two types of CSL states with background U(1) gauge flux φ = π − π/N and π − π/(2N ) (see Fig. 1(a)) are identified as the ground and lowest competing states for 6 ≤ N ≤ 9, respectively. For smaller N 's, various symmery-broken cluster/stripe states are obtained. We expect the large-N results are more reliable when N is large.
The rest of this paper is organized as follows. The general Heisenberg model of SU(N ) spins is introduced and simplified at the large-N saddle point in Sec. II. The self-consistent minimization algorithm is implemented to solve the reduced mean-field spinon Hamiltonian. The technique details of this algorithm are described in Sec. III. It is emphasized that the optimized solutions strictly satisfy the local constraints. In Sec. IV, both the ground states and lowest competing states for 2 ≤ N ≤ 9 are reported and analysed, especially for two types of CSL states. Finally, in Sec. V, this work conclude with a discussion about relevant systems and the experimental probes.

II. THE SU(N ) HEISENBERG MODEL IN THE LARGE-N APPROXIMATION
We begin with an SU(N ) Heisenberg model on the triangular lattice where each site comprises the fundamental representation of the SU(N ) group. This model can be obtained from the strong coupling limit of an SU(N ) Hubbard model with 1/N filling or one particle per site. The SU(N ) Heisenberg model is given as where J is the antiferromagnetic exchange interaction and the sum is taken over the nearest neighbor bonds and spin flavors. The SU(N ) spin operators can be simply expressed with the Abrikosov fermion representation S αβ (r) = f † rα f rβ , and α, β = 1, . . . , N . Hereafther, a summation over repeated indices in the form of Greek letters is supposed unless otherwise specified. A local constraint on the fermion occupation f † rα f rα = 1 is imposed to reduce the enlarged Hilbert space. In principle, an SU(N ) singlet could be formed by N sites. But the SU(N ) exchange rapidly transforms one N -site singlet to other sets of N -site singlets, resulting in a failure of the conventional understanding. Instead, the spin Hamiltonian Eq. (1) is solvable in the limit N → ∞ via a large-N saddle point approximation in the imaginary-time path integral formulation. The corresponding partition function can be expressed as The action S is given as where µ r is the Lagrangian multiplier to enforce the Hilbert space constraint, χ rr is the auxiliary field to decouple the fermion operators, and J ≡ N J. As the action S scales linearly with N , the large N limit leads to a saddle point approximation that results in the saddle point equations, χ rr = −J f † r α f rα /N , f † rα f rα = 1, and the saddle point or mean-field Hamiltonian for the fermionic spinons is In the following, we will search for the saddle point with the lowest mean-field energy E MF numerically and discuss the properties of ground states. Before that, we first discuss the lower bound of E MF and the bound saturation conditions for reference. An exact lower bound on E MF for generic lattices [26] was first obtained by Rokhsar for the half-filling, and was shown to be saturated by valence bond states with various spin singlet coverings. These states break the lattice translation but preserve the spin-rotation symmetry, and fluctuations beyond the mean field can break the high degeneracy among them [25]. The lower bound was generalized to the 1/N filling with [17] where N s is the number of lattice sites.
We set J max = J for each bond. The saturation for Eq. (5) is reached by a N -simplex VCS state composed by N -site simplices. That is, every site on the lattice is directly connected to the other N − 1 sites within the same cluster by a single bond with an exchange coupling J max . On a d-dimensional lattice, the N -simplex VCS with N > d + 1 is prohibited without fine-tuning of the exchange. Thus, there are only two-simplex and threesimplex VCS's on the triangular lattice. For N ≥ 4, possible cluster states are general N -site ones satisfying a stricter energy bound [17] is the total number of isolated bonds in the lattice.

III. THE SCM ALGORITHM
To determine the saddle-point solutions, we closely follow the numerical self-consistent minimization (SCM) algorithm developed in Refs. 16 and 17. The technical details of the algorithm are described in the following. During one energy optimization process for a given geometry and periodic boundary conditions, the algorithm starts from initializing the fields χ rr at each bond randomly with χ rr = |χ rr |e iφ rr , with a uniform distribution of amplitudes |χ rr | ∈ [0.02, 0.20] and phases φ rr ∈ [0, 2π]. The chemical potentials are set to be the default value µ r = 0 in the beginning. Obviously, the local constraints are violated here in general and the deviation of the local fermion density can be denoted as δn r = 1 − f † rα f rα . The expectation value is taken using the ground state of H MF with the unchanged µ r at this stage. To obtain the correct density, one must adjust the chemical potential µ r by δµ r at each site. It is found that to the lowest order, the desired adjustment of chemical potentials δµ r can be expressed as where G −1 (r − r , 0) is nothing but the inverse of the density-density correlation at zero frequency. The elements of the correlation G(r − r , 0) can be calculated at the single particle level from the mean-field Hamiltonian H MF , with the help of the standard linear response theory in principle. However, it follows that the inversion of G(r − r ) naively diverges. Physically, this is because a uniform adjustment of the chemical potential at every site is trivial and the fermion density will keep intact. This obstacle is overcome by following the treatment suggested in Ref. 17; that is diagonalizing G(r − r , 0) and only focusing on its non-zero eigenvalues g i . Note that the index i refers to the basis where G(r − r , 0) is diagonalized. In such a basis, the relationship Eq. (6) becomes Here, the index i should not be summed. For all indices with vanishing eigenvalues g i , µ i is taken to be zero concurrently. Then, the adjustment of the chemical potential δµ r is well-defined. A simple replacement of µ r → µ r + δµ r gives a new mean-field Hamiltonian H MF and related ground state. In consequence, it affects the local fermion density conversely, resulting in a new deviation δn r . These processes construct a self-consistent procedure and the problem of searching for an appropriate set of chemical potential deviation δµ r can be solved by iterating the procedure until the density is uniform. This is the core of the algorithm to preserve the local single-occupation constraints n r = f † rα f rα = 1. With the modified chemical potentials calculated in the previous step, an update set of auxiliary fields χ rr can be determined via χ rr = −J f † r α f rα /N , but the local constraints may be violated once again. The amended auxiliary fields and chemical potentials can be treated as a new and better starting point. The procedure described in the previous step is thus implemented iteratively until reaching a converge of the ground state energy within a given numerical error. It has been shown in Ref. 17 that the energy of the final state must be less or equal to that of the initial state after the optimization process. Therefore, we eventually obtain a local energy minimum. In order to reach the global minimum as much as possible, we start from at least 50 randomly initialized fields, and reap a collection of local minima satisfying the singleoccupation constraints. The lowest two are accepted as the best results of the ground state and lowest competing state energies.

IV. THE MEAN-FIELD RESULTS
We describe the ground states and lowest competing states of the mean-field Hamiltonian Eq. (4) from the SCM algorithm. Because different geometries can accommodate different candidate ground states especially cluster states with unknown dimensions, in the calculation, we consider all unit cells of a parallelogram geometry 1 × 2 with 1,2 ≤ 2N for 2 ≤ N ≤ 5 and 1,2 ≤ N for N ≥ 6, respectively. In some ordered cases for N = 4, 5, larger unit-cell sizes are also considered for confirmation. Each unit cell is repeated by L 1,2 (≥ 20) times along the triangular Bravais lattice vectors to form a superlattice. The superlattice itself has periodic boundary conditions. In practice, for each case, we ran the SCM procedure at least 50 times with different random seeds to avoid any missing of ground states caused by numerical problem. The results are listed in Table I and Table II. For N = 2 and 3, the lowest energies we found saturate the bound Eq. (5), meaning that the ground states are dimer and three-site simplex VCS, respectively. Actually, the obtained ground states are highly degenerate because any state with each site covered by one dimer/cluster unit has the same energy in the large-N limit. An ordered dimer or three-site simplex VCS state, as illustrated in Figs. 2(a-b), is expected to be selected beyond the mean field. Note that, the non-zero expectation values |χ rr | in Fig. 1(a) and (b) can only take 1/2 and 1/3, respectively. The true ground state for N = 2 is the wellknown 120 • Néel state and differs from the mean-field here. The N = 3 case was also shown to have a threesublattice magnetic order in previous numerics [27,28]. This reflects the deficiency of the large-N approximation for small N cases. In Table I, we further find that the lowest competing states are CSL states with φ = π−π/N . The N = 4 case is more compelling due to the rapid growth of experimental proposals [14,15,29,30] and large-N approximation could provide useful insight here. We find the ground state energy do not saturate the bound Eq. (5) as excepted, but saturate the stricter bound discussed previously. The ground states are foursite VCS depicted in Fig. 2(c), accompanied by a large degeneracy for the same reason as N = 2 and 3. The non-zero expectation values are |χ rr | = 1/4. A similar plaquette order was also reported in a recent work [15]. In Fig. 3, we depict the lowest competing state for N = 4. One can see that the lattice is covered by stripes with three different bond expectation values |χ rr | ≈ 0.030, 0.158, and 0.224. Along one of lattice vectors, there is a unit-cell doubling that breaks the lattice translation. The background fluxes through each plaquette are inhomogeneous as well and manifest the same unit-cell doubling pattern. Specifically, the average flux is a constant value φ avg = π − π/N where N = 4. The same ordered pattern has also been obtained through the DMRG study from the Kugel-Khomskii model [31]. They attribute the symmetry breaking to symmetry allowed Umklapp inter-actions in certain finite geometries. In our results, however, such a stripy state is not the ground state even in the two-dimensional limit.
From N = 5, neither the bound Eq. (5) nor the stricter one can be saturated. Thus, the ground states are no longer VCS. Our numerical calculation finds that the ground state for N = 5 is very similar to the lowest competing state for N = 4 (see Fig. 2(d)), except that the bond expectation values can only take |χ rr | ≈ 0.072, 0.145, and 0.179, and the average background flux φ avg = π − π/N shifts accordingly. It also breaks the lattice translation symmetry along one of lattice vectors and manifests itself as a stripe pattern with a doubled unitcell. The lowest competing state for N = 5 is a CSL with φ = 4π/5. With further increasing N , the frustration is enhanced. Eventually, the ground states become CSL states with φ = π − π/N when 6 ≤ N ≤ 9. Correspondingly, the lowest competing states we found share the identical form. They are CSL states as well except that the background magnetic fluxes shift to φ = π − π/(2N ).
With the two types of CSL states identified, we now discuss the properties of these topological liquid states. The CSL is characterized by the mean-field saddle point All bonds on the lattice have a uniform expectation value for |χ| but modulated by a U(1) gauge field a rr so that the flux φ on each plaquette is a constant. The CSL breaks both parity and time-reversal. The bond phase field a rr is treated as a fluctuating U(1) gauge field coupled to the fractionalized spinons. By checking the behavior of the mean-field Hamiltonian Eq. (4) at the CSL saddle points, we find that both types of CSL states have a fermion band structure with N bands where only the lowest is filled (see Fig. 1(b)). The Fermi level lies in the gap between the lowest two bands, thus all discussion can be applied to both CSL states. Furthermore, the first type of CSL with φ u,d = π − π/N on the triangular lattice can be mapped to the counterpart on a square lattice up to a time-reversal. If we regard the adjacent up and down triangles as a unit shown in Fig. 4, the phase of the shared bond has no contribution to the total flux, and we have the relationship where φ s is the background U(1) flux through each square plaquette. In Ref. [16], Hermele et al, found that, the CSL states are ground states for 5 ≤ N ≤ 10 on the square lattice in the mean-field calculation. As the spinon is gapped out by the emergent U(1) gauge flux pattern, the Chern-Simons term enters the theory for U(1) gauge fluctuations. After integrating the gapped spinons, one obtains a topological quantum field theory with a Chern-Simons term corresponding to the chiral Abelian topological order and anyonic statistics. The spinon is converted in anyons with a statistical angle π ± π/N . Gapless chiral states carrying spin degree of freedom are also supported by the CSL as edge modes and the low-energy theory is described by the SU(N) 1 WZW model.

V. DISCUSSION
To summarize, we study the Heisenberg antiferromagnets with SU(N ) symmetry on the triangular lattice. In the large-N approximation, a variety of ground states and lowest competing states are identified for different parameter N . At the mean-field level, the ground state for 2 ≤ N ≤ 4 is a N -site VCS state with a large degeneracy. Specifically, ordering patterns with doubled unit cell and average background flux φ avg = π − π/N are found in the N = 4 and 5 cases. These ordered states break the lattice translation symmetry along one of lattice vectors and become the lowest competing state and ground state for N = 4 and 5, respectively. The frustration from SU(N ) exchange interaction is enhanced when N > 5, resulting two types of CSL states as the lowest two states for 6 ≤ N ≤ 9. Among them, the CSL states with φ = π − π/N have a lower energy, and have its counterpart on the square lattice. Although the true ground states are not what we found for N < 4, our calculation can provide useful insight for cases N ≥ 4 where large-N approximation becomes more reliable. In a very recent DMRG study of an SU(4) spin model on the triangular lattice, phase diagrams for integer fillings were obtained and compared with conventional MFT ones [15]. The quite good agreement of the phase boundaries determeted by two methods suggests that N = 4 is perhaps large enough for the mean-field analysis we performed in this work.
Thanks to the development of ultracold experimental techniques, the SU(N ) Mott insulators have been realized with AEA on various optical lattices using the Pomeranchuk cooling [32], even the Mott crossover and SU(N ) antiferromagnetic spin correlations were recently observed with 173 Yb atoms [33,34]. Nontrivial physics of multicomponent fermions with broken SU(N ) symmetry were also proposed on this platform [35]. The emergency of the widely concerned SU(4) or even SU (8) symmetric interaction has been proposed in the twisted bilayer graphene and double moiré layers systems re-cently [12,14]. This moiré lattice system could be an novel platform to detect possible phases in this work. The ultracold atom systems may have many restrictions in the detection of anyonic spinon excitations and edge states. Nevertheless, spin-dependent Bragg spectroscopy may be used to detect the spinon continuum [16][17][18], singlet-triplet oscillation technique can detect the nearest-neighbor spin correlation [34] for CSLs, and the lattice potential could be adjusted to localize and manipulate the anyonic quasiparticles [16][17][18]. For solid-state systems, more techniques are available. These include but are not restricted to quantized thermal Hall transport [36] for the edge modes, scanning tunneling microscope of anyons at defects [37], or even angleresolved photon emission measurement for the spinon signatures [38].