Spin-orbital quantum liquid on the honeycomb lattice

In addition to low-energy spin fluctuations, which distinguish them from band insulators, Mott insulators often possess orbital degrees of freedom when crystal-field levels are partially filled. While in most situations spins and orbitals develop long-range order, the possibility for the ground state to be a quantum liquid opens new perspectives. In this paper, we provide clear evidence that the SU(4) symmetric Kugel-Khomskii model on the honeycomb lattice is a quantum spin-orbital liquid. The absence of any form of symmetry breaking - lattice or SU(N) - is supported by a combination of semiclassical and numerical approaches: flavor-wave theory, tensor network algorithm, and exact diagonalizations. In addition, all properties revealed by these methods are very accurately accounted for by a projected variational wave-function based on the \pi-flux state of fermions on the honeycomb lattice at 1/4-filling. In that state, correlations are algebraic because of the presence of a Dirac point at the Fermi level, suggesting that the symmetric Kugel-Khomskii model on the honeycomb lattice is an algebraic quantum spin-orbital liquid. This model provides a good starting point to understand the recently discovered spin-orbital liquid behavior of Ba_3CuSb_2O_9. The present results also suggest to choose optical lattices with honeycomb geometry in the search for quantum liquids in ultra-cold four-color fermionic atoms.


I. INTRODUCTION
The investigation of orbital physics in transition metal oxides has been recently boosted by the possibility to observe orbital excitations with Resonant Inelastic X-ray Scattering [1]. This has been demonstrated in cases where the crystal-field splitting is strong enough to select a unique orbital configuration in the ground state and push the orbital excitations to high energy, hence to separate them from magnetic excitations. However, this is not the only possibility. If the electronic configuration of the transition metal ion is such that several orbital occupations are consistent with the crystal field environment, a situation referred to as orbital degeneracy [2], orbital fluctuations are likely to be much softer. In most cases known until recently, a cooperative Jahn-Teller distortion occurs, resulting in orbital order and gapped orbital excitations, but this needs not be the case a priori, and the search for situations in which orbitals remain fluctuating in the ground state has been very active over the past decade [3][4][5][6][7][8][9]. To which extent this occurs in the triangular system LiNiO 2 [10,11] or in the spinel FeSc 2 S 4 [12] is still debated [13]. Interestingly, a new candidate has been recently put forward, Ba 3 CuSb 2 O 9 [14], a Cu oxide that lives on a decorated honeycomb lattice in which no trace of orbital order could be detected.
On the theory side, transition-metal oxides with orbital degeneracy are generally described by a Kugel-Khomskii model [15] in which spin and orbital degrees of freedom are coupled on each bond. A minimal model to investigate the possibility to stabilize an orbital liquid is the symmetric version of that model defined by the Hamiltonian where S i is a spin-1/2 operator and T i a pseudo-spin 1/2 operator that describes fluctuations of a two-fold degenerate orbital (a and b). Introducing the local basis | =|↑ a , | =|↓ a , | =|↑ b , | =|↓ b , the Hamiltonian exhibits the full SU (4) symmetry and can be written as H = i, j P i, j , where P i, j interchanges the states on sites i and j. The local basis states are often referred to as colors. In fact, the model is the straightforward generalization of the SU(2) symmetric Heisenberg model for S = 1/2 spins which, up to a constant and a factor 2, has the same form when expressed with P i, j operators. The first investigations of this model on various lattices have emphasized the role of 4-site plaquettes, the natural unit to build an SU(4) singlet [5,16,17]. The spontaneous formation of 4-site plaquettes has been proven for an SU(4) ladder [17], and plaquette coverings have been argued to provide the relevant variational subspace for the ground state properties of the SU(4) model on both the square and triangular lattices [5,16], with possibly plaquette long-range order on the square lattice [18,19]. In a variational study based on projected fermionic wave functions a gapless spin-orbital liquid state has been predicted in Ref. [8] on the square lattice. These previous conclusions have recently been challenged for the square lattice [20], for which spontaneous dimerization (as opposed to tetramerization) has been demonstrated on the basis of state-of-the-art iPEPS (infinite projected entangled-pair state) simulations. In the dimerized phase, the dimers are antisymmetric states built out of 2 of the 4 colors, and they form a columnar state. Each dimer is preferentially surrounded by dimers built out of the other 2 colors, leading to long-range color order [20]. In the context of orbital degeneracy, this phase is likely to be ordered. Indeed, one of the states selected by the dimerization consists of alternating pairs of a and b orbitals times a spin singlet. When coupled to the lattice, such a state is expected to undergo a cooperative Jahn-Teller distortion that will stabilize orbitals a and b, hence to lead to orbital long-range order.
In this paper, we consider the symmetric Kugel-Khomskii model on the honeycomb lattice. The first motivation is purely theoretical: since there are no 4-site plaquettes on this lattice, the ground state is unlikely to be a crystal of singlet plaquettes. The second one comes from experiments: the recent observation of a spin-liquid behaviour in Ba 3 CuSb 2 O 9 points to the honeycomb geometry as an outstanding candidate.
The main result of the present investigation is summarized in Fig. 1: The SU(4) symmetric Kugel-Khomskii model is shown to be a quantum spin-orbital liquid with short-range color correlations that follow the pattern illustrated in the top panel, and strong evidence is provided in favor of an algebraic spin-orbital liquid with typical conical singularities in the static structure factor, as shown in the bottom panel.
To reach these conclusions, we have used a variety of analytical and numerical methods: Linear flavor-wave theory (LFWT), infinite projected entangled-pair states (iPEPS), exact diagonalization of finite clusters (ED), and a variational approach based on the Monte Carlo sampling of Gutzwiller projected fermionic wave-functions (VMC). Details about each method can be found in appendix A. These methods are complementary and shed lights on different aspects of the model: LFWT is a good starting point to test for lattice symmetry breaking and color order. iPEPS is a variational approach for infinite systems that has proven to be very successful to check the presence of any kind of long-range order [20,21]. Exact diagonalizations reveal nearly exact information on short length-scale properties and are extremely useful to benchmark other approaches. The variational Monte Carlo of fermionic wave-functions have proven to provide a remarkably accurate description of algebraic quantum liquids. As a first test, we compare in Fig. 2  We note that the LFWT energy is not variational. ready be drawn from this comparison. First of all, among all the Gutzwiller projected fermionic wave-functions we considered, only one is really competitive, the wave-function based on the quarter-filled Fermi sea of standard fermions in the πflux state (see Sec. I C for details). Its energy is much lower than that of the 0-flux state, as well as that of the half-filled Fermi sea of Majorana fermions with 0 or π flux, and these alternative fermionic wave-functions will not be considered any further 1 . Secondly, the agreement between iPEPS, ED and Gutzwiller projected π flux state is quite remarkable. This suggests that all these methods constitute appropriate descriptions within their range of validity. 2

A. Absence of lattice symmetry breaking
Lattice symmetry breaking, be it dimerization or plaquette formation, leads to bonds of different strengths. At the classical level, which consists in minimizing the energy in the subspace of product wave-functions of the form |ψ = i |ψ i , all bonds are fully satisfied. Indeed, since the Hamiltonian of a bond H i j is a simple permutation, ψ i ψ j |H i j |ψ i ψ j = | ψ i |ψ j | 2 is minimal if neighboring states are orthogonal. On bipartite lattices such as the square or honeycomb lattices, it takes only 2 colors to achieve this, and the classical ground state for more than 2 colors is massively degenerate. This degeneracy however can be in principle lifted by zero-point fluctuations. The theory of harmonic fluctuations has been developed previously, and it goes under the name of flavor-wave theory (see appendix A 1 for details). At the harmonic level, the energy of a bond takes the smallest possible value if the two colors of the bond are different from the other nearest neighbors. For the honeycomb lattice, this can be achieved for all bonds simultaneously in an infinite number of ways. The configuration of Fig. 1(a) is the most symmetric one. In this configuration, all bonds have the same surrounding up to color permutations. Other configurations can be generated by exchanging the colors on a stripe of dimers (see Fig. 3(a-b)), leading to a degeneracy of order 2 √ N since, once a direction has been chosen, this can be done independently on all dimer stripes. In all configurations, the energy is the same on all bonds. This is a first hint that, by contrast to the square lattice, the lattice symmetry is not broken for the honeycomb lattice.
The same family of degenerate ground states is obtained with iPEPS if a unit cell with 4 × 4 = 16 different tensors and a small bond dimension D = 2 are used. Upon increasing D, more quantum fluctuations are taken into account, and the symmetric state of Fig. 3(a) is stabilized. In this state, all bonds have the same energy (see Fig. 3(c)). To test how robust this conclusion is, we have challenged it by performing iPEPS simulations using a 2 × 2 unit cell with only four different tensors. This leads to a dimerized state with two types of dimers A and B, which can be distinguished by their dominant colors, and different inter-and intra-dimer bond energies ( Fig. 3(d)). However, unlike on the square lattice, this dimerization vanishes in the infinite D limit, as shown in Fig. 4(a) where the difference in bond energies, , is plotted. Thus, both low-energy states found with iPEPS preserve the lattice symmetry in the large D limit.
We also tested with VMC the stability of the π-flux state towards the dimerization instability shown in Fig. 3(d) by strengthening/weakening the hopping amplitude of the bonds, with the conclusion that the variational energy is minimal in the absence of any dimerization. Similarly, we found no indication toward quadrumerization (SU(4) singlet formation), where the hoppings connected to say red sites (forming a 'tripod') in Fig. 1(a) are modified.
is not variational, but this is an indication that, in the present case, higher orders are likely to be important. Finally, in Fig. 4(b) we show the ED results for the connected bond energy correlations, which is a way to detect dimerization and other bond energy pattern formation tendencies. The correlations decay quite rapidly with distance, making dimerization or other patterns unlikely.
Thus, all methods consistently point towards a state which does not break the lattice symmetry.

B. Absence of SU(4) symmetry breaking
The color-ordered states predicted by LFWT and iPEPS with a small bond dimension (Fig. 3) break the SU(4) symmetry. Here we show that higher-order quantum fluctuations destroy this color order, i.e. that in the ground state the SU(4) symmetry is in fact unbroken.
In Fig. 4(c) we present the iPEPS result for the local ordered moment, where S β α = |α β| are the generators of SU(4) and α, β run over the four different flavors. A finite m implies that the SU(4) symmetry is spontaneously broken in the thermodynamic limit. For both low-energy states found with iPEPS the local ordered moment is strongly suppressed with increasing bond dimension, and most likely vanishes in the large D limit. This shows that quantum fluctuations, which are systematically taken into account by increasing D, eventually destroy the color order so that the SU(4) symmetry is restored. This is in contrast to the same model on the square lattice [20] where the local ordered moment was found to remain finite in the infinite D limit. Consistent results for the flavor correlation function are obtained with ED and VMC for the π−flux state shown in Fig. 4(d), which is decaying rapidly with increasing distance, indicating absence of long range order. The very good qualitative and quantitative agreement between the ED and the VMC results provides substantial evidence that the π−flux state correctly describes the short range physics of the ground state of the Hamiltonian (1). In the next section we will show that the decay predicted by VMC is algebraic, i.e. that the state described by this wave function is an algebraic spin-orbital The orange, outermost hexagon shows the extended Brillouin zone of the triangular lattice (including sites at the centers of the hexagons in the honeycomb lattice), the structure factor is maximal and has a cusp at M = (π, π/ √ 3) and the symmetry related points. K is given by (4π/3, 0), K is (2π/3, 2π/3 liquid.

C. Algebraic spin-orbital liquid
A standard way to describe spin liquids for SU(2) models is based on the fermionic representation of the spin operators [22][23][24][25][26] using a variational wave function [27,28], where the multiply occupied sites are projected out from a suitable chosen, noninteracting Fermi-sea. In the generic case, the Fermisea has a finite Fermi surface. However this is not the only possibility. For the SU(2) Heisenberg model on the square lattice, Marston and Affleck have shown that introducing a πflux per elementary plaquette leads to the formation of Diracnodes [29]. At half filling, the Fermi surface of this π-flux state shrinks to points, and its energy is lower than that of the state with equal hopping amplitudes and a finite Fermi surface. In such a spin-liquid the structure factor is singular at momenta related to the difference between Fermi-points, leading to the algebraic decay of spin correlations. In one dimension, this type of approach leads to an accurate description of the algebraic decay of the correlations for the SU(2) case [30], and as well as for SU(4) using the representation . On the honeycomb lattice, a Dirac-node is already present at the middle of the band without any flux, and the Fermi surface reduces to points at half filling. So the zero flux state would be a good starting point to describe an algebraic spin liquid for the SU(2) Heisenberg model [31][32][33]. However, for the SU(4) Heisenberg model, the band must be quarter-filled, and the equivalent of the Affleck-Marston approach requires to have a Dirac node at the Fermi energy of the quarter-filled system. It turns out that this is achieved in the π-flux state, as shown in Fig. 5(c). As for the square lattice, this state leads to a lower energy than the zero-flux state, as already stated above.
Starting from the noninteracting wave function, with a band populated up to the Dirac-node at ε D = − √ 3t for any of the four flavored fermions, we implemented the Gutzwillerprojection using a variational Monte-Carlo (VMC) sampling. The energy of this wave function, E = −0.894 per site, compares remarkably well with that of iPEPS (see Fig. 2), especially considering that no variational parameter was used. Let us also mention that the state (and the ones related by symmetry) shown in Fig. 1(a) has the maximal weight in the variational wave function.
To investigate the physics of this wave function, we have calculated the spin-spin correlation function as a function of distance. The results clearly demonstrate an algebraic decay P i j − 1/4 ∼ |r i j | −α , with an exponent α between 3 and 4, as shown in Fig. 6. If one considers the honeycomb lattice as built from zigzag chains, these correlations correspond to even distances along one of the zigzag chains, and the exponent should be compared to that of the dominant correlations with wave vector π/2 of a single chain [34]. This exponent is equal to 3/2, a number actually very accurately reproduced by VMC. So color-color correlations decay faster on the honeycomb lattice than on a chain, but still algebraically. This is a rather peculiar situation in view of the standard paradigms: the development of long-range order, as in weakly coupled SU(2) chains in square geometry, or the spontaneous formation of local singlets and exponentially decaying color-color correlations, as e.g. in the SU(4) ladder [17].
This Gutzwiller projected π flux state is actually a prototypical wave function for a phase which should be called algebraic spin-orbital liquid in the present context, in analogy to the algebraic spin liquids which have been discussed in the spin liquid literature [22,24,25,35]. These states are characterized by the algebraic decay of a number of correlations, with wave vectors corresponding to differences between the Dirac cone locii. In our case for example the algebraic spin-orbital correlations are modulated with wave vector M in the standard honeycomb Brillouin zone, and this corresponds to the distance between two Dirac cones in the π flux state. Another important aspect of this wave function is that it has been shown that it can describe an extended phase in parameter space and not just an unstable fine-tuned point [24,36]. Upon adding perturbations of suitable strength, many different phases can be found in the vicinity of an algebraic spin-orbital liquid [24], making the present model an interesting starting point for further explorations of exotic phases in spin-orbital systems.

II. CONCLUSIONS
In conclusion, the results reported in this paper provide very strong evidence that the SU(4) symmetric Kugel-Khomskii model is a quantum spin-orbital liquid, and build a case in favor of an algebraic spin-orbital liquid. Clearly, the present results do not allow to exclude that the quantum spin-orbital liquid is of another type, for instance some kind of resonating valence-bond (RVB) liquid with resonances between 4-site cluster singlets, but our results suggest that the corresponding gap would be quite small. In particular, experience with highly frustrated magnets have shown that, as good as its energy might be, a fermionic variational wave function might fail to capture the correct low-energy physics. This seems for instance to be the case for the SU(2) Heisenberg model on the kagome lattice, for which the variational energy of the algebraic spin liquid wave function [25] is close to the best numerical estimates [37,38], yet DMRG results have given strong evidence in favor of a gapped Z 2 spin liquid [37]. In the present case, the projected fermionic wave function has not just been tested for its energy, but correlations have been shown to be in remarkable agreement with ED up to intermediate distances. So we believe that the case is strong, but not closed.
In any case, the fact that the ground state is a quantum spinorbital liquid is quite firmly established. This result is quite interesting in view of the liquid behavior reported recently in Ba 3 CuSb 2 O 9 . A detailed comparison will clearly require some adaptation of the present model, in particular to take into account the additional magnetic Cu sites present in the system on top of the honeycomb lattice. Still, the absence of orbital order reported in the present paper is consistent with the experimental results.
Finally, we note that the SU(N) Heisenberg model is a rather accurate effective model for the 1/N-filled Mott insulating phase of alkaline-earth metal atoms with N internal degrees of freedom loaded in an optical lattice [39][40][41]. Currently the main issue in that field is to reach low enough entropies to observe correlations typical of long-range order, but the next step will definitely be to realize exotic quantum states. In that respect, the N = 4 case on the honeycomb lattice appears to be a very strong candidate.

III. ACKNOWLEDGMENTS
The ED simulations have been performed on machines of the platform "Scientific computing" at the University of Innsbruck -supported by the BMWF, and the iPEPS simulations on the Brutus cluster at ETH Zurich. We thank the support of the Swiss National Science Foundation, MaNEP, and the Hungarian OTKA Grant No. K73455.
Appendix A: Methods

Linear flavor-wave theory
The linear flavor-wave theory (LFWT) is a method to treat harmonic quantum fluctuations on top of a mean-field (or Hartree) solution based on a site-factorized variational wave function [42,43]. It starts from a Schwinger boson representation of the SU(4) operators with 4 types of bosons. In a mean-field ground state, each site has a well defined color. The corresponding boson is assumed to condense, and the resulting Hamiltonian is a bosonic quadratic form. On a nearestneighbor bond with color α on site i and color β on site j, it is given by H fw = A † i j A i j −1, with A † i j = b † α, j +b β,i where b † and b are bosonic operators. In a given mean-field ground state, the LFWT Hamiltonian is the sum of independent Hamiltonians that describe the motion of bosons on the connected clusters spanned by pairs of colors. The zero-point energy per bond tends to increase with the cluster size, and is minimal on a 2-site cluster, when the ground state energy of the Hamiltonian is equal to −1, i.e., there is no zero-point contribution to the energy. By contrast, larger clusters lead to finite frequencies, hence to strictly positive contributions to the zero-point energy.

Infinite projected entangled-pair states (iPEPS)
A projected entangled-pair state (PEPS) [44,45], also called tensor product state, is a variational ansatz where the wave function of a two-dimensional system is efficiently represented by a product of tensors, with one tensor per lattice site. It can be seen as a two-dimensional generalization of matrix product states -the class of variational states underlying the famous density matrix renormalization group (DMRG) method [46]. On the square lattice each tensor T p i jkl has a physical index p carrying the local Hilbert space of a lattice site with dimension d, and four auxiliary bond indices i, j, k, l with dimension D which connect to the four neighboring tensors. Thus, each tensor consists of dD 4 variational parameters, and by changing the bond dimension D the accuracy of the ansatz can be systematically controlled. A D = 1 PEPS simply corresponds to a site factorized wave function (a product state), and upon increasing D quantum fluctuations can be systematically taken into account.
An infinite PEPS (iPEPS) [47,48] consists of a unit cell of L x × L y = N T tensors which is periodically repeated in the lattice to represent a wave function in the thermodynamic limit. We use the iPEPS method developed for the square lattice described in Refs. [20,47,49] to simulate the model on the honeycomb lattice by mapping it onto a brick-wall lattice as illustrated in Fig. 7(a). The bond dimension of the auxiliary bonds connecting two sites which do not directly interact (dotted lines) can be chosen as D = 1. This ansatz is equivalent to an iPEPS with only three auxiliary indices on the honeycomb lattice. The advantage of this mapping is that the codes developed for the square lattice only require minor modifications to simulate models on the honeycomb lattice.
Describing the iPEPS method in full detail is beyond the scope of this paper and we therefore only mention the most important technicalities with corresponding references and details on the simulation parameters in the following. For an introduction to PEPS and iPEPS we refer to Refs. [45,49].
The optimization of the tensors (i.e. finding the best variational parameters) is done through imaginary time evolution, first with the so-called simple update (see Refs. [49,50]) which is equivalent to the time-evolving block decimation method in one dimension [51,52]. The solution is used as an initial state for an imaginary time evolution using the full update [47,49], which is computationally more expensive, but more accurate than the simple update, since the full wave function is taken into account for the truncation of a bond index. We use a second order Trotter-Suzuki decomposition with time steps down to dτ = 10 −3 . For large values of D a larger time step dτ = 10 −2 is used, where the estimated Trotter error is small compared to symbol sizes (e.g. below 0.5% for the ordered moment m (2)).
Expectation values of observables can be computed by contracting the tensor network, i.e. by computing the trace of the product of all tensors. For the approximate contraction of the iPEPS we use the corner transfer matrix method described in Refs. [49,53,54]. The accuracy of this contraction can be controlled by the so-called boundary dimension χ, where we used values up to χ = 500 (typically up to χ = 350) for large D. Observables like the energy and the ordered moment are extrapolated in χ, with an extrapolation error being small compared to symbol sizes. An example is shown in Fig. 7(b) for the ordered moment m (2).
To increase the efficiency of the method we implemented a Z q symmetry in the tensors (see Refs. [55,56]) which is a discrete subgroup of SU (4). This implies that the tensors have a block structure (similar to a block diagonal matrix) which reduces the computational cost considerably.
Since iPEPS is an ansatz for an infinite system, symmetries such as SU(4) or translational symmetry may be spontaneously broken. In order to test different types of translational symmetry breaking we compared the variational energies obtained with different unit cell sizes. We found two competing low-energy states with unit cell sizes 2 × 2 and 4 × 4, shown in Fig. 3, which have similar energies for large bond dimension. We note that broken symmetries can be an artifact of a finite bond dimension D, and that the symmetry can be restored if D is sufficiently large. In other words, a classical or a low entanglement solution (small D) might exhibit order, but this order can be destroyed by quantum fluctuations which are systematically included upon increasing D. Therefore, it is important to study order parameters as a function of D to verify if they are finite in the large D limit.

Exact diagonalization
We have performed exact diagonalizations of the Hamiltonian (1) for selected finite size samples of up to N = 24 sites. The energies of the samples with 8,14,18 and 24 sites are reported with star symbols in Fig. 2. Note that only the samples with 8 and 24 sites can form a SU(4) singlet ground state, and this explains the relatively high energy per site for the samples with 14 and 18 sites. Note also that due to computational limitations we were only able to calculate the energy and eigenfunction of the ground state in the completely symmetric representation of the spatial symmetry group of the N = 24 sample (the Hilbert space in this symmetry sector contains 4'008'417'658 states, including a cyclic color permutation symmetry). Since this sector is the absolute ground state for the N = 8 sample, we expect this sector to host the ground state for N = 24 as well.

Fermionic Variational Monte Carlo
The variational wave function for the algebraic spin-orbital liquid is a Gutzwiller projected noninteracting Fermi sea at quarter filling: where j α l denotes the position of the l-th fermion with color α, and the summation is over the all N!/((N/4)!) 4 possible distributions of the fermions so that each site is single occupied (i.e. { j α } ∩ { j β } is an empty set for α β). The weight of each configuration is given by the Slater-determinant, where ξ k ( j) is the amplitude of the fermion at site j in the k-th eigenfunction of the hopping Hamiltonian The expectation values of operators with this variational wave function are evaluated by a Monte Carlo sampling [27,28]. The variational wave function with Majorana fermions is described in detail in Ref. [8].
The different trial states correspond to different choices of the χ i, j hopping parameters. In the π-flux state |χ i, j | = 1 and the phases are chosen so that an electron hopping around the hexagon picks up a minus sign, 6 k=1 χ i k ,i k+1 = −1, where the product is over the bonds of a hexagon. This can be realized in many ways. Here we choose real hopping amplitudes, where every hexagon has a single bond with a χ i j = −1, while the rest of the bonds have +1, as shown in Fig. 5(a). Furthermore, we allow for antiperiodic boundary conditions when a degeneracy of quarter filled Fermi sea needs to be removed for a given cluster.
We use Monte Carlo sampling of the projected wave function |Ψ to evaluate physical quantities. The elementary step is the exchange of two randomly chosen fermions with different colors. To speed up the convergence we used importance sampling, with acceptance ratios defined according to the Metropolis algorithm: the new configuration is always accepted if its weight w new is higher than the weight w old of the old configuration, but for w new < w old it is accepted only with a probability w new /w old . The configurations are thus represented with a probability p { j} proportional to their weight in the wave function: where . . j α N/4 in the projected fermionic state |Ψ , see Eq. (A2). We set the number of elementary steps between two measurements large enough to avoid autocorrelation effects. See Table I for details. We also performed a binning analysis as a further test for the independence of the measurements. The statistical error for different bin sizes did not show any change, verifying once more that the sampling distances were large enough. Furthermore, for each system we made several (5 -10) runs with randomly chosen starting configurations to independently verify the error bars obtained from the binning analysis.  We measured diagonal and off-diagonal operators. The spin-spin correlation function, the average of the off-diagonal P k,l can be expressed using the diagonal n β k n β l operator (where n β k is the occupation number on site k for the fermion of color β) as P k,l = 20 n β k n β l − 1, supposing that the ground state is a singlet wave function -as it is the case when the hopping Hamiltonian is independent of the colors. The measurement of the diagonal n β k n β l correlation functions is quite simple using the importance sampling, where N({k, l} ⊂ { j β }) denotes the number of times both k and l sites were occupied with β fermion among the N MC measured configurations.
With a little more effort one can directly calculate the off-diagonal P k,l correlation functions as well. Using the fermionic representation for the P k,l exchange operator, the convenient form that follows the convention of the fermion ordering in the wave function is given as In this case one follows the same importance sampling as before, although the measurement itself is more complicated: where { j } is the configuration that leads to { j} by exchanging the color of fermions on sites k and l, and the sum in the last equation is over the measured configuration { j} MC . Following Eq. (A7), the sign s k,l ({ j}) is 1 if the colors of fermions on sites k and l in the configuration { j} are the same, and −1 if the colors are different. We have explicitly verified that the Eq. (A5) holds. Similarly one can calculate P i j P kl as well, here for each { j} configuration |{ j} = P i j P kl |{ j } . For the sign one should take s k,l ({ j})s i, j ({ j}), here we assumed that the sites i, j,k, and l are all different.