Competition between threesublattice order and superfluidity in the quantum 3state Potts model of ultracold bosons and fermions on a square optical lattice
Abstract
We study a quantum version of the threestate Potts model that includes as special cases the effective models of bosons and fermions on the square lattice in the Mott insulating limit. It can be viewed as a model of quantum permutations with amplitudes and for identical and different colors, respectively. For it is equivalent to the SU(3) Heisenberg model, which describes the Mott insulating phase of 3color fermions, while the parameter range can be realized in the Mott insulating phase of 3color bosonic atoms. Using linear flavor wave theory, infinite projected entangledpair states (iPEPS), and continuoustime quantum MonteCarlo simulations, we construct the full phase diagram, and we explore the properties for . For dominant antiferromagnetic interactions, a threesublattice longrange ordered stripe state is selected out of the ground state manifold of the antiferromagnetic Potts model by quantum fluctuations. Upon increasing , this state is replaced by a uniform superfluid for , and by an exotic threesublattice superfluid followed by a twosublattice superfluid for . The transition out of the uniform superfluid (that can be realized with bosons) is shown to be a peculiar type of KosterlitzThouless transition with three types of elementary vortices.
pacs:
67.85.d, 75.10.Jm, 02.70.cI Introduction
The possibility to realize Ncolor Fermi or BoseHubbard models with ultracold atoms in optical lattices has attracted increasing interest recently. Wu et al. (2003); Honerkamp and Hofstetter (2004); Cazalilla et al. (2009); Fukuhara et al. (2009); Gorshkov et al. (2010); Taie et al. (2010); Sugawa et al. (2011) With fermionic alkalineearth atoms (and with the alkalineearth like atomYb), it is possible to realize systems with up to N=10 different flavors (or colors) of fermions by exploiting different nuclear spin states, Gorshkov et al. (2010) while 3color bosonic systems can be created with bosonic alkali atoms with a hyperfine spin , e.g. Na, K, and Rb. A good starting point to describe these systems is the SU(N) Hubbard model given by
(1)  
where denotes a sum over all pairs of nearest neighbors, and the operators and create and annihilate a bosonic or fermionic particle with color on site , respectively. The second term describes an onsite interaction with amplitude between particles with different colors , with , whereas the third term corresponds to the onsite interaction between particles with the same color. For fermions, the last term vanishes because of the Pauli exclusion principle. For femionic alkalineearth atoms, this description is very accurate because, since the interactions between the particles do not depend on the nuclear spin, these systems exhibit a full SU(N) symmetry. In the bosonic case, there is in general an additional onsite interaction proportional to the square of the total spin that can be antiferromagnetic (as in Na) or ferromagnetic (as in Rb). Law et al. (1998) This interaction is small  it vanishes if the scattering lengths of the spin0 and spin2 channels are equal  and we will not consider it in the present paper.
At filling (one particle per site) and for large enough repulsion the system is in a Mott insulating phase, and color fluctuations can be described by the effective Hamiltonian:
(2) 
where is defined by . The values of the coupling constants and depend on the statistics of the particles. To second order in and , they are given by for fermions, and by and for bosons. Kuklov and Svistunov (2003); Altman et al. (2003) Thus, the experimentally accessible parameter range is given by for fermions, and for bosons. Throughout the paper we will use the notation and to parametrize the Hamiltonian by a single parameter .
In the fermionic case with equal couplings () the Hamiltonian corresponds to the antiferromagnetic (AF) SU() symmetric Heisenberg model. For , we recover the wellknown AF Heisenberg model. The case of larger has been the subject of many theoretical studies, in which a large variety of different ground states have been found. For example, for , a three sublattice Néel ordered state has been predicted on the square and triangular lattice,Tóth et al. (2010); Bauer et al. (2012) whereas on the honeycomb and the kagome lattice the ground state can be understood as a generalized valencebond solid state. Arovas (2008); Hermele and Gurarie (2011); Corboz et al. (2012a); Zhao et al. (2012); Corboz et al. (2013) Even more exotic ground states have been predicted for larger , such as a dimerNéel ordered state on the square lattice, Corboz et al. (2011a) an algebraic spinorbital liquid on the honeycomb lattice Corboz et al. (2012b) (for ), and chiral spin liquids. Hermele et al. (2009); Hermele and Gurarie (2011); Szirmai et al. (2011)
The anisotropic case of the model of Eq. (2) () can be seen as a generalization of the anisotropic Heisenberg (XXZ) model for , for which the phase diagram on the square lattice is well known: if the diagonal terms in the Hamiltonian are dominant and the ground state is an antiferromagnetic (ferromagnetic) state for (), whereas if the offdiagonal terms are dominant, and the ground state is a superfluid with uniform (2sublattice) superfluid order parameter for (). However, for larger the ground state phase diagram is unknown.
In this article we focus on the case , and study the full ground state phase diagram of the model (2) on a square lattice, using linear flavor wave theory (LFWT, Sec. II), infinite projected entangledpair states (iPEPS, Sec. III) and continuoustime worldline quantum MonteCarlo (CTQMC, Sec. IV). With the latter method we also explore the finite temperature phase diagram for , in which case there is no negative sign problem.
Ii Linear flavor wave theory (LFWT)
ii.1 The method
The LFWT is a generalization of the spinwave expansion of models to .Papanicolaou (1984) The large parameter of LFWT is the filling denoted by (it is usually called spin and denoted by for ). The representation used on each site corresponds to a Young tableau with boxes arranged horizontally.
First, one chooses a classical state to which quantum fluctuations will be added. These states (the product states, described in Sec. II.2) are the ground states of the model and are the subject of the next subsection. We then perform a local SU() transformation such that the image state has bosons of color on each site. Next the transformed Hamiltonian (expressed in terms of new operators) is expanded in powers of and its truncation to quadratic operators is solved by a Bogoliubov transformation. The results for our model are given in Sec. II.3. The possible phase transitions are discussed in Sec. II.4.
ii.2 Product states
Let us define as an operator creating a bosonic particle of color at site and as the vacuum of particles. A state at site can be fully described by a unitary component complex vector such that . The local U(1) phase of the vector has no importance. Thus the state corresponds to a single point of the complex projective space CP(2), but the vectorial notation will be conserved to keep notations simple.
Product states are simply tensor products of single site states, i.e. nonentangled states. They are the generalization of the classical spin states to the color models (for , product states are fully magnetized states characterized by a unitary threedimensional real vector at each site). For a lattice of sites, a product state is described by complex numbers, to be compared with for any state of the full Hilbert space. In this reduced space, the minimal energy is variational. Since the square lattice with nearestneighbor interactions is bipartite, the lowest energy state minimizes the energy of each link which depends only on and :
(3) 
Depending on , there are four phases of product ground states for (see the inner circle of Fig. 1):

the 2sublattice superfluid (2subSF) phase (). There are three distinct ground state manifolds. The states of one of them are characterized by a phase such that the vectors on the two sublattices are and . The two remaining manifolds are obtained by a cyclic permutation of the vector components. This leads to a set of ground states homeomorph to . The ground state energy per bond is .

the antiferromagnetic (AF) phase (): the vectors are chosen among , and and set on the lattice in such a way that two neighboring sites have different vectors. These states are the ground states of the classical antiferromagnetic Potts model, which has an infinite degeneracy, and a ground state energy per bond of .

the uniform superfluid (USF) phase (): all vectors are identical on the lattice and the components fulfill . The ground state manifold is homeomorph to a torus (). The states can be parametrized by with , with an energy per bond.

the ferromagnetic (F) phase (): all vectors are identical on the lattice and equal to either , or . These states are the ground states of the classical ferromagnetic Potts model. In fact, these states are the exact quantum ground states in this range of . The ground state manifold is homeomorph to , and the energy is per bond.
Two values of do not break the SU(3) symmetry (red points in Fig. 1):

: This is the AF SU(3) symmetric point. Every state with perpendicular neighboring vectors is a ground state and there exist an extensive number of such states. At this point, the 2 sublattice AF and the 2subSF phases of Fig. 1 are linked by a global SU(3) transformation.

: This is the F SU(3) symmetric point. The ground states are the ferromagnetic states, i.e. those with uniform . The ground state manifold is homeomorph to . At this point, the F and the USF phases of Fig. 1 belong to this manifold, and they are linked by a global SU(3) transformation.
For , since the square lattice is bipartite, we have the same properties for as for and this change of parameters can be compensated by a multiplication of the vectors of one sublattice by the matrix (rotation of the spins of around the axis in the corresponding spin model). It is this property that solves the sign problem of the Heisenberg spin model on the square lattice. Then and are too SU(2) symmetric points for . For we have a reduced SU(2) symmetry at these angles: states with one vector component equal to zero everywhere and their transforms have the same energy for and , thus the subset of such ground states of a SU(3) symmetric point are sent to degenerate subsets for opposite . It gives a continuum of ground states for including both the 2subSF and F states. At , we do not have such a continuum since the vectors of the ground states have three non zero components.
ii.3 Quantum fluctuations at the harmonic level
Harmonic quantum fluctuations lift some of the degeneracies found in the product state analysis of the previous section, however, not all of them: when degenerated ground states are connected via a symmetry of the Hamiltonian, then the harmonic (and higher order) fluctuations cannot lift their degeneracy. This is the case at the SU(3) symmetric point of the FUSF boundary (), because the two product states can be transformed into each other by a global SU(3) transformation.
The degeneracy is not lifted either at the F2subSF phase boundary (), but for a different reason: the F and 2subSF states both belong to a continuum of states linked by a set of transformations homeomorph to SU(2), but that are not in the group of the Hamiltonian symmetries (the origin of this reduced symmetry was given in the previous section). The states with vectors and on both sublattices have the same energy for any complex numbers and such that . The ground state manifold is thus extended to SU(2) and contains both the F and the 2subSF state. Since such states are exact quantum ground states (they have the same energy per bond as an isolated bond), the degeneracy remains at any order.
At all other boundaries of the LFWT phase diagram and for any value of in the AF phase (since there is an extensive degeneracy of productstates), the term of order will lift the classical degeneracy. To calculate the first order quantum correction, we first apply a local SU() transformation to the operators that become such that the vector is transformed into . We then perform the following replacement by bosonic operators in the modified Hamiltonian:
(4)  
( and are colors different from ) and solve the Hamiltonian on the square lattice neglecting the terms.
Among the extensive number of AF states, we will only consider the 2 sublattice order (noted 2subAF) and the 3 sublattice stripe orderBauer et al. (2012) (noted 3subAF) depicted in Fig. 2. Each state of Sec. II.2 is stable up to the first order correction in the interval of where is is the ground state. Moreover, the subAF are both stable until . For a product state, stability means that the derived quadratic Hamiltonian has a ground state (the Bogoliubov transformation is possible). At the classical boundary between two orders, we can compare the corrections of order of both phases and determine which one has the lowest energy to this order. Since simply taking does not give an accurate estimate of the energy (one would need further orders), one cannot be quantitative regarding the boundary, and we just indicate in which direction it shifts due to harmonic fluctuations (arrows on Fig. 1).
The corrections for the 2subAF and 2subSF phases are the same at the AF SU(3) symmetric point because of the symmetries, but since 2subAF is degenerate with other product states with different symmetry properties, the first order corrections of these other states will be different and possibly more favorable. Thus, the boundary between AF and 2subSF phases can still be pushed towards larger ’s.
We now give the numerical results and the evolutions for all boundaries. At the AFUSF boundary, the correction per site is for the 3subAF, for the 2subAF and for the USF phase. Thus we expect that the USF phase is stable until .
At the AF2subSF boundary, the correction is for the 3subAF and for the 2subAF and 2subSF. It suggests that the 3subAF survives for . But the degenerate 3subAF states obtained by global SU(3) transformations at the SU(3) symmetric point do not have the same classical energy as soon as we are not exactly at the point. For smaller angles, the better states are those described previously, with vectors among , and . For larger angles, we can reach on a link with vectors among , and , with and . It will thus be this 3 sublattice superfluid stripe order dressed by quantum fluctuations that will survive for , noted 3subSF in Fig. 1 (the ground state manifold is homeomorph to ).
All these predictions will be confirmed by the iPEPS and/or QMC measurements.
ii.4 Symmetries and phase transitions
From the broken symmetries of the productstates of Sec. II.2, we can have an idea of the finite temperature phase transitions. Equivalently, we can look at the topology of the ground state manifold. The MerminWagner theoremMermin and Wagner (1966) states that a connected ground state manifold cannot give rise to a finite order phase transition in two dimensions, except at zero temperature. This is the case for the USF phase.
When the ground state manifold consists of several disconnected regions (when some points cannot be linked by a continuous path in the manifold), then there is at least one first or second order phase transition. The F and 2subSF phases will give rise to such a phase transition with a order parameter. The 3subAF and 3subSF phases break a symmetry, that can be restored through one or more transitions ( is the group of permutations for three elements). Moreover, it is not excluded that 2subAF becomes more stable than 3subAF at finite temperature,Tóth et al. (2010) and that it gives rise to its own transitions.
Besides the finite order phase transitions, infinite order phase transitions associated to topological defectsMermin (1979) can occur when the components of the ground state manifold are not simply connected (some paths are not continuously shrinkable to a point). This is the case for the USF, 2subSF, and 3subSF phases, possessing pairs of defects at low temperature. The defects encountered in the 2subSF phase are in (the first homotopy group of is ), whereas those of the USF and 3subSF phases are in . Thus, a KosterlitzThouless (KT) transitionKosterlitz and Thouless (1973) could occur in the 2subSF phase and a KT type phase transition in the USF and 3subSF phases. The defects encountered in these two phases are the same as the dislocations of triangular two dimensional crystals and the phase transition is similar to its melting. Nelson (1978)
The case of the USF phase is detailed below. Among the three phases with defects, USF is the only one in which defects do not coexist with a discrete symmetry breaking and for which the occurence of a topological phase transition makes no doubt. In the two other phases, 2subSF and 3subSF, we have both defects and discrete symmetry breaking. It can lead to a variety of phase transitions, with the possibility of fractional vortices, as in the case of the AF XY model on the triangular lattice. Korshunov and Uimin (1986); Korshunov (2002); Hasenbusch et al. (2005) The larger discrete symmetry and the more complicated topological defects for the 3subSF phase require further investigations about the nature of the phase transitions, which go beyond the scope of the present paper.
We concentrate here on the USF topological phase transition. If such a transition occurs in the 3subSF phase, the same description can be easily adapted. We neglect the component fluctuations and construct a subset of product states interacting via the same Hamiltonian and possessing the same superfluid properties. More precisely, we only consider product states whith . We can fix the gauge such that and parametrize a the state of a site state with only two real parameters and . Since , we note so that the sign of the different quantities that will follow is obvious. The energy on a bond then reads:
(5)  
where and . An expansion in and gives:
(6) 
with the energy of a uniform state . Note that a uniform state is labelled by a point of the order parameter space, whose topological properties will be discussed below. Taking the continuum limit of this model on the square lattice, we obtain
(7) 
If we look for the local minima of the energy, we obtain the two conditions
(8) 
These conditions are verified if and are constant (global minimum of energy, per surface unit), but also for configurations with local singularities called vortices (metastable states). This model looks very similar to two uncoupled XY models , that also leads to Eq. (7) in the continuum limit, but they are different because of the nature of the topological defects. The topological defects are linked to the set of equivalent points in the covering group of the order parameter space: the plane .
In the model, the covering group of the order parameter space is the real line , and two points are equivalent if their distance is a multiple of . A closed path in the order parameter space corresponds to a path in between two equivalent points, characterized by a relative integer , the distance between these two points divided by , also called the winding number.
In the two uncoupled models, the covering group of the order parameter space is the real plane , and equivalent points are on a square lattice of spacing . A closed path in the order parameter space is now characterized by two relative integers .
In our model (two coupled models), the covering group of the order parameter is still a plane, but equivalent points form a triangular lattice of spacing (related to the Burger’s vectors of a triangular crystalNelson (1978)). The most symmetric way to characterize the topological properties of a closed path in the order parameter space is to use three relative integers defined by
(9)  
(10)  
(11) 
linked by and represented on Fig. 3.
We now want to calculate the energy of a configuration with a vortex of winding numbers at position and satisfying Eq. (8). Since the problem has a spherical symmetry, and only depend on the modulus . Thus and . We obtain from Eq. (7)
(12) 
We needed two cutoff length: the lattice spacing and the lattice size . Six elementary vortices exist, with an energy : , , and their antivortices with opposite ’s.
We have here detailed the case of a three color USF phase, but similar color USF phases exist. Then, the number of elementary vortices would be . The transition associated to these vortices is numerically characterized through the winding number (see Sec. IV).
Iii Infinite projected entangledpair states
iii.1 Method
An infinite projectedentangled pair state (iPEPS) is a variational tensor network ansatz for twodimensional ground state wave functions in the thermodynamic limit,Verstraete and Cirac (2004); Jordan et al. (2008) where the accuracy can be systematically controlled by the socalled bond dimension . It consists of a unit cell of rank5 tensors with one tensor per lattice site, which is periodically repeated on the lattice. Corboz et al. (2011b) The first index of a tensor carries the local Hilbert space of the corresponding lattice site, and the four remaining indices, called auxiliary bonds, connect the tensor to its four nearest neighbors. Details on the iPEPS method can be found in Refs. Jordan et al., 2008; Corboz et al., 2010, 2011b, and also in Ref. Bauer et al., 2012, where we used the same approach to study the AF SU(3) symmetric point of the current model.
For the experts we note that the optimization of the tensors is done through imaginary time evolution using the socalled simple update,Vidal (2003); Orús and Vidal (2008); Jiang et al. (2008); Corboz et al. (2010) and we verified some simulations results also with the full update.Corboz et al. (2010) The cornertransfer matrix method Nishino and Okunishi (1996); Orús and Vidal (2009); Corboz et al. (2010) is used to contract the tensor network, where the error of the approximate contraction can be controlled by the boundary dimension . The error due to the finite is small compared to the symbol sizes. To simulate the state with 3sublattice (diagonal) order we used tensors with symmetry Singh et al. (2011); Bauer et al. (2011) to improve the efficiency. In the superfluid phases we cannot exploit this symmetry since it is spontaneously broken in the thermodynamic limit.
By performing simulations with different (rectangular) unit cell sizes in iPEPS we can represent states with different type of lattice symmetry breaking. For example, a twosublattice state requires a unit cell with two different tensors and arranged in a checkerboard order, whereas a threesublattice state requires a unit cell with three tensors arranged in a stripe order. In practice we perform simulations with different unit cell sizes and check which ansatz yields the lowest variational energy. Finally we note that the optimization scheme (imaginary time evolution) requires an ansatz with at least two different tensors and , and this is why the smallest unit cell we consider is (containing two different tensors), even for the uniform phases.
iii.2 Results
Our iPEPS results are summarized in Fig. 4: In Fig. 4(a) we present the local color densities of the three colors, given by piecharts, on each site in the iPEPS unit cell for the different states in the phase diagram. Figure 4(b) shows the ground state energy obtained for different values of in the whole parameter range of , and in Fig. 4(c) we plot the following order parameters,
(13)  
(14) 
averaged over all sites in the unit cell. A finite value of indicates that there is longrange color (diagonal) order, whereas a finite implies offdiagonal longrange order (a superfluid phase).
In the following we discuss the results for the individual phases and then explain how we accurately determined the phase boundaries.
In agreement with the findings in the previous sections, iPEPS predicts a ferromagnetic (F) phase between . In this phase each site is occupied by the same color, as shown in Fig.4(a). Since the ferromagnetic state is a product state it can be trivially represented by an iPEPS, i.e. with a bond dimension . This is why the energies obtained with and the shown in Fig.4(b) are identical in this phase.
In the range iPEPS confirms the uniform superfluid phase (USF). The states in this phase do not have color order, as can be seen in Fig.4(a) where all color densities on all sites are the same. However, we find a finite, uniform offdiagonal order parameter , indicating a superfluid phase.
In the range we clearly find the lowest variational energy with a 3sublattice iPEPS ansatz. This range also includes the SU(3) symmetric point , where iPEPS previously found a 3sublattice Néel ordered state.Bauer et al. (2012) For we find a 3sublattice diagonal order (3subAF in Fig. 4(a) with ) where each site is occupied by one dominant color, and the three colors are arranged in a stripe pattern. The three sublattice stripe structure persists for , however, in this range it is the offdiagonal terms in which are finite and the diagonal terms vanish (3subSF phase).
We also confirm the 2subSF phase with iPEPS in the range . Fig. 4(a) shows that the state has a uniform diagonal order (each site is equally populated by only two colors with the third color being absent), however, the offdiagonal terms exhibit a twosublattice structure (not shown).
Next we discuss the phase transition between the USF and 3subAF phase. We find clear evidence of a first order phase transition, with a jump of the order parameter at the transition value. In order to accurately determine the critical value we make use of the hysteresis effect in the vicinity of a first order phase transition: When we initialize a state in the USF phase and change to a value slightly above then the state will remain in the USF phase (because it is metastable). Similarly we can obtain the energy of the 3subAF state slightly below . We then determine from the intersection of the energies of the two states, as shown in Fig.4(d). The error bar of the phase boundary includes the uncertainties when extrapolating the energies to the infinite limit.
In a similar way we have determined the location of the first order transition between the 3subSF and the 2subSF phase, as shown in Fig. 4(e).
We did not find that the phase boundaries of the ferromagnetic phase change when going to larger bond dimensions. Thus, the location of the phase transitions is the same as predicted classically, i.e. and .
Iv Continuous time world line MonteCarlo with cluster updates
The iPEPS numerical simulations give access to informations at zero temperature. Using quantum Monte Carlo (QMC) methods, finite temperature observables can be measured. Moreover, QMC methods work in the full Hilbert space and thus is a non biased numerical method. They give interesting complementary informations to the iPEPS and LFWT results.
Unfortunately, the socalled sign problem occurs as soon as some nondiagonal elements of the Hamiltonian are positive. To study the model of Eq. (2), we chose a basis of states consisting of the product states , where is the color of the boson on site and is the number of sites of the lattice. These are eigenstates of all the . The sign of the nondiagonal elements of the Hamiltonian matrix is the sign of . The sign problem only allows us to simulate this model for .
In some specific cases, the sign problem can be overcome. On an infinite or open one dimensional chain, the physical properties of the color model are unchanged by the transformation , whatever . This allows us to use the QMC method whatever the value of , although mainly the SU() symmetric point was considered.Frischmuth et al. (1999); Messio and Mila (2012) In two dimensions, the same transformation still conserves the physical properties for (see discussion in Sec. II.3), but not any more for , where qualitatively different phases can be obtained for two opposite values of , as can be seen in the phase diagram of Fig. 1. QMC is thus limited to .
iv.1 The algorithm
We use a continuous time worldline algorithm with cluster updates,Prokof’ev et al. (1996); Beard and Wiese (1996) adapted to our color model. We express the partition function as a path integral over the configurations , where is the imaginary time going from 0 to ( is the temperature) and is a basis state. The functions that contribute to can be represented by and by a set of world line crossings at which the colors of two neighboring sites and at a time are exchanged. A local configuration on a link at time is represented by
(15) 
Cluster algorithms are well known for 2color (spin) models. Here we generalize them by chosing randomly two different colors and out of the and by constructing clusters on which only these two colors are encountered. The steps to construct a cluster are the following. We first randomly place the graphs depicted in the first column of Tab. 1 on each link of the configuration using a Poisson distribution whose time constant is given in the last columns of Tab. 1 and accepting them if (if a color being neither nor appears in the local configuration, the graph is rejected). Then we assign graphs to the worldline crossings between and colors with a probability proportional to . At the places where no graph has been attributed, we follow the path with the same color. Finally we follow each constructed cluster and exchange and on it with a probability . This completes one Monte Carlo step. During the whole simulation, steps are performed.
We use periodic square lattices of linear size . The goal is now to calculate averages of observables such as the correlations over the sampled configurations.
We define the diagonal correlation on site as
(16) 
and the associated structure factor
(17) 
with the position vector of site . This structure factor is normalized in such a way that .
iv.2 Results
We have studied the sector of Fig. 1 using the previously detailed algorithm. Confirming LFWT and iPEPS results, we found two different phases at low temperature: the USF and the 3subAF stripe order. We studied the thermal phase transitions generated by these two phases (see Fig. 5). For USF, it is the KT type phase transition linked to the winding number described in Sec. II.4. The 3subAF stripe order gives rise to a first order transition due to the symmetry breaking.
The clearest sign of the symmetry breaking is the structure factor in the reciprocal space. Bragg peaks appear at two points of coordinates or depending on the direction of the stripes (see Fig. c). Their evolution with the size of the lattice as indicates long range order. The transition temperature increases with more negative ’s (see Fig.5) and corresponds to a first order transition: coexistence of ordered and unordered phases occurs during the simulations. When , no trace of the 3subAF phase is found at any temperature. At large temperature, the structure factors have peaked maxima at point (similar to Fig. b), but the finite size effects indicate no long range order associated to these wave vectors.
For the KT type phase transition of the USF phase, we looked at the evolution of the stiffness with the temperature.Pollock and Ceperley (1987) The critical temperature is deduced from the intersection of the curves with the line ,Nelson and Kosterlitz (1977); Harada and Kawashima (1997) extrapolated to the thermodynamic limit. At it is (see Fig. 7). Even if this phase cannot be characterized by the structure factor as there is no long range order, it looks quite different from the high temperature phase. The peak at is much broader, as illustrated on Fig. a.
V Results and conclusions
We have studied a quantum version of the state Potts model on the square lattice which can be seen as a generalization of the well known XXZ model. We have presented the full phase diagram as a function of the coupling parameter , where and are the exchange amplitudes between respectively identical and different colors. The model corresponds to the infinite repulsion limit of the fermionic Hubbard model in the case, and of the bosonic Hubbard model for .
Three complementary methods have been used that give coherent results: linear flavor wave theory, iPEPS and quantum Monte Carlo. The resulting ground state phase diagram as a function of is given in Fig. 1 and includes 5 phases: three superfluid phases (USF, 2subSF, and 3subSF), a ferromagnetic phase (F), and a threesublattice colorordered phase (3subAF). The temperature dependent properties have been determined for and are presented in the phase diagram of Fig. 5. The 3subAF phase gives rise to a first order phase transition and the USF phase to a KT type phase transition.
The phase diagram for has many similarities with the one from the XXZ model (): A uniform superfluid and a sublattice superfluid characterized by phase angles, a ferromagnetic phase that can take different colors, and a sublattice antiferromagnetic, stripeordered phase. Each of these phases for can be seen as a generalization of the corresponding phase for . However, the case has one additional phase: the 2subSF phase, which is not strictly speaking a generalization of a phase of the case, but is in fact equivalent to the 2sublattice superfluid phase of the case (since only two colors are present in this phase). Thus it is interesting that for a 2color superfluid is competing with a 3color superfluid, and that both phases are stable in a certain parameter range. One finds also a competition between a 2sublattice and a 3sublattice color ordered state (which are degenerate at the classical level). In this case, however, it is the 3sublattice colorordered state which has always a lower energy than the 2sublattice state.
Finally, coming back to possible experimental realizations, three phases of this phase diagram should be accessible with cold atoms loaded in optical lattices: the USF phase and the F phase with bosonic particles, and the 3subAF phase at the symmetric AF point with fermionic particles.
v.1 Acknowledgment
We thank S. E. Korshunov for useful discussions about the superfluid transitions, and the Swiss National Fund for financial support.
References
 C. Wu, J.P. Hu, and S.C. Zhang, Phys. Rev. Lett. 91, 186402 (2003).
 C. Honerkamp and W. Hofstetter, Phys. Rev. Lett. 92, 170403 (2004).
 M. A. Cazalilla, A. F. Ho, and M. Ueda, New J. Phys. 11, 103033 (2009).
 T. Fukuhara, S. Sugawa, M. Sugimoto, S. Taie, and Y. Takahashi, Phys. Rev. A 79, 041604 (2009).
 A. V. Gorshkov, M. Hermele, V. Gurarie, C. Xu, P. S. Julienne, J. Ye, P. Zoller, E. Demler, M. D. Lukin, and A. M. Rey, Nat Phys 6, 289 (2010).
 S. Taie, Y. Takasu, S. Sugawa, R. Yamazaki, T. Tsujimoto, R. Murakami, and Y. Takahashi, Phys. Rev. Lett. 105, 190401 (2010).
 S. Sugawa, K. Inaba, S. Taie, R. Yamazaki, M. Yamashita, and Y. Takahashi, Nature Physics 7, 642 (2011).
 C. K. Law, H. Pu, and N. P. Bigelow, Phys. Rev. Lett. 81, 5257 (1998).
 A. B. Kuklov and B. V. Svistunov, Phys. Rev. Lett. 90, 100401 (2003).
 E. Altman, W. Hofstetter, E. Demler, and M. D. Lukin, New Journal of Physics 5, 113 (2003).
 T. A. Tóth, A. M. Läuchli, F. Mila, and K. Penc, Phys. Rev. Lett. 105, 265301 (2010).
 B. Bauer, P. Corboz, A. M. Läuchli, L. Messio, K. Penc, M. Troyer, and F. Mila, Phys. Rev. B 85, 125116 (2012).
 D. P. Arovas, Phys. Rev. B 77, 104404 (2008).
 M. Hermele and V. Gurarie, Preprint (2011), arXiv:1108.3862 .
 P. Corboz, K. Penc, F. Mila, and A. M. Läuchli, Phys. Rev. B 86, 041106 (2012a).
 H. H. Zhao, C. Xu, Q. N. Chen, Z. C. Wei, M. P. Qin, G. M. Zhang, and T. Xiang, Phys. Rev. B 85, 134416 (2012).
 P. Corboz, M. Lajkó, K. Penc, F. Mila, and A. M. Läuchli, arXiv:1302.1108 (2013).
 P. Corboz, A. M. Läuchli, K. Penc, M. Troyer, and F. Mila, Phys. Rev. Lett. 107, 215301 (2011a).
 P. Corboz, M. Lajkó, A. M. Läuchli, K. Penc, and F. Mila, Phys. Rev. X 2, 041013 (2012b).
 M. Hermele, V. Gurarie, and A. M. Rey, Phys. Rev. Lett. 103, 135301 (2009).
 G. Szirmai, E. Szirmai, A. Zamora, and M. Lewenstein, Phys. Rev. A 84, 011611 (2011).
 N. Papanicolaou, Nuclear Physics B 240, 281 (1984).
 N. D. Mermin and H. Wagner, Phys. Rev. Lett. 17, 1133 (1966).
 N. D. Mermin, Rev. Mod. Phys. 51, 591 (1979).
 J. M. Kosterlitz and D. J. Thouless, Journal of Physics C: Solid State Physics 6, 1181 (1973).
 D. R. Nelson, Phys. Rev. B 18, 2318 (1978).
 S. Korshunov and G. Uimin, Journal of Statistical Physics 43, 1 (1986).
 S. E. Korshunov, Phys. Rev. Lett. 88, 167007 (2002).
 M. Hasenbusch, A. Pelissetto, and E. Vicari, Journal of Statistical Mechanics: Theory and Experiment 2005, P12002 (2005).
 F. Verstraete and J. I. Cirac, Preprint (2004), arXiv:condmat/0407066 .
 J. Jordan, R. Orús, G. Vidal, F. Verstraete, and J. I. Cirac, Phys. Rev. Lett. 101, 250602 (2008).
 P. Corboz, S. White, G. Vidal, and M. Troyer, Phys. Rev. B 84, 041108 (2011b).
 P. Corboz, R. Orus, B. Bauer, and G. Vidal, Phys. Rev. B 81, 165104 (2010).
 G. Vidal, Phys. Rev. Lett. 91, 147902 (2003).
 R. Orús and G. Vidal, Phys. Rev. B 78, 155117 (2008).
 H. Jiang, Z. Weng, and T. Xiang, Phys. Rev. Lett. 101, 090603 (2008).
 T. Nishino and K. Okunishi, J. Phys. Soc. Japan 65, 891 (1996).
 R. Orús and G. Vidal, Phys. Rev. B 80, 094403 (2009).
 S. Singh, R. N. C. Pfeifer, and G. Vidal, Phys. Rev. B 83, 115125 (2011).
 B. Bauer, P. Corboz, R. Orús, and M. Troyer, Phys. Rev. B 83, 125106 (2011).
 B. Frischmuth, F. Mila, and M. Troyer, Phys. Rev. Lett. 82, 835 (1999).
 L. Messio and F. Mila, Phys. Rev. Lett. 109, 205306 (2012).
 N. Prokof’ev, B. Svistunov, and I. Tupitsyn, Journal of Experimental and Theoretical Physics Letters 64, 911 (1996).
 B. B. Beard and U.J. Wiese, Phys. Rev. Lett. 77, 5130 (1996).
 E. L. Pollock and D. M. Ceperley, Phys. Rev. B 36, 8343 (1987).
 D. R. Nelson and J. M. Kosterlitz, Phys. Rev. Lett. 39, 1201 (1977).
 K. Harada and N. Kawashima, Phys. Rev. B 55, R11949 (1997).