SU$(3)_1$ Chiral Spin Liquid on the Square Lattice: a View from Symmetric PEPS

Quantum spin liquids can be faithfully represented and efficiently characterized within the framework of Projected Entangled Pair States (PEPS). Guided by extensive exact diagonalization and density matrix renormalization group calculations, we construct an optimized symmetric PEPS for a SU$(3)_1$ chiral spin liquid on the square lattice. Characteristic features are revealed by the entanglement spectrum (ES) on an infinitely long cylinder. In all three $\mathbb{Z}_3$ sectors, the level counting of the linear dispersing modes is in full agreement with SU$(3)_1$ Wess-Zumino-Witten conformal field theory prediction. Special features in the ES are shown to be in correspondence with bulk anyonic correlations, indicating a fine structure in the holographic bulk-edge correspondence. Possible universal properties of topological SU$(N)_k$ chiral PEPS are discussed.

Introduction -Quantum spin liquids are long-range entangled states of matter of two dimensional electronic spin systems [1][2][3]. Among the various classes [4], spin liquids with broken time-reversal symmetry, i.e., chiral spin liquids (CSL) [5,6], exhibit chiral topological order [7]. Intimately related to Fractional Quantum Hall (FQH) states [8], CSL host both anyonic quasi-particles in the bulk [9] and chiral gapless modes on the edge [10]. It was early suggested that, in systems with enhanced SU(N ) symmetry, realizable with alkaline earth atoms loaded in optical lattices [11], CSL can naturally appear [12]. Later on, many SU(N ) 1 CSL with different N were identified on the triangular lattice [13], while the original proposal on the square lattice [12] remains controversial.
In recent years, Projected Entangled Pair States (PEPS) [14] have progressively emerged as a powerful tool to study quantum spin liquids. As an ansatz, PEPS provide variational ground states competitive with other methods [15,16], and equally importantly, offer a powerful framework to encode topological order [17] and construct non-chiral [18,19] and chiral -both Abelian [20] and non-Abelian [21] -SU(2) spin liquids. Generically, SU(2) CSL described by PEPS exhibit linearly dispersing chiral branches in the entanglement spectrum (ES) well described by Wess-Zumino-Witten (WZW) SU(2) k (with k = 1 for Abelian CSL) conformal field theory (CFT) for one-dimensional edges [22]. However, to our knowledge, there is no known example of more general SU(N ) PEPS with unambiguous chiral edge modes. Thus it remains unclear whether symmetric PEPS can describe higher SU(N ) CSL faithfully. In order to address these issues, we propose and investigate a frustrated SU(3) symmetric spin model on the square lattice with a symmetric PEPS ansatz, thereby taking the first step towards describing general SU(N ) k CSL with PEPS.
Model and exact diagonalization -On every site, we place a three-dimensional spin degree of freedom, which transforms as the fundamental representation of SU (3). The Hamiltonian, defined on a square lattice, includes the most general SU(3)-symmetric short-range three-site interaction: where the first (second) term corresponds to two-site permutations over all (next-)nearest-neighbor bonds, and the third and fourth terms are three-site (clockwise) permutations on all triangles of every plaquette. We have chosen J 2 = J 1 /2 so that the two-body part (J 1 and J 2 ) on the interacting triangular units becomes S 3 symmetric, hence mimicking the corresponding Hamiltonian on the triangular lattice [23] and further parameterized the amplitude of each term as: J 1 = 2J 2 = 4 3 cos θ sin φ, J R = cos θ cos φ, J I = sin θ. We have performed extensive exact diagonalization (ED) calculations [24,25] on various periodic N s -site clusters to locate the CSL phase in parameter space. We expect (1) to host a SU(3) 1 CSL equivalent to the 221 Halperin FQH state [26], whose spectral signatures on small tori can be established precisely [27,28]. A careful scan in θ and φ reveals that there is a small region where, for N s = 3p (p ∈ N + ), there are three low-lying singlets below the octet gap, reflecting the expected topological degeneracy of the CSL. For N s = 3p, the lowenergy quasi-degenerate manifold reflects perfectly the anyon content of the CSL. In both cases, the momenta of the lowenergy states match the heuristic counting rules of the 221 Halperin FQH state with 0, 1 or 2 quasi-holes [28-31] (see Fig. 1 and supplemental material (SM) [32]). In the following, we shall focus on angles θ = φ = π 4 , for which clear evidence of a gapped CSL is found. [33] Symmetric PEPS ansatz -Representing the CSL in terms of a symmetric PEPS allows us not only to obtain short-range properties such as energy density efficiently, but also to reveal its topological properties by examining the entanglement structure [34,35]. This can be accomplished by using SU(3)symmetric tensors, analogously to the SU(2) case [36]. The simplest virtual space available here is V = 3 ⊕ 3 ⊕ 1 such that (i) a symmetric maximally entangled SU(3) singlet |Ω can be realized on every bond by pairing two neighboring virtual particles and (ii) four virtual particles around each site can be fused into the 3-dimensional physical spin with an on-site projector P [37, 38], see Fig. 2(a). The so-called bond dimension is thus D = 7. In addition to the continuous rotation symmetry, full account of the discrete C 4v point group symmetry (shown as purple arrows in Fig. 2(a)) can be taken [36], and tensors can be classified according to the corresponding irreps. By linearly combining on-site projectors of two different irreps with opposite ±1 characters w.r.t. axis reflections, one can construct a complex PEPS ansatz breaking both parity (P) and time-reversal (T) symmetries while preserving PT, as required for a CSL ground state. Details about this construction, used for SU (2) CSL [20,21], are provided in the SM. For the ease of the following discussion, it is convenient to define the tensor A by absorbing the adjacent singlets on, e.g., the right and down bonds around each site into the on-site projector, forming an equivalent way to express the wavefunction.
The fact that center of the SU(N ) group is isomorphic to the Z N group allows one to associate a Z 3 charge Q = +1 to the physical space of the tensor A, while the virtual space carries Z 3 charges Q = {+1, −1, 0}, i.e. it contains a regular representation of Z 3 [39]. Hence, the tensor A bears an important Z 3 gauge symmetry associated to local charge conservation: where the action on virtual indices reads as left, right, up and down, and ω = e i2π/3 , Z = diag(ω, ω, ω, ω 2 , ω 2 , ω 2 , 1) is the representation of the Z 3 generator in V. This built-in gauge symmetry is central to topological properties, such as topological degeneracy on the torus and anyonic excitations [17]. Let us emphasize that the Z 3 gauge symmetry naturally appears from the physical SU(3) and point group symmetry, and is not a symmetry we imposed ad hoc.
Variational optimization -The best variational ground Ns, respectively. The PEPS energy is optimized at χ = D 2 and extrapolated to χ → ∞ (red circles). Blue squares stand for DMRG data on several finite-width cylinders (Nv = 3,4,5,6), and the ED results on tori with Ns = 12, 15, 18, 21, 24 sites and different geometries are indicated by stars. The dotted (dash-dotted) line is an exponential fit of the DMRG (ED) data to the thermodynamic limit. state is obtained by taking the ansatz P = N1 a=1 λ a is the number of linearly independent projectors in the B 1 (B 2 ) class, and optimizing the (few) variational parameters {λ a 1 , λ b 2 } ∈ R with a conjugate-gradient method [40]. For a given tensor the energy is obtained via the corner transfer matrix renormalization group (CTMRG) method, computing an effective environment of bond dimension χ surrounding an active 2 × 2 region embedded in the infinite plane (so-called iPEPS) [41][42][43][44]. The gradient is then simply obtained by a finite difference approach [45]. U(1) quantum number is also used occasionally to speed up the computation [46][47][48]. The exact contraction scheme corresponds to the limit χ → ∞.
To establish the relevance of our symmetric PEPS ansatz for the model (1), we compare the PEPS energy density with that obtained by ED on several different tori up to size N s = 24, and by the density matrix renormalization group (DMRG) method [49, 50] on various finite cylinders. In DMRG, for each cylinder width N v , we compute the ground state energies at two different cylinder lengths to subtract the contribution from the edges [51]. A detailed description of the DMRG method and additional data can be found in the SM. As shown in Fig. 2(b), the PEPS energy density obtained on the infinite plane turns out to lie close to the energy density in the thermodynamic limit, estimated from a finite-size scaling of the ED and DMRG data.
Entanglement spectrum -To get further insight into the nature of the CSL phase, we now explore the properties of our symmetric PEPS, where the Z 3 gauge symmetry implies topological degeneracy on closed manifolds. On finite-width cylinders, quasi-degenerate ground states can be constructed by restricting the virtual boundary of PEPS to fixed Z 3 charges Q = 0, ±1, with or without inserting Z 3 flux line through the cylinder. Here we shall focus on states without Z 3 flux line, and briefly discuss those with flux line in the SM. The topological properties can be most easily obtained through a study of the entanglement spectrum, which is defined to be minus log of the spectrum of reduced density matrix (RDM) of subsystem, say the left half of a cylinder [34]. For a PEPS on an infinite long cylinder, the RDM can be constructed from the leading eigenvector of the transfer operator [35]. Since the onsite tensor carries charge +1, the cylinder width N v must be multiple of three. In our current setting with bond dimension D = 7, exactly contracting the transfer operator is not feasible for large enough N v . Instead we use the iPEPS environment tensors computed with CTMRG to construct the approximate leading eigenvector [52], where large enough environment dimension χ is needed to get converged results [21]. The constructed RDM is fully invariant under translation and SU(3) rotations, which allows to block diagonalize it, introducing appropriate quantum numbers. In practice, we use the Z 3 charges (associated to the Z 3 gauge symmetry), two U(1) quantum numbers, and the momentum quantum number to do ED. The results with N v = 6, χ = 343 are shown in Fig. 3 for the three different charge sectors, i.e., Q = 0, ±1.
Linearly dispersing chiral modes well separated from the high energy continuum are seen with the same velocity, one mode in the Q = 0 sector and three modes in the Q = +1 sector. The Q = ±1 sectors have identical spectra: As both the bare tensor and the bond |Ω are PT-symmetric, so is the wavefunction, but after the reflection, the bonds are at the other side of the entanglement cut, and since the bonds exchange 3 ↔ 3, this maps between the Q = ±1 spectra. Interestingly, for all different χ we have considered, the lowest level in the  sectors. The spectrum in Q = −1 sector (not shown here) is found to be identical to that in Q = +1 sector but with conjugated SU(3) irreps, and is shown in the SM. For convenience, the lowest eigenvalue is subtracted in each plot. One chiral branch is seen in (a) starting at momentum K0 = −π/3 and three branches are seen in (b) starting at momenta K±1 = −π/3, π/3 and π. In each sector, the irreps encircled by the red boxes (or the blue boxes and the arrows) agree with the level counting of the SU(3)1 WZW CFT (shown on the plot vertically). Q = 0 sector appears at finite momentum K 0 = −π/3, while the three branches in the Q = ±1 sectors start at momenta K ±1 = −π/3, π/3, and π. We believe the momentum shift is due to a quantum of magnetic flux trapped in the cylinder, and is an intrinsic property of the optimized PEPS, that constrains us to choose N v = 6p, p integer (For N v = 3(2p + 1), where K 0 and K ±1 do not belong to the reciprocal space, see SM).
Reconstructing the SU(3) irreps from the two U(1) quantum numbers (the Young tableaux for the relevant SU(3) irreps are provided in the SM), we found that the level contents follow the prediction of the Virasoro levels of the SU(3) 1 WZW CFT [22,53]. However, we observe a tripling of the branches in the Q = ±1 sectors that we shall discuss later.
Bulk correlations -The above entanglement spectrum provides strong evidence of SU(3) 1 chiral topological order. However, it has been shown that in PEPS describing chiral phases, certain bulk correlation lengths computed from the transfer matrix spectrum diverge [20,21,52,[54][55][56][57][58][59]. Nevertheless, a priori it is not known which type of correlation is quasi long-ranged, and how critical bulk correlations are related to the observed chiral edge modes. Here we address this question with our symmetric PEPS ansatz, where both the SU(3) symmetry and the associated Z 3 gauge symmetry will provide valuable insights.
Within the PEPS methodology, correlation lengths of different types of operators, including the anyonic type, can be obtained from two complementary methods. On one hand, correlation functions of usual local operators, e.g., spin-spin can be obtained directly by applying the local operator on the physical indices.
Here the spin operators are the eight generators of su(3) algebra in fundamental representation (see SM for explicit expression), and the dimer operator is defined as D x i = S i · S i+dex . The Z 3 gauge symmetry enables to define topologically nontrivial local excitations like spinon, vison and their bound state [17,60,61] . A local excitation in the spinon sector can be created by applying an operator X satisfying XZ = ωZX on the virtual indices of the local tensor such that it carries zero Z 3 charge instead of the original charge 1. Similarly, X 2 can create a charge −1 spinon, since X 2 Z = ω 2 ZX 2 . A pair of vison excitations can be created by putting a string of Z (or Z 2 ) operators on the virtual level, whose end points correspond to the vison excitations. Parafermions, bound states of a spinon and a vison, can be created by putting spinons at the end points of the Z string. All these real space correlations can be obtained using the CTMRG environment tensors, see SM for further details. The correlations of the different types are shown in Fig. 4(a), computed with χ = 392.
On the other hand, correlation lengths can also be extracted from the spectrum of the transfer matrix of the environment tensors provided by CTMRG (see SM), also termed channel operator, whose eigenvalue degeneracies carry information about the types of correlation. Correlation lengths along horizontal and vertical directions are found to be the same, as expected. Denoting the distinct transfer matrix eigenvalues , we extract the correlation lengths using exponential fits, which are shown in (b) (using the same symbols), along with those extracted from the transfer matrix spectrum with or without flux inserted (shown as lines), with g the degeneracy of the eigenvalue. Both approaches agree for the spinon, vison and dimer correlation lengths which show no saturation with increasing χ.
as t a (a = 0, 1, ...) with |t 0 | > |t 1 | > |t 2 | > ..., it turns out t 0 is non-degenerate, suggesting that there is no longrange order in the variational wave function (confirming the ED results). The sub-leading eigenvalues t a (a = 1, 2, 3) are six-fold degenerate, followed by a non-degenerate t 4 . These eigenvalues give direct access to series of correlation lengths ξ (a) = −1/log(|t a /t 0 |), which therefore carry the same degeneracies. We have also computed the correlation length with a ±1 Z 3 flux by inserting a string of Z (or Z 2 ) operators, where the leading eigenvalue of the corresponding transfer matrix is denoted as t Z,1 [62]. From t Z,1 , which is nondegenerate, one can obtain the leading correlation length in the flux sector ξ (1) Z = −1/log(|t Z,1 /t 0 |). A summary of various correlation lengths versus χ from both methods is shown in Fig. 4(b). We find that the largest one in all sectors, ξ (1) Z , is equal to the correlation length found between a pair of visons; it is non-degenerate, in agreement with the fact that visons carry no spin. In the sector without flux, the leading correlation length ξ (1) is in perfect agreement with the correlation length ξ spinon extracted from placing a spinon-antispinon pair. Moreover, since PT symmetry maps spinons placed on reflected bonds to anti-spinons, we expect the spinon correlations to have a degeneracy structure 3 ⊕ 3 which indeed agrees with the six-fold degeneracy in ξ (1) , and is further supported by checking the U(1) quantum numbers of the t 1 multiplet. The U(1) quantum numbers further suggest that t 2,3 , which are also six-fold degenerate, also carry SU(3) representation 3 ⊕ 3. Thus, ξ (1,2,3) all correspond to spinon correlation length. This, in fact, is in correspondence with the three linear dispersing branches in the ES in the Q = ±1 charged sectors, as we shall discuss later. Further examining the correlation length, we find ξ (4) is identical to dimer correlation length, where non-degeneracy agrees with the fact that the dimer operator is invariant under SU(3) rotation. Depending on the parafermion type, the ξ parafermion have different values, both of which are smaller than the spinon correlation length. Interestingly, all of these correlation lengths, except the spin correlation length ξ spin , have no sign of saturation with increasing χ, in agreement with our expectation that the state is not in the Z 3 quantum double phase.
Degeneracy structure of topological chiral PEPS -A remarkable feature of our results is the correspondence between the leading four eigenvalues of the transfer matrix and the different sectors in the ES: The Q = 0 sector has one branch, while Q = ±1 each have three almost degenerate branches. This is in direct analogy to the unique leading eigenvalue t 0 which has trivial spin, and the approximate three-fold degeneracy of t 1 , t 2 , t 3 , which have perfectly degenerate spins 3 and 3, matching the perfect degeneracy between Q = ±1. A similar correspondence between (approximate) degeneracy of the (2D) transfer operator and of the ES branches was observed for chiral PEPS with SU(2) 1 counting, where it could be explained as arising from the symmetry of the tensors, and subsequently used to remove the degeneracy in the non-trivial sector in the vicinity of a (fine-tuned) perfectly degenerate point [59]. Furthermore, we checked that the same correspondence also holds in the PEPS description of non-Abelian SU(2) 2 CSL [21]. It is suggestive that such a correspondence in the (approximate) degeneracy structure is a general feature of chiral PEPS and will also hold for general SU(N ) k models; indeed, both ES and the eigenvalues t i are extracted from the same objects, namely the left/right (or up/down) fixed points of the CTMRG environment. It would be interesting to see whether such a correspondence holds in a context of general chiral models, and whether it could possibly even be used to further characterize the precise nature of a chiral theory.
Conclusion and outlook -In this work, we have proposed a model for a SU(3) 1 CSL on square lattice and unambiguously identified the relevant parameter space based on ED technique. Guided by ED and DMRG results, we have focused on constructing and optimizing a symmetric PEPS ansatz for the CSL, whose variational energy is remarkably good. For the first time, linear dispersing branches in all three sectors of SU(3) 1 WZW CFT can be obtained with PEPS. A comparison between edge spectrum and bulk correlations reveal a fine structure in the bulk-edge correspondence, which will be tested in further study of SU(N ) k PEPS CSL. Certain unresolved issues, e.g., the number of variationally degenerate ground states on torus, and anyon statistics of chiral topological order remain open, which we hope to uncover in the future.
Acknowledgements  on the other layer has the same spectrum as the one without flux, and the spectrum with Z or Z 2 on only one layer share the same absolute values but with different phase ω or ω 2 . To complement the main findings in the manuscript, we provide several relevant details in this supplementary material, organized as follows: basic knowledge about SU(3) group and its irreducible representations (irreps) in Sec. I, exact diagonalization study on various small tori in Sec. II, DMRG study on finite cylinders in Sec. III, construction of SU(3) symmetric PEPS and the tensor classification scheme in Sec. IV, the specific corner transfer matrix renormalization group (CTMRG) method we use and the optimization procedure in Sec. V, additional data for entanglement spectrum (ES) in Sec. VI, topological excitations and correlation functions of symmetric PEPS in Sec. VII, and in the end we also list the nonzero elements of the tensors in Sec. VIII.

I. BRIEF OVERVIEW OF SU(3) IRREPS
Since the theory of SU(3) group and its irreps can be found in many textbooks, e.g. Ref. 1, here we only list the relevant known results without derivation.
As a special case of the general representation theory of SU(N ) group, to each irrep of SU(3) we can associate a Young tableau containing a maximum of two rows (see Fig. S1). Denoting by p (q) the number of columns in the first (second) row, with p ≥ q, the dimension of the corresponding irrep is (1/2)(p + 2)(q + 1)(p − q + 1).

FIG. S1. A generic Young tableau characterizing an irrep of SU(3).
Unlike the SU(2) case where the states of a given multiplet are labeled by a unique U(1) quantum number (eigenvalue of S z ) and related to each other by a unique ladder operator S − (or S + ), multiplets of SU(3) should rather be seen as two-dimensional objects where states are characterized by two U(1) quantum numbers S z = (s z 1 , s z 2 ) and related by two ladder operators (S − 1 , S − 2 ) (or (S + 1 , S + 2 )). Note that a given couple (s z 1 , s z 2 ) is no longer necessarily associated to a unique state.
We note that, in irreps of SU(N ), there is some arbitrariness in defining the N − 1 diagonal generators (so-called Cartan subalgebra). Without basis change, it is possible to linearly combine them. The two U(1) quantum numbers shown in Tab. S1 indeed correspond to the eigenvalues of S 7 and The generator for the center of SU(3) can also be expressed in terms of the two diagonal generators as: The two-site permutation operator P ij used in defining the Hamiltonian, can be expressed with su(3) generators in the irrep 3 as:

II. EXACT DIAGONALIZATION ON VARIOUS CLUSTERS
For exact diagonalization, we have used the Lanczos (respectively Davidson) algorithm to compute the ground-state (respectively low-energy excitations) of our model on various finite-size clusters of N s sites with periodic boundary conditions, see Tab. S3. Since we are looking for a quantum spin liquid state, all clusters are adequate even though they possess different momenta in their Brillouin zone. Even more, we have considered some clusters which are not perfectly square to get additional signatures: we define the eccentric-ity of a cluster by the ratio of the two smallest inequivalent loops of the nearest-neighbor graph around the torus, which is a measure of the "two-dimensionality" of the cluster, where a value close to one is considered fully two-dimensional. In order to reduce the size of the Hilbert space, we have used all space symmetries (translation and point-group) as well as color conservation, which is equivalent to fixing the values of the two U(1) quantum numbers S z : namely, we diagonalize the Hamiltonian in a subspace with a given number of particles per color (N 1 , N 2 , N 3 ) with a constraint of singleoccupation, i.e.
Ns eccentricity point group 12 In the seminal work by Halperin [2], a SU(2) spin-singlet fractional quantum Hall (FQH) state was introduced for hardcore spin-1/2 bosons at filling ν = 2/3. As stated in the main text, the SU(3) 1 CSL that we are investigating is a lattice realization of such a phase [3]. The simplest signatures of a FQH phase are given by the ground-state degeneracy and the quasi-hole properties (quasi-degeneracy and momentum quantum numbers) [4,5]. Moreover, there is a simple generalized exclusion principle [6,7] with (2, 1) clustering properties in our case: for instance, in the spinful bosonic language, when N s = 3p, the three (quasi) degenerate ground-states are given by occupations: (↑, ↓, 0, ↑, ↓, 0, . . .) and its translations. These occupations have to be understood as a function of N s orbitals which are obtained when folding the Brillouin zone [4,5]. This exclusion rule simply enforces that there are no more than 2 particles in 3 consecutive orbitals and that a ↓ particle must necessarily be followed by a hole.
As a result, for all clusters N s = 3p ranging from N s = 12 to N s = 24, we have confirmed that our model does indeed show a quasi three-fold degeneracy of the ground-state, and their quantum numbers are given by the above generalized exclusion principle. This is a non-trivial feature and it is different from what would be expected for a charge-density wave phase. For example, in the cluster 18a, the three states are found at momentum Γ = (0, 0), while in the cluster 18b, one state is found at Γ and the two others at ±(2π/3, 2π/3) as predicted. For N s = 21, in all three considered clusters, we also FIG. S2. Low-energy spectra obtained from ED on additional Ns = 21 clusters: (a) 21b and (b) 21a. Each Brillouin zone is shown as inset. On these clusters, the ground-state is a global SU(3) singlet. We confirm the presence of two additional low-energy singlet states with the expected momenta, see text.
find three quasi-degenerate states with the correct momenta, see Fig. 1a in the main text and Fig. S2. Moreover, as shown in the main text, the ground-state energy shows a rather quick saturation with N s , compatible with a gapped phase.
Regarding the quasi-hole case, it can be obtained from clusters N s = 3p − 1. In such a case, the counting in the sector (N 1 , N 2 , N 3 ) = (p, p, p − 1) (which are three-fold degenerate due to SU(3) symmetry) is given by the generalized exclusion principle for spinful particles with N ↑ = p and N ↓ = p − 1 (number of holes being p). For example, for all clusters with N s = 20, we predict N s low-energy quasi-hole states, more precisely one (three-fold degenerate) per momentum sector using the heuristic rule [4,5], which is indeed observed in Fig. S3(a) (other data not shown), and each low-energy state transforms as a3 irrep.
For clusters N s = 3p−2, we expect that it would be best described by two quasi-hole states since a single quasi-hole excitation is rather localized. As a result, we would observe lowenergy excitations both in sectors (N 1 , N 2 , N 3 ) = (p, p, p−2) as well as (p, p − 1, p − 1) (and their equivalent sectors). This is indeed what we have found: for instance in both clusters N s = 19, the counting in sector (7, 7, 5) predicts 3 states per momentum and the one in sector (7, 6, 6) leads to 7 states per momentum, which is indeed observed in Fig. S3(b) (other data not shown). Reconstructing the SU(3) irreps from the states quantum numbers, one then finds four 3 and three 6 multiplets in each momentum sector, as predicted.

III. RELEVANT DETAILS OF DMRG METHOD
For DMRG, we have computed the ground-state wavefunction on various cylinders N s = L x × L y (with open/periodic boundary conditions in the long/short direction). We have used explicitly the two U(1) quantum numbers to ease convergence. Using up to m = 4 000 states, we can obtain reliable energies (discarded weight below 5e-5) up to L y = 6. In order to stabilize a global SU(3) singlet ground-state, we have chosen the system size N s multiple of 3, more specifically we took L x = 3p as the integer closest to 2L y .
By computing the total energy for cylinders L x × L y and (L x + 3) × L y , we can obtain an accurate estimate of the ground-state energy density (per site) by subtraction, providing the data plotted in Fig. 2b in the main article. Quite remarkably, there is a very fast convergence since all data for L y = 4, 5 or 6 are compatible with a ground-state energy density e 0 = −2.05(1), very close to the ED estimate.
In Fig. S4, we plot the bond strengths P ij on nearestneighbor bonds which do not show any modulation at all. Moreover, we have measured the local Cartan (s z 1 , s z 2 ) average values and found that they are vanishing (below 1e-6). All these measurements are indicative of a featureless phase.
Since our model possesses one SU(3) fermion per unit cell (equivalent to 1/3 filling), a trivial featureless gapped groundstate is impossible according to the Lieb-Schultz-Mattis theorem for SU(N) spin systems and its generalization to two dimensions [8][9][10][11]. Therefore, our ED and DMRG data are suggestive of a gapped topological phase.

IV. SU(3) SYMMETRIC PEPS ON THE SQUARE LATTICE
Here we present details about the symmetric PEPS construction, following the same spirit of Ref. 12-14. To construct a faithful PEPS representation of a chiral spin liquid wave function, we could encode the symmetry property of the desired wave function into local tensors. On the microscopic lattice scale, the symmetries that we need to take into account are: (1) the wave function |ψ is invariant under global SU(3) rotations, i.e., it is a SU(3) singlet; (2) under onesite translation and π/2 lattice rotation, |ψ is invariant up to a phase; (3) under lattice reflection P or time-reversal action T, |ψ is transformed into its complex conjugate |ψ → |ψ (also up to a possible phase), but is invariant under their combination PT. These symmetry requirements can be fulfilled by taking a suitable unit-cell of tensors, where these tensors satisfy certain symmetry constraints.
To implement global SU(3) symmetry, the PEPS wave function can take a form where a virtual SU(3) singlet formed by two virtual spins, denoted as |Ω , is put on every bond, and the on-site tensor P does a projection from four virtual spins on every site into the physical spin. Using translation symmetry, we can put the same virtual singlet on every bond, and the same local projector on every site, so as to work directly in the thermodynamic limit. The wave function then takes a simple form: where i stands for site index, and lrud, s is for left, right, up, down virtual, and physical spin on every site, respectively. The lattice point group symmetry imposes strong constraints on these local tensors, namely, up to a phase, the local projector P should be invariant under π/2 lattice rotation, and becomes complex conjugate under reflection, and the virtual singlet |Ω should be invariant under reflection. One noticeable difference between SU(3) (more generally SU(N )) and SU(2) group is that SU(2) group is self-conjugate such that one can form a singlet with two spins carrying the same irrep. Such a property is absent with general SU(3) spins, and we need to combine two spins carrying irreps with opposite Z 3 charges to form the singlet. The virtual space we use in this work, V = 3 ⊕ 3 ⊕ 1 with bond dimension D = 7, satisfies this SU(3) symmetry requirement, and allows us to construct the virtual singlet: where the labeling of basis for each irrep follows Tab. S1. This (unnormalized) maximally entangled virtual singlet is indeed symmetric under reflection. Unlike the SU(2) case, where using on-site unitary transformation on one sublattice, one can transform bond singlet into an identity matrix without changing the on-site projector, the same trick cannot be applied to SU(3) virtual singlet |Ω . Nevertheless, we can absorb the two neighboring bond matrices into the on-site projector P, e.g. the right and down one, forming the tensor A, without enlarging the unit cell. This strategy is taken in the numerical calculation in this work.
To systematically construct the on-site projector, we first did a classification of the rank-5 tensor. According to the fusion rules of SU (3) then determined. For each occupation number channel, the highest weight states (corresponding to S z = (1/2, 0)) are expressed in the tensor product basis of V ⊗4 . A point group analysis (C 4v ) is then performed (see Tab. S5), and highest weight states are symmetrized accordingly. Lower weight states (namely S z = (−1/2, 1/2) and S z = (0, −1/2)) are determined using lowering operators expressed in the tensor product basis of V ⊗4 .
As a result of the classification, the local projectors are now classified according to irreps of the square lattice point group C 4v , denoted as A 1 , A 2 , B 1 , B 2 and E. One can then construct the on-site projector P by linearly combining different classes of tensors, such that it is invariant under π/2 lattice rotation but becomes its complex conjugate upon reflection (up to an irrelevant phase). One choice we considered is: where {λ a 1 , λ b 2 } are real coefficients, as mentioned in the main text. Here N 1 = 6, N 2 = 5 is the number of tensors in the B 1 and B 2 classes, respectively. We note that, one could also use A 1 and A 2 classes to build chiral PEPS [14], whose energy turns out to be significantly higher than Eq. (S5) (data not shown). Thus we do not examine the detailed property of the later.
The expressions for the classes of tensor considered in this work are provided in Sec. VIII. See also Fig. S5 for a pictorial illustration.

V. CTMRG METHOD AND VARIATIONAL OPTIMIZATION
For completeness, here we briefly describe the specific CTMRG method we used in this work, which follows Ref. 15 and is further simplified in Ref. 16.
In our setting, for tensor network of the wave function norm, the unit cell contains only one tensor, denoted as E (see  Fig. S6(a)), which is obtained by contracting tensor A and its complex conjugate over the physical index. The CTMRG method allows us to approximately contract the whole network on the infinite plane by computing the effective environment tensors surrounding the unit cell. In our case, the environment tensors are composed by corner tensors and edge tensors {C i , T i }(i = 1, 2, 3, 4), see Fig. S6(b) for graphic notation. The accuracy of CTMRG method is controlled by the environment bond dimension, denoted as χ, and typically we choose χ = kD 2 (k ∈ N + ). In the CTMRG procedure, we dynamically increase χ by a small amount to keep the complete SU(3) multiplet structure. To further speed up the CTMRG procedure, we have explicitly kept track of the first U(1) quantum number of the SU(3) multiplets [17].
We note that, although the wave function has certain lattice symmetry, we do not use them in the CTMRG procedure, since after absorbing the bond singlet into on-site projector to construct tensor A, the tensor A is not invariant under π/2 lattice rotation. As a result, the four corner (edge) tensors {C i }(i = 1, 2, 3, 4) ({T i }(i = 1, 2, 3, 4)) are not necessarily the same. Nevertheless, we have checked that the physical observables, e.g., correlation lengths, along the horizontal and vertical directions are the same, as expected.
For a given set of variational parameters {λ a 1 , λ b 2 }, we can now compute the energy density with the environment tensors, simply by inserting the identity operator or local Hamiltonian terms in the central region, see Fig. S6(b). Energy gradient can then be easily obtained by finite difference method, which is feasible due to the significantly reduced number of variational parameters (compare to general PEPS ansatz). The conjugate-gradient method [18] is then utilized to find the variational optimal parameters. In practice, this optimization procedure is carried out with χ = D 2 . Then we evaluate the energy density of the optimized ansatz with several larger χ = kD 2 (k = 2, ..., 6) and eventually extrapolate to the χ → ∞ limit.

VI. ADDITIONAL DATA FOR ENTANGLEMENT SPECTRUM
The entanglement property of PEPS can be most easily characterized by studying the entanglement spectrum on finite width cylinders, which is defined to be minus log of the spectrum of reduced density matrix (RDM) of subsystem [19], say the left half of the cylinder. For PEPS on an infinitely long cylinder, the RDM can be constructed from the leading eigenvector of the transfer operator through the relation: where U is an isometry relating the physical degrees of freedom to the virtual ones [20], and we have adopted the convention that the first index of σ L,R is in the bra layer. This RDM further shares the same spectrum as ρ = σ L σ T R , which we diagonalize to get information about edge properties.
As mentioned in the main text, for our case with bond dimension D = 7, it is not feasible to compute σ L,R exactly, except for cylinders with small width. Instead, we can use the environment tensors computed from CTMRG method to approximate the σ L,R , see Fig. S7 for illustration. This is justified by the fact that CTMRG method is essentially approximating the fixed-point of certain transfer operator with matrix product state formed by environment tensors. One advantage of this approach is that, using σ L,R constructed in this way one can find RDM in all different charge sectors simultaneously, while with exact contraction one has to find them separately.
The RDM has both translation symmetry and SU(3) symmetry, allowing us to block diagonalize it with momentum quantum number K, the two U(1) quantum number, and Z 3 quantum number (charge Q). A typical result of full diagonalization is shown in Fig. S8, where linear dispersing chiral FIG. S7. Using CTMRG environment tensors to construct RDM for the left half. (a) shows the transfer operator on a width Nv = 6 cylinder, whose right leading eigenvector σR is approximated by a ring of T2 tensors, shown in (b). Similarly for the left leading eigenvector σL. The RDM is then obtained by contracting a ring of T4 and T2 tensors, shown in (c). modes can be seen in the low energy spectrum. The degeneracy between the charge Q = +1 and Q = −1 sector can be identified with Fig. S8(b) and (c), where same energy levels with same momenta but conjugated SU(3) irreps appear in the Q = +1 and Q = −1 sector separately, confirming the degeneracy mentioned in the main text. In all three Z 3 sectors, an entanglement gap [19] separating the chiral mode from the high energy continuum can be identified, although the magnitude of the gap in charge Q = 0 sector is much larger than the gap in charge Q = ±1 sectors. In the Q = 0 sector, the chiral branch starts at momentum K 0 = −π/3, while the three quasi-degenerate chiral branches in the Q = ±1 sectors start at momenta K ±1 = −π/3, π/3 and π, individually.
Further examining the level counting of the chiral modes confirms that they satisfy SU(3) 1 Wess-Zumino-Witten (WZW) conformal field theory (CFT) prediction [21,22]. See Tab. S6 for a list of the tower of states in SU(3) 1 WZW CFT. It is interesting to see that the Virasoro level contents also exhibit degeneracy between the charge Q = +1 and Q = −1 sectors, which is perfectly recovered by the numerically computed entanglement spectrum. Due to this degeneracy, in the following, we have plotted the Q = +1 sector and Q = −1 sector together, using open symbols and filled symbols respectively, to stress this symmetry (see Fig. S9 and Fig. S10).
It should be noted, when using CTMRG environment tensors to construct approximate RDM, the environment bond dimension χ is the only tuning parameter. Since χ controls the accuracy of CTMRG procedure, we expect that the level contents of chiral CFT mode becomes more complete with increasing χ. This finite χ effect is shown in Fig. S9. Certain features of the ES, e.g., momentum shift in all sectors and three branches in charged sectors, are present for all different χ we have considered. This is reasonable since the low energy spectrum converges first with increasing χ, and suggests that these features are intrinsic properties of the optimized PEPS wave function, i.e., not artifacts of the approximation.
The momentum shift observed in ES on N v = 6 cylinder has dramatic consequence for ES on N v = 3 and N v = 9 cylinders, shown in Fig. S10. In the charge Q = 0 sector for both N v = 3 and 9, a linear dispersing mode can be vaguely identified. However the content of each level is doubled, i.e., two singlets for n = 0 level, two 8 for n = 1  In each sector, the (in)complete Virasoro levels have been indicated by (blue) red boxes when necessary, and the missing levels can be found in the higher energy spectrum, marked as blue arrows. Their contents are shown in red vertically (see also Tab. S6). (Since the Q = +1 sector is degenerate with the Q = −1 sector, only Virasoro levels in the former sector are marked out.) With increasing χ, more Virasoro levels become complete, which is evidently seen in the charge Q = 0 sector. This trend is not monotonic in charged sectors, due to the large number of branches -three -on a relatively small width (Nv = 6) cylinder, which mix each other in the higher energy spectrum. Nevertheless, the CFT spectrum in all sectors get more separated from the high energy continuum with increasing χ.
level, since the finite momentum of ground state K 0 = −π/3 is incommensurate with N v = 3, 9. This scenario is further confirmed by exact contraction for N v = 3 case (data not shown). In the charged sectors with N v = 3, the low-est level can appear at different momenta, which typically depends on χ, see Fig. S10(c) and (d). This is in agreement with the three quasi-degenerate branches with different momenta Finally, we close this section by briefly discussing the ES in the flux sector. In this work we have mainly focused on ES in the topological sectors without flux, and succeeded in finding signatures for chiral topological order. The Z 3 gauge symmetry implies that we can also construct topological sectors with flux, see Fig. S11 for illustration. However, previous studies in the SU(2) case [23,24] suggest the ES in flux sectors does not follow a simple CFT description. Therefore, we do not

l a t e x i t s h a 1 _ b a s e 6 4 = " 6 + D Q A A 1 b J + C n G D + 2 7 5 D r v L l m U z g = " > A A A C z n i c j V H L T s J A F D 3 U F + I L d e m m g Z i Y m J C i x s e O 4 M Y l G n k k Q E h b B m w o b T O d k h A k b v 0 A 3 e p n G f 5 A / 8 I 7 Q z E a Y v Q 2 b c + c e 8 + Z u X e s w H V C Y R i T h L a w u L S 8 k l x N r a 1 v b G 6 l t 3 c q o R 9 x m 5 V t 3 / V 5 z T J D 5 j o e K w t H u K w W c G b 2 L Z d V r d 6 l z F c H j I e O 7 9 2 K Y c C a f b P r O R 3 H N g V R 9 f t i 6 6 b B T a / r s l Y 6 a + Q M F f o 8 y M c g W 8 g 0 D p 8 m h W H J T 7 + h g T Z 8 2 I j Q B 4 M H Q d i F i Z C e O v I w E B D X x I g 4 T s h R e Y Y x U q S N q I p R h U l s j 7 5 d W t V j 1 q O 1 9 A y V 2 q Z d X H o 5 K X X s k 8 a n O k 5 Y 7 q a r f K S c J f u b 9 0 h 5 y r M N 6 W / F X n 1 i B e 6 I / U s 3 q / y v T v Y i 0 M G 5 6 s G h n g L F y O 7 s 2 C V S U 5 E n 1 7 9 1 J c g h I E 7 i N u U 5 Y V s p Z 3 P W l S Z U v c v Z m i r / r i o l K 9 d 2 X B v h Q 5 5 S X f C F j N O v 6 5 w H l a N c / j h 3 c p 3 P F o q Y R h J 7 y O C A 7 v M M B V y h h L K a + D N e 8 K q V t I E 2 1 h 6 m p V o i 1 u z i R 2 i P n 9 P 6 l v E = < / l a t e x i t > hB L | < l a t e x i t s h a 1 _ b a s e 6 4 = " S i a f W W A V 9 N 3 U m A J f A v F Z Z X w u 3 8 4 = " > A A A C z 3 i c j V H L S s N A F D 2 N r 1 p f V Z d u Q o s g C C V V 8 b E r d e P C R Q v 2 A W 0 p S T q t w W k S k o l S q u L W v b j V v 5 L + g f 6 F d 6 a p K E X 0 h i R n z r 3 n z N y 5 l s + d U B j G K K H N z M 7 N L y Q X U 0 v L K 6 t r 6 f W N a u h F g c 0 q t s e 9 o G 6 Z I e O O y y r C E Z z V / Y C Z f Y u z m n V 1 K v O 1 a x a E j u d e i I H P W n 2 z 5 z p d x z Y F U c 0 m N 9 0 e Z 3 q x f X 7 b T m e N n K F C n w b 5 G G Q L m e b u 0 6 g w K H n p N z T R g Q c b E f p g c C E I c 5 g I 6 W k g D w M + c S 0 M i Q s I O S r P c I c U a S O q Y l R h E n t F 3 x 6 t G j H r 0 l p 6 h k p t 0 y 6 c 3 o C U O r Z J 4 1 F d Q F j u p q t 8 p J w l + 5 v 3 U H n K s w 3 o b 8 V e f W I F L o n 9 S z e p / K 9 O 9 i L Q x b H q w a G e f M X I 7 u z Y J V K 3 I k + u f + t K k I N P n M Q d y g e E b a W c 3 L O u N K H q X d 6 t q f L v q l K y c m 3 H t R E + 5 C n V g E 9 k H H 6 N c x p U 9 3 L 5 / d x B O Z 8 t F D G O J L a Q w Q 7 N 8 w g F n K G E C n n 7 e M Y L X r W y d q P d a w / j U i 0 R a z b x I 7 T H T y 0 6 l w 8 = < / l a t e x i t >
FIG. S11. On infinitely long cylinders, topologically quasidegenerate states can be constructed by choosing virtual boundary state |BL , |BR belonging to fixed Z3 charge sector, with or without nontrivial Z3 flux insertion (shown as blue squares). explore it here but leave it to further study.

VII. TOPOLOGICAL EXCITATIONS AND CORRELATIONS IN SYMMETRIC PEPS
As discussed in the main text, Z 3 gauge symmetry, generated by Z(Z 3 = I D ), implies topological excitations on infinite plane, whose type can be labeled by the group element and group irreps. Here we will not describe the full details of the theory, but refer to Ref. [25] for the interested reader. Spinon, one of the topologically nontrivial elementary excitations, can be created by modifying a local tensor such that it belongs to different irreps of Z 3 . This can be achieved by acting on the virtual level of the local tensor A with an operator X, see Fig. S12. Its anti-particle can then be similarly created with an operator X 2 also acting on the virtual index.
Apart from basic algebraic relation of X and Z: with ω = e i2π/3 , which generalizes the anti-commutation relation between Pauli matrix σ x and σ z to Z 3 , the choice of X is not unique, due to the internal SU(3) symmetry in the ansatz. The specific X we use to compute spinon-antispinon correlation function is: In principle, the operator X can be put on any of the four virtual indices of A. However, unlike the Pauli matrix σ x in the Z 2 case, X in the Z 3 case cannot be chosen to be symmetric. Thus the order of indices matters when computing anyonic correlation functions using X. In practice, we always put X or X 2 in the ket layer, with the first index contracted with down index of local tensor A. Apart from spinons, the Z 3 gauge symmetry also allows us to construct vison excitations, which are end points of Z or Z 2 string operator acting on the virtual level. Their bound states, so-called parafermions, can be created by attaching spinons to the end points of virtual string.
With these topological excitations at hand, the calculation of their correlation functions is straightforward, shown in Fig. S13(a) and (b). The corresponding transfer matrix can also be constructed, without or with flux, see Fig. S13(c) and (d).

VIII. NONZERO ELEMENTS OF SYMMETRIC TENSORS
For the sake of completeness, we now present the resulting tensors from classification. According to the irreps of C 4v group, the on-site projectors can be classified into four real classes, and two complex classes. Since only the real classes are used, we list their nonzero elements below, denoted as A 1 , . See also Tab. S1 for the U(1) quantum numbers of each basis in both physical and virtual spaces.
The expressions of the three components of each real tensors are provided in tables S17-S27, where the tensor indices are in [up, left, down, right] order.