Dimer description of the SU(4) antiferromagnet on the triangular lattice

In systems with many local degrees of freedom, high-symmetry points in the phase diagram can provide an important starting point for the investigation of their properties throughout the phase diagram. In systems with both spin and orbital (or valley) degrees of freedom such a starting point gives rise to SU(4)-symmetric models. Here we consider SU(4)-symmetric “spin” models, corresponding to Mott phases at half-filling, i.e. the six-dimensional representation of SU(4). This may be relevant to twisted multilayer graphene. In particular, we study the SU(4) antiferromagnetic “Heisenberg” model on the triangular lattice, both in the classical limit and in the quantum regime. Carrying out a numerical study using the density matrix renormalization group (DMRG), we argue that the ground state is non-magnetic. We then derive a dimer expansion of the SU(4) spin model. An exact diagonalization (ED) study of the effective dimer model suggests that the ground state breaks translation invariance, forming a valence bond solid (VBS) with a 12-site unit cell. Finally, we consider the effect of SU(4)-symmetry breaking interactions due to Hund’s coupling, and argue for a possible phase transition between a VBS and a magnetically ordered state. Copyright A. Keselman et al. This work is licensed under the Creative Commons Attribution 4.0 International License. Published by the SciPost Foundation. Received 18-11-2019 Accepted 21-04-2020 Published 13-05-2020 Check for updates doi:10.21468/SciPostPhys.8.5.076


Introduction
Frustrated quantum antiferromagnets may possess non-magnetic ground states, avoiding spin order through short or long range entanglement of spins. In the late 1980s and early 1990s, a dominant approach to this physics was based on generalizing the SU(2) group of spin rotations to SU(N) or Sp(2N) [1][2][3][4]. In the limit N → ∞, models with such enlarged symmetry may be solved exactly by a fully symmetric saddle point of a path integral representation of the partition function, and consequently possess non-magnetic ground states. More recently, it has become possible to study models with SU(N) symmetry for finite N using computational methods. In particular, models describing SU(N) fermions at half filling were shown to host a variety of non-magnetic states [5][6][7][8][9][10], indicating that features of the N → ∞ solutions survive to N of order one.
Interest in models of this type has also been stimulated by their possible experimental realization in cold atoms [11,12], in materials with orbital degeneracy [13,14], and, more recently, in moiré superlattices with valley degeneracy [15][16][17][18]. Here we investigate in particular the SU(4) antiferromagnet in the self-conjugate six dimensional representation. This is the rep-resentation corresponding to two electrons distributed amongst four degenerate spin/orbital states on each site. This representation thus occurs naturally in systems with spin and a twofold orbital or valley degeneracy. Like the familiar S=1/2 SU(2) spins, pairs of spins in this representation may form a singlet "valence bond", so that there is a natural "dimer" picture upon which non-magnetic states may be based.
Indeed, it has been shown by Rokhsar [4] that at N = ∞, dimerized states (i.e. products of singlet bonds) are the ground states in the self-conjugate representation for a very wide class of lattices and exchange interactions (including almost all those of physical interest). For N large but finite, it is therefore expected that a quantum dimer model [19][20][21][22][23], which describes the projection of the Hamiltonian to the singlet subspace, should capture the physics of the problem. A dimer model consists of a Hilbert space in which dimers, which stand in for singlets, realize a covering of the lattice, with each site covered by one, and only one dimer. Such models are known to predict quantum spin liquids and various valence bond solid orders, depending on the exact model, dimensionality, and lattice. Notably, the simplest quantum dimer model on the triangular lattice was argued to possess a 2 spin liquid ground state [22].
With these motivations, we study the aforementioned SU(4) model on the triangular lattice both analytically and numerically. In contrast to the previous numerical works mentioned above, which focused on bipartite lattices, the triangular lattice considered here does not to our knowledge admit a sign-free Monte Carlo approach for generic parameters (though for a specific choice of parameters a sign-free algorithm exists [24]). Hence we attack the problem differently, using the density matrix renormalization group (DMRG), exact diagonalization, and an analytic dimer expansion. First, we describe a reformulation which takes advantage of the fact that SU(4) is a double cover of SO (6), through a mapping of the self-conjugate representation of SU (4) to the vector representation of SO (6), which we define. We determine the classical ground states of the model, and find that they generalize the three-sublattice coplanar states of the S=1/2 Heisenberg model on the triangular lattice. However, carrying out a numerical study of the model using the DMRG method, we argue that this order is absent in the quantum limit. We note that this result is in agreement with a recent pseudofermion functional renormalization group study [25] of the same model. We then extend the "overlap expansion" approach -originally developed by Rokhsar and Kivelson [23] to obtain the quantum dimer model Hamiltonian from the SU(2) spin Hamiltonian -to derive an analytic overlap expansion for the quantum dimer model relevant to the SU(4) case. We show that the parameter for the overlap expansion is x = 1/6 in the SU(4) case, which should be compared to x = 1/2 for SU (2) spins: hence the expansion is expected to be much more accurate for the present problem. Our exact diagonalization (ED) study of the effective dimer model suggests that the ground state is a twelve-site valence bond solid (VBS), although larger system sizes are required to conclusively rule out a proximate quantum spin liquid state. We remark that the related model of Ref. [24] shows such VBS order when generalized from SO (6) to SO(N) with N>12, which provides some further indications for this VBS state in this family of model Hamiltonians.
In addition, we study a generalization of the model which breaks SU(4) symmetry down to SU(2)×SU(2), by including an additional "atomic" Hund's coupling J H on the sites. In the limit of large J H , the model reduces to the SU(2) symmetric spin S=1 Heisenberg Hamiltonian, and the ground state has long-range three sublattice order. The generalized model thus exhibits a quantum phase transition from a paramagnetic to antiferromagnetically ordered state at zero temperature by varying J H . We present signs for this transition in numerics.

From Hubbard to Heisenberg
In the limit of strong interactions, electrons localize and the appropriate description of their physics becomes that of their spins, localized at lattice sites. Relevant models may be derived from Hubbard models. Here we proceed with this approach and consider two electrons per site hopping on the triangular lattice with hopping parameter t and subjected to an on-site Hubbard interaction U as well as a Hund's coupling J H , which tends to enforce an alignment of the spin degrees of freedom and thus breaks SU(4) symmetry. More precisely, we consider the following Hamiltonian: where c † i,a , c i,a are flavor a electron creation and annihilation operators at site i, In the final term, we used a notation based on the physical origin of the four electronic states from the spin-1/2 of the electron s z = ±1/2 and a two-fold orbital degeneracy τ z = ±1/2. We may associate a = 1, 2, 3, 4 (note the use of sans-serif font for this purpose) to the states with (2s z , 2τ z ) = (1, 1), (−1, 1), (1, −1), (−1, −1), respectively. The σ µ = σ µ ⊗ 2 (µ = x, y, z) Pauli matrices act in spin space. We specialize to the case of filling n i = 2, and carry out the standard degenerate perturbation theory in t/U to derive an effective "spin" model. A basis for the states with n i = 2 on a single site is: where, on the right-hand-side of the equations, |a b〉 = c † a c † b |0〉, where |0〉 is the vacuum. We start by analyzing the SU(4) symmetric model, i.e. we take J H = 0, and return to the effects of a finite J H later in Sec. 3.3.2. At second order in small t/U, Eq. (1) becomeŝ with J ∼ t 2 /U andT a b = P 6 c † a c b P 6 , with P 6 the projection operator onto the six-dimensional vector space defined above. TheT a b are 6×6 matrices which are related to the generators of SU(4). In Eq. (3) we extracted a constant which sets the zero of energy at a convenient value. To see this, we bring out the analogy to SU(2) spins by extracting the trace from thê T a b matrices: With this definition TrT a b = 0 (note that the trace here is over the 6-dimensional SU(4) space). The Hamiltonian in Eq. (3) can now be written aŝ One can also check that aT aa = 0, so that there are clearly only 15 such independent SU (4) matrices, which comprise a basis for the generators of SU(4) in the 6-dimensional representation. Now, regardless of the precise microscopic Hamiltonian H, we may consider a spin model, determined on the basis of, and constrained by symmetry. To proceed to the derivation of the most general SU(4) model, it is useful to make use of the following.

Map to the vector representation of SO(6)
SU(4) is a double cover of SO (6) and there exists a convenient map (which is faithful) from the six-dimensional representation of SU(4) to the fundamental (vector) representation of SO(6) [26]. By using the following basis, where each basis state transforms under the vector representation of SO(6), i.e. O : |n〉 → O|n〉, where O is an SO(6) matrix, it is straightforward to write all the SU(4) invariant two-site operators: HereP i j is the singlet projector over sites i j, where a (normalized) singlet over sites i j is written whileΠ i j is the permutation operator over sites i j andÎd i j is the identity. Then the general SO(6) invariant Hamiltonian with nearest-neighbor interactions is a sum of these terms: for α, β, γ ∈ . The "Heisenberg" model Eq. (3) is realized for −α = β = J, γ = 0. One can readily check then that for two sites with the SU(4) (or SO(6)) singlet in Eq. (10), In this SO(6) basis, we may also define the symmetricŜ mn and antisymmetric operatorŝ A mn , as well as the Hermitian (and still traceless), versions of the latter,Â mn = iÂ mn A mn = iÂ mn .
TheÂ mn operators can be considered as the generators of SO (6), and their square, Tr[Â ·Â T ] = 5Îd 6 is the quadratic Casimir operator, up to a normalization constant. Here we have defined the matrix of operatorsÂ such that (Â) mn =Â mn . Using these operators, the "Heisenberg" HamiltonianĤ becomeŝ where the trace, · and transpose operations act on the superscripts of theÂ l matrices ofÂ mn l operators.

Magnetic order
In this section, we first examine the classical ground states of the SU(4) model, which are "magnetically" ordered, i.e. they break the SU(4) symmetry and have a non-zero expectation value of the "spin" operator matrixT ab i orÂ mn i on each site. Having identified the type of magnetic order which is most favored, we next describe numerical studies which search for it. We find that this magnetic order is in fact absent, and that the ground state appears to lack any form of SU(4) symmetry breaking, i.e. is non-magnetic.

Classical limit
In order to look for a product ground state we first ask about the definition of the classical limit of SU(4) (SO(6)) spins. Like for SU(2) spins, we should replace, in the Hamiltonian, each of the fifteen SO(6) generatorsÂ mn , which are 6×6 matrices, by a single classical number.
To do so, we interpret the classical limit as a variational problem in the subspace of states consisting of direct products of single-site wavefunctions. Within any such state, the expectation value of any product ofÂ mn is replaced by a product of expectation values of eachÂ mn , which are c-numbers, as desired. A general single-site wavefunction is given by where v i is a complex six-dimensional unit vector, i.e. p |v i p | 2 = 1, so that |ψ〉 i is normalized. Upon going to the classical limit, The matrix A is now an antisymmetric 6 × 6 matrix of scalar matrix elements A mn , and the Heisenberg Hamiltonian becomeŝ Note that, while each operatorÂ mn for fixed m, n is Hermitian, the matrix A is real and antisymmetric. Moreover, while 6 m,n=1Â mnÂmn = 5Îd 6 , one can show that 0 ≤ TrAA T ≤ 1 (see Appendix A). In solving a classical SO (6) in this representation model, one should find matrices A which verify the above constraints (much like SU(2) S=1/2 (resp. S=1) classical spins are described by a three-dimensional vector with unit norm |S| = 1 (resp. with 0 ≤ |S| ≤ 1).

Product variational states
We now specialize to the triangular lattice and nearest-neighbor "Heisenberg" Hamiltonian, and look for the ground state of the corresponding classical model. The SU(2)-invariant spin-1/2 model on the triangular lattice is one of the best-studied models of frustrated magnetism. Its ground states are the so-called 120 • -ordered states. They triple the unit cell and each elementary triangular unit is such that the spins point at 120 degrees of one another.
We may rewrite the classical version of the Heisenberg Hamiltonian in Eq. (18), as where the sum runs over all unit triangles of the triangular lattice. The energy is clearly minimized when the first term in the square brackets vanishes and the second one is maximized.
The magnitude of the second term is maximized when the upper bound on TrA i A T i is saturated for all i ∈ t. As shown in Appendix A, the upper bound is saturated when the complex vector v describing the single-site wavefunction is given by v = (x + iy)/ 2, with x, y real, six-dimensional orthogonal unit vectors, i.e. |x| = |y| = 1 and x · y = 0. The first term vanishes for three-sublattice states that satisfy A 1 + A 2 + A 3 = 0 6 .
To minimize the first term in the square brackets in Eq. (19), and simultaneously maximize the magnitude of the second term, we can choose corresponding to where x, y, z are three orthonormal unit vectors such that x · y = x · z = y · z = 0. For x, y, z along each of the first three basis vectors of 6 , we get for example Note that once the spins on two nearest-neighbor sites are fixed, the remainder are fully determined by the condition A 1 + A 2 + A 3 = 0 6 , which can be successively applied to the spins on triangles sharing two of the sites which have already been fixed, to cover the entire lattice. This implies that all classical ground states are of the three-sublattice type.

Numerical analysis using DMRG
To probe the presence of magnetic order in the system we study the model numerically, using DMRG [27,28]. To this end, we consider finite cylinders in a geometry that allows for the formation of a 120 • -ordered state. Denoting the basis vectors of the triangular lattice by a 1 = (1, 0), a 2 = (1/2, 3/2), we consider cylinders such that the sites of the lattice modulo R = N y a 2 are identified, and N y is a multiple of three. This geometry is depicted in Fig. 1(c,d). Due to the large single-site Hilbert space dimension in this problem we are limited to narrow cylinders with N y = 3 and N y = 6, and we only study very short cylinders for the latter. We note that a different geometry, namely one in which lattice sites modulo R = N y ( a 2 − a 1 /2) are identified, is also compatible with a 120 • -ordered state for N y = 4. However, in this case, we find indications that the system behaves as a quasi-1D system with localized modes at the ends of the cylinder. We thus leave out these results from the discussion of the 2D limit presented here.
In the following, we first discuss the flavor gap in the system, and show that it remains finite, suggesting the absence of a low-energy Goldstone mode that would be expected if the system formed a magnetically ordered state. We then probe the presence of long range order in the ground state by looking at the static response of the system to polarizing fields applied at its boundary. We show that the expectation value of the magnetization decays rapidly away from the boundary in the presence of SU(4) symmetry, implying a lack of long range order.
Our DMRG simulations were performed using the ITensor library [29].

Flavor gap
We calculate the flavor gap only for cylinders of width N y = 3, as extracting the gap requires finite length scaling, and we are limited to very short systems for cylinders of width N y = 6, as mentioned above. Note that for N y = 3, the length of the system N x has to be even to allow for an SU(4)-singlet ground state.
Similarly to the conservation of s z which is often used when studying SU(2) spins, we can employ the conservation of three U(1) quantum numbers for the SU(4) case: t 3 ≡ n 1 − n 2 , t 8 ≡ n 1 + n 2 − 2n 3 , and t 15 ≡ n 1 + n 2 + n 3 − 3n 4 , where n a=1,..,4 denote the occupations of the four flavors as before. We calculate the gap of a ∆ t 3 = 2 excitation. To this end, we first obtain the ground state, which we expect to be an SU(4) singlet, and hence lie in the (t 3 , t 8 , t 15 ) = (0, 0, 0) sector, and then calculate the lowest energy state in the (t 3 , t 8 , t 15 ) = (2, 0, 0) sector. The latter state belongs to the 15 dimensional irreducible representation of SU(4), as can be verified by calculating the quadratic Casimir operator abT ba . The resulting energy gap is plotted in Fig. 1(a) as function of inverse system length. Even though we present data for relatively short systems, it is clear that the gap remains finite in the infinite system size limit, and we can estimate it to be larger than 2.5J. The maximal bond dimension in our simulations was M = 4000, resulting in a truncation error of 10 −5 (10 −4 ) for the largest system size in the t 3 = 0 (t 3 = 2) sector.

Probing long range magnetic order in the ground state
For the analysis of long range magnetic order in the ground state, it is instructive to consider the effect of SU(4) symmetry breaking by a Hund's coupling term as introduced in Eq. (1). More specifically, the Hamiltonian we consider iŝ As was mentioned previously, a finite J H > 0 breaks the SU(4) symmetry down to SU(2)×SU (2), pairing the two electrons on each site into a spin-triplet state. Using the definitions Eq. (2,6), the projection on the on-site spin-triplet and spin-singlet subspaces is given bŷ so that the Hund's coupling term can be written simply asŜ 2 i = 2P i,S=1 . In the large J H limit, the spin model Eq. (22) reduces to an SU(2) spin-1 Heisenberg model, for further details). The latter is known to form a 120 • -ordered state on the triangular lattice [30]. Below, we study the model in Eq. (22) as J H is increased from J H = 0 (the SU(4) symmetric point), where a three-sublattice order is predicted by our classical analysis, to a large J H J, where the 120 • order is known to form also in the quantum limit.
Once SU(4) symmetry is broken down to SU(2)×SU(2), only two U(1) quantum numbers are conserved: the z components of the spin and valley degrees of freedom, namely 2s z = n 1 + n 3 − n 2 − n 4 and 2τ z = n 1 + n 2 − n 3 − n 4 . Employing the conservation of these two quantum numbers, we now look for the ground state in the sector (s z , τ z ) = (0, 0). Once again, we consider finite cylinders of geometry and size compatible with the 3-sublattice order of the 120 • state. In particular, we consider cylinders of width N y = 3 and N y = 6 in the same geometry as before.
To facilitate the formation of a long range ordered state, we follow the approach introduced in Ref. [31] and apply pinning fields at the boundaries of the cylinder. We then calculate the expectation value of the spin component parallel to the field in the bulk, far from the boundary for different ratios of the length of the cylinder to its circumference. A complementary analysis, where we calculate the spin-spin correlations in the absence of pinning fields is presented in Appendix C.1 and gives similar results. To retain the conservation of s z , we apply the pinning fields only along the s z axis. More specifically, the field applied is −s z on the A sublattice, and +s z /2 on the B and C sublattices as depicted in Fig. 1(c,d). Note that in the SU(4)-symmetric case this corresponds to a field along A = A 21 (see Sec. 3.2).
The expectation value of s z on a site in the middle of the system, as function of J H , for different system sizes is shown in Fig. 1(b). The expectation value of s z remains small close to J H = 0, even when the ratio of the length of the cylinder to its circumference is unity, suggesting the absence of magnetic order in this case. As J H is increased, a finite expectation value develops as expected. The range of system sizes accessible by our simulations is not large enough to perform finite size scaling, but a relatively sharp increase in the magnetization around J H /J 1 suggests a phase transition occurs in the vicinity of this value. In Figs. 1(c,d) we plot the expectation values of s z on all the sites of a 6 × 6 cylinder, for J H = 0 and J H /J = 2 respectively. While in the former case, the magnetization decays rapidly away from the boundary where the pinning fields are applied, in the latter case the magnetization is finite and uniform across the system.
In these simulations the maximal bond dimension for cylinders of width N y = 3 was M = 2000, resulting in a truncation error smaller than 5 · 10 −5 . For cylinders of width N y = 6 the maximal bond dimension was M = 8000 for J H = 0 and M = 4000 for J H > 0, resulting in a truncation error of ∼ 2 · 10 −3 for values of J H /J < 1.5 at which no long-range ordering is observed, and a truncation error of ∼ 5 · 10 −4 or smaller for J H /J ≥ 1.5 at which a 120 • order develops.
To summarize, our numerical study suggests that the SU(4)-symmetric Heisenberg model does not have magnetic long-range order. A transition into a 120 • -ordered state can be driven by a Hund's coupling term which breaks SU(4) symmetry.

Singlet projection
The short-range nature of the spin correlations observed in DMRG motivates an approach focusing on SU(4) singlets. As we saw explicitly in Eq. (10), the six-dimensional representation of SU(4) considered in this work allows for the formation of a singlet on a pair of sites. Hence we can build many singlet states for the entire system by partitioning the sites into pairs, and placing each pair of corresponding spins into a singlet state. Following the pioneering work of Rokhsar and Kivelson [23] who considered the projection of the usual SU(2) Heisenberg model (in the S=1/2 representation) to a nearest-neighbor singlet manifold, we study the projection of the SU(4) Heisenberg model onto the subspace of nearest-neighbor SU(4) singlet "dimer" coverings of the lattice.
In this section, we start with a simple analytic comparison between energies of the singlet states and those of the classical ones discussed earlier, showing that the dimer states are superior in a variational sense. Then we provide further numerical justification for the projection to the singlet subspace. We next discuss the projection of the SU(4) Hamiltonian to the nearest-neighbor singlet coverings subspace and derive an effective dimer model. Finally we study the resulting dimer model using exact diagonalization.

Crude estimate of energy competition between singlet and ordered states
We first estimate the energy of such a singlet state, and compare to that of an ordered state. The optimal ordered product states were found in Sec. 3.2. They comprise 3-sublattice ordered states which spontaneously break SU(4) symmetry analogously to the 120 • ordered states for classical SU(2) Heisenberg spins. In those states, the energy per bond is the same for all bonds and is equal to Hence Now we consider a singlet state which is the product of two-site singlet "dimers". Specifically, a singlet covering is given by a partition of the set of N sites i into pairs C = {(i 1 j 1 ), (i 2 j 2 ), · · · (i N /2 , j N /2 )}, where (i, j) denotes a pair of nearest-neighbor sites. Such a state can be visualized by drawing a dimer -a colored bond -between the pairs of sites (i a j a ). We define using normalized singlets |s〉 i j as in Eq. (10). Note that in contrast to the SU(2) case, in the SO(6) representation the singlet state has a purely positive wavefunction, and is without any sign ambiguity. Thus there is no need to define the directionality of a singlet which is required to determine the sign of the wavefunction in the SU(2) case. For a crude estimate, we consider the variational energy of a single dimer covering, Unlike for the classical state, all the bond expectation values are not equal. As shown in Eq. (12), the singlet |s〉 i j is an eigenstate of H i j , with energy −5J. Hence 〈C|H i j |C〉 = −5J for those bonds covered by dimers. For bonds that are not covered by singlets, the two spins on the bond are uncorrelated, and one has 〈C|H i j |C〉 = 0 for those bonds. Hence the variational energy of the dimer state is −5J per bond times the fraction of bonds occupied by singlets, which is 1/6. Thus, since N bonds = 3N sites Comparing Eq. (28) and Eq. (25), we see that the dimer state has lower energy. This gives some simple understanding of the avoidance of magnetic order.
It is instructive to compare to the SU(2) case, with spin S spins. In this case for the usual Heisenberg model the classical product ground state with 120 • order has J〈S i ·S j 〉 class = −JS 2 /2, so the classical energy is For a spin singlet bond, we can write , so that 〈JS i · S j 〉 singlet = −JS(S + 1). Thus the dimer energy is Comparing Eq. (29) and Eq. (30), we see that the energies are equal for S = 1/2 (E = −J N bonds /8), with the classical state superior for all larger S. In summary the simplest possible variational dimer state of a single singlet covering is already better than a classically ordered state for the SU(4) problem, which is distinctly different from the SU(2) case. In the following sections we will refine the approach to singlet states, and consider superpositions of many terms, each with the form of Eq. (26).

Numerical justifications for the projection onto the singlets subspace
We define a nearest-neighbor singlet subspace as the Hilbert space spanned by superpositions of all nearest-neighbor singlet coverings of the form of Eq. (26). In this subsection we compare the low energy spectrum of the Hamiltonian in the full Hilbert space with that of its projection onto this nearest-neighbor singlet space. To this end we perform a numerical study on systems with size of up to 18 sites, using ED for systems with less than 12 sites, and Matrix Product State (MPS)-based simulations for larger systems, as described in detail in Appendix C.2.
We first compare the flavor gap in the SU(4) spin model with the gap in the projected problem. For cylinder of width N y = 3, the flavor gap was discussed in Sec. 3.3 and estimated to be larger than 2.5J for an infinitely long cylinder. The gap obtained for the projected problem is 0.281J, 0.203J for system sizes of N x = 4 and N x = 6 respectively. For cylinders of width N y = 4, in the same geometry, we find the flavor gap to be very weakly dependent on system size already for small system sizes, and larger than 3.8J. The gap in the projected problem is Table 1: Wavefunction overlaps between the ground state of the SU(4) spin model and the ground state of the Hamiltonian projected onto the subspace of nearest-neighbor singlet coverings for different system sizes (N x × N y ). "OBC" indicates open boundary conditions along both x and y, while "Cylinder" indicates periodic boundary conditions along y and open boundary conditions along x. The error indicated in brackets is estimated from the DMRG truncation error for the ground state of the spin model. Values for which no error is indicated were obtained using ED.  (1) 1.738J, 1.724J for system sizes of N x = 3 and N x = 4 respectively. Thus, we find that in both cases the gap of the projected Hamiltonian is smaller than the flavor gap, suggesting that the low energy physics is governed by the singlets.
In addition, we calculate the overlaps between the ground state of the SU(4) spin model and that of the projected Hamiltonian. These are summarized in Table 1 for a number of system sizes and different boundary conditions. We find that the overlaps decrease with increasing system size as expected. However, given the immense reduction in the dimension of the Hilbert space upon the projection, we find surprisingly large overlaps even for systems with N 10 − 20 sites.

Derivation of the effective dimer model
We now turn to the analytic derivation of the projected Hamiltonian. Rokhsar and Kivelson [23] constructed an expansion to express the effective projected Hamiltonian as a sum of local terms of increasing length of dimer re-arrangements. We obtain a similar expansion here for the SU(4) ∼ SO(6) case. We follow specifically a reformulation of the expansion by Ralko et al. [20].
We seek the best variational state of the form The wavefunction ψ C is required to minimize where The minimum of the variational energy is given by the condition ∂ E/∂ ψ * = 0. This gives where E 0 = min ψ E(ψ) is the best variational energy. This is a generalized eigenvalue problem for E 0 . We can convert it to a conventional one by defining Ψ = S 1/2 ψ, which leads to Therefore the variational ground state energy (and from it ultimately the variational ground state wavefunction) is obtained from the ground state of H eff , which is the desired effective quantum dimer Hamiltonian.
To obtain H eff , we expand both H and S in a series of increasingly small terms, which are related to the number of dimer rearrangements forming "loops". The small parameter of this expansion is the overlap x in the smallest such non-trivial loop: two dimers cyclically permuted on four sites. More generally, the inner product of a sequence of dimers pairing sites x = 1/6. In a full calculation of H and S, products of such overlaps appear, resulting in multiple factors of x, which determines the order of these terms in the expansion. Details of this quite technical procedure, which we formulate for SU(4) on a general lattice, will be presented in a separate publication. Starting from the general SO (6) where the prime on the sum indicates a sum over all the symmetry-equivalent plaquettes shown in the bras and kets, throughout the lattice. All the coefficients are given in terms of α, β, and x and are summarized in Table 2. Previous studies of this model on the triangular lattice [21,22,32,33] found that the ground state for v/t = 0 is a 12 × 12 VBS state. At large enough negative v/t, the ground state is a columnar ordered

Numerical study of the dimer model
state, while for positive v/t, first a phase transition into the 2 RVB spin liquid phase occurs at v/t 0.83, followed by a transition at v/t = 1 into a staggered ordered phase.

Geometry and sectors
To understand how higher order terms in the expansion affect the ground state, we now study the model in Eq. (37) numerically, using ED. We consider systems with periodic boundary conditions along both directions (i.e. systems on a torus), and focus on two types of clusters which keep all the symmetries of the infinite lattice, following Ref. [21]. Denoting the basis vectors of the triangular lattice by a 1 = (1, 0), a 2 = (1/2, 3/2), these clusters are defined by identifying the sites of the lattice modulo the vectors ( R 1 , R 2 ), where ( R 1 , R 2 ) = m( a 1 , a 2 ) for clusters of type A, and ( R 1 , R 2 ) = (m a 1 + m a 2 , −m a 1 + 2m a 2 ) for clusters of type B. These two types of clusters are shown in Fig. 2. Note that the number of sites in cluster of type A (B) is m 2 (3m 2 ). tors defined by the parities of the number of dimers intersected by closed loops winding around the torus along the two axes. We will denote these sectors by TS(p x , p y ), with p x , p y = 0(1) for even (odd) parity along x and y respectively. As pointed out in Ref. [21], on a cluster with C 6 symmetry, three of these topological sectors are always degenerate since they can be related by C 6 rotations of the lattice. Which three sectors are degenerate depends on the parity of m/2, but in order to understand the spectrum of the problem it is enough to consider the two sectors TS(0, 0) and TS(1, 1), as these two sectors are never related by C 6 rotations.
We consider system sizes of up to 36 sites, i.e. clusters of type B with m = 2 (12 sites) and clusters of type A with m = 4 and m = 6 (16 and 36 sites respectively). We note that to allow for a 12 × 12 order, the number of sites in the system must be a multiple of 12. In addition, for m/2 odd, only the topological sector TS(1, 1) can accommodate this ordering without defects (see Fig. 2).

Exact diagonalization results
We study the successive approximations obtained by working to increasingly higher order in x, denoting by H n the sum of all terms in the effective Hamiltonian up to and including O(x n ). More explicitly, H 0 consists solely of the kinetic term on a plaquette with t = 1, while H 1 contains in addition the potential energy term v as well as the kinetic terms corresponding to hopping on loops of length six. The values of the coefficients in H 1 are given by t = 5/6, and v = t 6,a = 1/6 (note that t 6,b = 0). The Hamiltonian H 2 contains all the terms in Eq. 37 with the corresponding values given in Table 2. We calculate the lowest energy state in each topological sector of H n=0,1,2 , for the physical situation x = 1/6. The values obtained are summarized in Appendix C.3.1. We find that the correction due to second order terms is indeed small compared to the first order ones.
We then interpolate between the Hamiltonians H 0 and H 2 , calculating the low energy spectrum of H(η) = (1 − η)H 0 + ηH 2 . We find that there are no level crossings in the low energy spectrum, and the ground state remains in the topological sector TS(1, 1) for system sizes which can accommodate the 12× 12 order (see Fig. 5 in Appendix C.3.1). The smooth continuity suggests that H 0 and H 2 describe the same phase of matter. Furthermore, we find that, in each topological sector, the wavefunction overlap between the lowest energy state of H 0 and that of H 2 is very close to one, in particular in the topological sector TS(1, 1). More specifically, for the 6 × 6 system, the overlaps are 0.88 and 0.97 for TS(0, 0) and TS (1,1) respectively.
In addition, we compare the dimer-dimer correlations in these states. We find that the correlations in the lowest energy states of H 2 become slightly more uniform compared to those in the lowest energy states of H 0 , but overall display the same features (see Appendix C.3.2). In Fig. 3 we plot the real space dimer-dimer correlations, as well as their Fourier transform, calculated in the lowest energy state of H 2 in the two topological sectors for a 6 × 6 lattice. As can be clearly seen, for the state in TS(1, 1) sharp peaks at k = ±(π/(2 3), π/2) are present, suggesting breaking of translational invariance compatible with the formation of a 12-site unit cell. We note that the six-fold rotational symmetry expected in the ground state is broken in Figs. 3(b,d) by the choice of the set of bonds used in the calculation of the dimerdimer correlation function. The structure factor shown in Figs. 3(b,d) is obtained for the bonds parallel to the lattice basis vector a 2 . When the correlation function is calculated with respect to a set of bonds related by a π/3 rotation on each site, the peaks in the structure factor appear at momenta related by the corresponding π/3 rotation.
Although a better finite size scaling analysis is required to make a conclusive statement regarding the nature of the ground state of the dimer model, we believe that these observations -(i) the similarity of the ground state correlations to those of the "standard" dimer model at v/t = 0 for small system sizes, and (ii) the smooth evolution of the spectrum upon interpolation between the two models -strongly suggest that the ground state remains a 12 × 12 VBS ordered state.

Conclusion
In this work, we considered SU(4) spins in the six-dimensional (self-conjugate) representation, on the triangular lattice, with nearest-neighbor antiferromagnetic interactions. Our DMRG study suggests that the ground state is non-magnetic, but remains inconclusive as to the exact nature of the ground state. We developed and carried out a dimer expansion, which we argued is capable of capturing the low energy properties of the model. The study of the the associated dimer model led us to conjecture that the ground state of the SU(4) model may be a 12-site valence bond solid (VBS).
As the mapping to the dimer model involves an uncontrolled projection, we do not know how to systematically improve it. Hence, a fully conclusive study should return to the original SU(4) spin model. This, however, remains numerically challenging due to the large on-site Hilbert space dimension. As a first step in this direction, we carried out preliminary calculations in addition to those reported in this paper, using the infinite DMRG (iDMRG) method on width-four cylinders. By choosing appropriate boundary conditions, this geometry is com-patible with the 12-site VBS order. However, we did not find signatures of this order in our iDMRG simulations. One possible interpretation is that the non-observation of VBS order is simply due to the effects of finite size or finite bond dimension. Another possibility is that the VBS order is truly absent, indicating some type of spin liquid state without broken symmetries. The proximity of a 2 spin liquid phase in the effective dimer model suggests this as an intriguing possibility. Regardless, this conundrum highlights the challenges of a direct simulation of the original SU(4) problem.
In the study of the effective dimer model, we focused on the parameters corresponding to the SU(4) Heisenberg model, α = −1, β = 1. In the future it would be interesting to explore the full phase diagram of the general dimer model derived, understand if it can realize the 2 spin liquid phase, and identify the nature of the interactions in terms of the SU(4) spins required for this.
In addition, it would be desirable to study in more detail the evolution of the ground state with increasing J H . If the ground state of the spin model at J H = 0 is indeed a 12-site VBS, and if there is, as suggested by our numerics, a direct transition to a three-sublattice ordered state with increasing J H , then this is a Landau-forbidden quantum phase transition. If this is realized via a continuous quantum critical point, then it must be an example of deconfined quantum criticality. It would be interesting to understand the nature of this critical point and test it in numerics.

or in matrix notations
. We next note that Further let v = (x + iy)/ 2 with x, y real six-dimensional vectors with unit norm. Then It is now easy to see that the upper bound on Tr AA T is reached when x ⊥ y.
In the classical limit the Hamiltonian is given by: On the triangular lattice we can rewrite: For antiferromagnetic coupling, J > 0, to minimize the energy, we would like the first term to vanish, and the second to be as negative as possible. Let us denote by x, y, z three real, orthogonal, six-dimensional unit vectors, and define v(θ ) = 1 Then, the matrices A l=1,2,3 = A 2πl 3 satisfy 3 l=1 A l = 0 and Tr A l A T l = 1, thus minimizing the energy on a triangle.
Note that choosing x = e 1 , y = e 2 , z = e 3 , with e n denoting the unit vector along the nth dimension in 6 we obtain a state v(θ ) that belongs to the spin-triplet valley-singlet subspace on a given site (see also Appendix B below and in particular Eq. (58) therein). The classical ground state corresponding to the states v l=1,2,3 = v( 2πl 3 ) on each triangle of the lattice is then, in this case, exactly the 120 • ordered state of the SU(2) spin-ones.

A.2 SU(4) formulation
Here we derive the classical energy function and constraints using the SU(4) formulation, i.e. starting from the Hamiltonian Eq. (5), and using the basis states on the right-hand-sides of the equalities in Eq. (2).
We have in turn T † = T and TrT = 0.
In summary, a classical SU(4) "spin" T satisfies: In this formulation, in the classical limit the Hamiltonian is given by: On the triangular lattice we can rewrite: For antiferromagnetic coupling, J > 0, to minimize the energy, we would like the first term to vanish, and the second to be as negative as possible. Let us denote by z, x two real, orthogonal, three-dimensional unit vectors, and define n(θ ) = cos θ z + sin θ x, Then, the matrices T l=1,2,3 = T 2πl 3 satisfy 3 l=1 T l = 0 and Tr T 2 l = 1, thus minimizing the energy on a triangle. The corresponding state ψ l can be chosen to be where θ l = 2πl/3, so that Indeed, choosing z = (0, 0, 1) and x = (1, 0, 0), we have ψ 0 = 1 2 (iσ y ) ⊗ 1 0 0 0 = |13〉, and ψ l is obtained from ψ 0 through the rotation ψ l = R T l ψ 0 R l , with R l = σ 0 r l , where where y = (0, 1, 0).

A.3 Mapping between the SO(6) and SU(4) formulations
Here we describe the mapping between the SO(6) and SU (4) formulations and show that the classical ground state obtained in the two formulations is indeed the same state. The six basis states |n〉 in Eq. (6) correspond to the following 4 × 4 antisymmetric matrices ψ: Using this mapping one can translate the classical states that optimize the energy on the triangular lattice corresponding to v(θ ) in Eq. (51) to the corresponding ψ(θ ). More explicitly, for x = e 1 , y = e 2 , z = e 3 with e n denoting the unit vector along the nth dimension in 6 we obtain a state where n(θ ) = (sin θ , 0, cos θ ) and τ = (τ x , τ y , τ z ). Thus, ψ(θ l ), with θ l = 2πl/3 + π reproduce the states in Eq. (53) up to an overall phase.

B Large Hund's coupling limit
In the large Hund's coupling limit, i.e. J H /J 1, the term −J H i S 2 i in Eq. (22) requires the total spin at each site to be in the S = 1 representation of SU (2). The associated vector space is spanned by and thus the operatorP i,S=1 = 3 n=1 |n i 〉〈n i | projects the state on site i onto the S = 1 subspace. Note also that |S = 1, s x = 0〉 = |1〉 and |S = 1, s y = 0〉 = |2〉, and thus the S = 1 spin operators can be written as Denoting byP S=1 = iPi,S=1 , where i runs over all lattice sites, to lowest order in J/J H the SO(6) "Heisenberg" Hamiltonian becomes:

C.1 Probing magnetic order
To complement the analysis presented in Sec. 3.3 of the main text, indicating that a finite Hund's coupling, J H , is required to drive the system into a magnetically ordered state, we calculate flavor-flavor correlations in the absence of pinning fields at the ends of the cylinder. More specifically, the correlations calculated are 4 a,b=1 〈T ab 0T ba r 〉, where 0 denotes the origin which we choose to be at the left end of the cylinder, and we consider positions r on the lattice which correspond to the same sub-lattice as the site at 0 when 120 • -order is present. When more than one site on the lattice correspond to the same distance | r|, a symmetrization is performed and an average value for the correlations is used. Resulting correlations are shown in Fig. 4 for cylinders of circumference N y = 3 as J H is increased and the bond dimension is varied. For the maximal bond dimension used of M = 4000, the truncation error was of order 10 where the first (second) term in the product above projects out the anti-symmetric (symmetric) representation. Given a nearest-neighbor covering C = {(i k , j k )} k=1,..,N /2 , as was defined in the main text, we can obtain the corresponding singlet covering MPS by applying the matrix product operator (MPO) representation of the product of projectors N /2 k=1P i k j k to a random initial MPS. Note that to allow for an SU(4) singlet covering state on a system of width N y a bond dimension of 6 N y is required for the MPS. Once the MPS representations of the singlet coverings are obtained, both the overlap matrix, required to solve the generalized eigenvalue problem, and the matrix elements of the projected Hamiltonian can be computed. For the latter, an MPO representation of the original spin Hamiltonian is used. We then solve the generalized eigenvalue problem (since the dimension of the projected Hamiltonian is greatly reduced compared to the one of the original spin Hamiltonian, it can be easily diagonalized using standard sparse diagonalization), both to find the ground state of the projected Hamiltonian in terms of the singlet coverings, and to calculate the gap in the projected problem.
To calculate the overlap of the ground state of the projected Hamiltonian with the ground state of the original spin Hamiltonian, we obtain an MPS representation of the latter using DMRG. For the results presented in Table 1 in the main text, bond dimensions used for the calculation of the ground state were between M = 1000 and M = 2000 depending on system size, resulting in truncation errors ε smaller than 2 · 10 −3 in all cases. A finite truncation error gives rise to an error in the calculation of the overlap that we estimate to be of order ε.

C.3 Exact diagonalization of the dimer model C.3.1 Energies and excitation spectrum of the interpolated Hamiltonian
In Table 3 we summarize the energies of the lowest energy states of the dimer Hamiltonians H n=0,1,2 (where n denotes the order of the expansion in x) obtained using ED. We list the energies of the lowest energy states in the topological sectors TS(0, 0) and TS(1, 1), for three different system sizes with N = 12, 16 and 36 sites. We note that the energies obtained for H 0 reproduce the ones presented in [21] for v/t = 0. Table 3: Energies of the lowest energy states in the two topological sectors, as well as the gap ∆E = E (0,0) − E (1,1) , obtained using ED of the dimer models H n (where n = 0, 1, 2 is the order of the expansion in the parameter x = 1/6) for three different system sizes.

C.3.2 Dimer-dimer correlations
In Fig. 6 we present side by side the real space dimer-dimer correlations for the lowest energy state of H 0 and H 2 respectively, in the two topological sectors TS(0, 0) and TS(1, 1). As was mentioned in the main text the correlations in TS(0, 0) become more uniform for H 2 , while for TS(1, 1) the correlations remain practically unchanged.