Quantum and Thermal Phase Transitions of the Triangular SU(3) Heisenberg Model under Magnetic Fields
Daisuke Yamamoto, Chihiro Suzuki, Giacomo Marmorini, Sho Okazaki, Nobuo Furukawa
QQuantum and Thermal Phase Transitions of the Triangular SU(3) Heisenberg Modelunder Magnetic Fields
Daisuke Yamamoto, Chihiro Suzuki, Giacomo Marmorini, Sho Okazaki, and Nobuo Furukawa
Department of Physics and Mathematics, Aoyama Gakuin University, Sagamihara, Kanagawa 252-5258, Japan (Dated: April 8, 2020)We study the quantum and thermal phase transition phenomena of the SU(3) Heisenberg modelon triangular lattice in the presence of magnetic fields. Performing a scaling analysis on large-size cluster mean-field calculations endowed with a density-matrix-renormalization-group solver, wereveal the quantum phases selected by quantum fluctuations from the massively degenerate classicalground-state manifold. The magnetization process up to saturation reflects three different magneticphases. The low- and high-field phases have strong nematic nature, and especially the latter isfound only via a nontrivial reconstruction of symmetry generators from the standard spin andquadrupolar description. We also perform a semi-classical Monte-Carlo simulations to show thatthermal fluctuations prefer the same three phases as well. Moreover, we find that exotic topologicalphase transitions driven by the binding-unbinding of fractional (half) vortices take place, due to thenematicity of the low- and high-field phases. Possible experimental realization with alkaline-earth-like cold atoms is also discussed.
Introduction. — In solid-state physics, lattice Hamilto-nians symmetric under the special unitary group of de-gree N = 2, denoted by SU(2), have been intensivelystudied since the electron – the main actor in solids –has two internal (spin) degrees of freedom. Higher de-gree of symmetry, or N >
2, can be accessed only withfine-tuning of parameters in some models, e.g., of spinliquid crystals [1–3] and transition metal oxides [4–7],or as a consequence of exotic emergent phases [8–10].However, recent advances in experiments with cold gasesof alkaline-earth(-like) atoms, such as
Yb [11–16] and Sr [17, 18], have provided a new platform and strongmotivation in studying the enhanced continuous sym-metry of SU( N > I ( I = 5 / / N ≤ I +1) extension of the Hubbard model [19] andits strong-coupling limit, namely the SU( N ) Heisenbergmodel [20, 21]. In such higher symmetric systems, theground states often form a massively (quasi-)degeneratemanifold. Therefore, of particular interest are the quan-tum and thermal fluctuations selecting one of the many-body states and the emergence of exotic phase transitionphenomena [23].The SU(3) Heisenberg model on triangular lattice hasbeen theoretically studied as a special symmetric pointof the spin-1 bilinear-biquadratic (BLBQ) model [24–26]. Since the number of colors ( N = 3) is compat-ible with the tripartite structure of the triangular lat-tice, the SU(3) Heisenberg model with antiferromagneticcouplings exhibits no (apparent) geometrical frustration,unlike the SU(2) case [27]. Indeed, the ground state isuniquely determined ( up to trivial degeneracy) to bea simple three-color three-sublattice order at the levelof the classical, mean-field (MF), analysis [24, 25], andit has been confirmed by numerical investigations [26]. Whereas the ground state may not be so exciting, theproperties under the presence of magnetic field remain aninteresting open problem since the MF analysis yields anaccidental continuous degeneracy [24]. It is thus requiredto study the degeneracy lifting induced by quantum andthermal fluctuations.In this work, we explore the effect of quantum andthermal fluctuations on the phase transition phenomenaof the triangular SU(3) Heisenberg model in magneticfields. High magnetic field experiments have been play-ing a central role in understanding the properties of mag-netic materials [28], one of the fundamental reasons be-ing that a magnetic field, in combination with latticegeometry, topological features, fluctuation effects, etc.,stimulates the emergence of a rich variety of nontrivialmagnetic states such as magnetization plateaus [29–31],nematic states [32, 33], and field-induced quantum spinliquids [31, 34]. This is naturally expected to occur forgeneral SU( N ) systems. In experiments with N -colorcomponents of cold atoms, the “magnetic field” in thesense of spin models could be implemented by introduc-ing a difference among the chemical potentials for differ-ent color components. The recent quantum simulationof the SU(2) Heisenberg model has demonstrated thatthe magnetic-field effect can indeed be obtained by in-troducing finite population imbalance between the twohyperfine states of Li [35].First, we employ the cluster mean-field plus scal-ing (CMF+S) method [36–38] with two-dimensional(2D) density matrix renormalization group (DMRG)solver [39] to reveal the quantum phases selected fromthe nontrivial classical ground-state manifold. We findthat the quantum order-by-disorder mechanism stabilizesthree different phases with the three-sublattice structuredepending on the field strength, until the system reachesthe magnetic saturation. Of particular interest is thelow-field phase, whose field-dependence exhibits nontriv- a r X i v : . [ c ond - m a t . s t r- e l ] A p r ial preservation of the uniform (total) scalar nematic or-der parameter at zero value, despite of having a finitevalue on each sublattice. In addition, the intermediate-field phase also preserves the total magnetization as well,resulting in a characteristic plateau structure in the mag-netization curve at 2 / semi-classical multi-color Monte-Carlo (sMC) simulations [40] by intro-ducing a “relaxation acceleration” technique, and discussthe thermal phase transition phenomena. In addition tothe stabilization of the same three phases by thermalfluctuations, we find particular topological phase transi-tions characterized by the binding-unbinding of fractional(half) vortices. The SU(3) Heisenberg model in magnetic fields. — TheSU(3) Heisenberg model is given byˆ H SU(3) = 2 J (cid:88) (cid:104) i,j (cid:105) (cid:88) A=1 , , ··· , ˆ T A i ˆ T A j ( J > , (1)where ˆ T A i = ˆ λ A i / S i = ( ˆ S xi , ˆ S yi , ˆ S zi ) for A = 1 , , Q i = ( ˆ Q x − y i , ˆ Q z i , ˆ Q xyi , ˆ Q yzi , ˆ Q xzi )for A = 4 , · · · , λ A i , instead of the standard Gell-Mann matrix basis. The quadrupolar operators are( ˆ S xi ) − ( ˆ S yi ) , √
3( ˆ S zi ) − / √
3, ˆ S xi ˆ S yi + ˆ S yi ˆ S xi , ˆ S yi ˆ S zi +ˆ S zi ˆ S yi , and ˆ S zi ˆ S xi + ˆ S xi ˆ S zi , respectively. In this spin-1representation, the Hamiltonian (1) is equivalent to theBLBQ model [23–26, 41, 42] with equal positive coeffi-cients, acting on spin states σ = − , , H SU(3) = J (cid:88) (cid:104) i,j (cid:105) (cid:16) ˆ S i · ˆ S j + ˆ Q i · ˆ Q j (cid:17) . (2)Below, we discuss the system under magnetic (Zee-man) fields: ˆ H ≡ ˆ H SU(3) + ˆ H Z with ˆ H Z = − H (cid:80) i ˆ S zi .The magnetic field explicitly breaks the SU(3) symme-try down to U(1) × U(1); specifically, the global rotationsaround the ˆ S z and ˆ Q z axes remain since (cid:80) i [ ˆ S zi , ˆ H ] = (cid:80) i [ ˆ Q z i , ˆ H ] = 0. Hereafter, we specifically write the tworotational symmetries as U(1) S z and U(1) Q z . Withinthe site-decoupling MF approximation [24], the specificspin and quadratic orders in the ground state for 0 3; red triangles), obtained by the CMF+S analysis.The MF value of M is plotted together (dashed line). Leftand right axes are shifted by 2/3. The inset shows cluster-sizescalings of the critical fields. Quantum order-by-disorder. — In order to discuss thelifting of the accidental degeneracy by quantum fluctu-ations and reveal the quantum orders, we perform theCMF+S calculations [36–38] with 2D DMRG solver [39].We employ a triangular-shaped cluster of N C sites, inwhich the quantum intersite correlations are treated ex-actly within the cluster, whereas the couplings with theoutside spins are replaced by MF interactions. Under thethree-sublattice ( µ = A , B , C) ansatz, the self-consistentequations (cid:104) ˆ S µ (cid:105) = N C (cid:80) i µ ∈ C (cid:104) Ψ N C | ˆ S i µ | Ψ N C (cid:105) and theanalogous expressions for (cid:104) ˆ Q µ (cid:105) are solved by calculat-ing the ground state of the N C -site cluster, | Ψ N C (cid:105) , with2D DMRG in an iterative way until convergence [39].The scaling parameter ζ ≡ N B / (3 N C ), with N B beingthe number of bonds inside the cluster, serves as anindicator of the extent to which quantum correlationsare taken into account , interpolating the classical (MF; ζ = 0) and exactly-quantum ( N C → ∞ ; ζ = 1) lim-its. Here we perform the calcuations for N C = 10 , , ζ = 3 / , / , / 7) and make the linear extrapolation ofthe results towards ζ → N C = 36( ζ = 7 / 9) is also considered for the determination of thephase boundaries (see the inset of Fig. 1).We plot the quantum magnetization curves M ( H ) ≡ (cid:80) µ (cid:104) ˆ S zµ (cid:105) / H c = 3 . 40 and H c = 4 . 38 in the range of0 < H < H s . The low-field (LF) phase is characterizedby (cid:104) ˆ S zA (cid:105) = (cid:104) ˆ S zB (cid:105) (cid:54) = (cid:104) ˆ S zC (cid:105) ≈ (cid:104) ˆ Q x − y A (cid:105) = −(cid:104) ˆ Q x − y B (cid:105) , and (cid:104) ˆ Q x − y C (cid:105) = 0, modulo a global rotation in the ( Q x − y , Q xy ) plane and sublattice exchanges; the other compo-nents are all zero [see Fig. 2(a)]. Although the spin sec-tor ( S x , S y , S z ) forms a collinear structure along the field FIG. 2: (a) Nonzero components of the spin and quadrupolarmoments, obtained by the CMF+S analysis in a fixed gaugewith (cid:104) ˆ Q xyA (cid:105) = (cid:104) ˆ Q yzA (cid:105) = 0. The inset shows the three-sublatticestructure. (b) Spherical plots of |(cid:104) S | ˆ Q x − y | S (cid:105)| and its π/ π rotations about U (1) S z with | S (cid:105) being the spin coherentstate pointing in the S direction [23]. axis, the transverse quadrupolar moments ( Q x − y , Q xy )break the rotational symmetry around S z . It is partic-ularly interesting that the a π -rotation around the S z axis is sufficient for ( Q x − y , Q xy ) to return the initialstate as illustrated in Fig. 2(b) due to the nematic na-ture, reflecting the factor 2 in the commutation relation[ ˆ Q x − y , ˆ Q xy ] = 2 i ˆ S z . Thus, it is concluded that the LFphase breaks the (U(1) S z / Z ) × Z (i.e., half of the orig-inal rotational and three-fold translational) symmetries.Consequently, the remaining U(1) Q z symmetry guaran-tees the preservation of the uniform nematic scalar or-der parameter S u ≡ (cid:80) µ (cid:104) ˆ Q z µ (cid:105) / 3, resulting in the plateauformation at zero value in the field dependence (Fig. 1),whereas each sublattice value (cid:104) ˆ Q z µ (cid:105) is finite [Fig. 2(a)].At H = H c , the transverse quadrupolar moments van-ishes and the U(1) S z symmetry is restored. Thus, inthe intermediate-field (IF) phase, both M and S u ex-hibit plateau behavior in the range of H c < H < H c .The longitudinal spin moments (cid:104) ˆ S zµ (cid:105) have the valuesof approximately (1 , , 0) (not exactly, due to quantumdepletion) and thus M = 2 / 3. Such a plateau for-mation has been reported in the spin-1 BLBQ modelwhen the quadrapolar coupling is larger than the dipolarone [24]. Our results showed that the plateau is stabi- lized by purely quantum effects even for equal bilinear-biquadratic (SU(3)-symmetric) coupling.In the high-field (HF) phase, the spin ( S x , S y , S z ) sec-tor forms a “2:1” structure of the V shape, similar tothe SU(2) case [29]. Therefore, it apparently seems to bea standard non-nematic spin order. However, we noticethat the curves of M and S u / √ S z × U(1) Q z :the linear combination of generators ˆ P z + ≡ ˆ S z + √ ˆ Q z is broken, while ˆ P z − ≡ √ ˆ S z − ˆ Q z is preserved. TheU(1) P z ± action produces a rotation of the system in theplane of ˆ P x ± ≡ ( ˆ S x ± ˆ Q xz ) / √ P y ± ≡ ( ˆ S y ± ˆ Q yz ) / √ (cid:104) ˆ P x − ,µ (cid:105) = (cid:104) ˆ P y − ,µ (cid:105) = 0 in theHF phase, which indicates the preservation of the U(1) P z − symmetry. As for the broken U(1) P z + , a π -rotation is suf-ficient for ( ˆ P x + , ˆ P y + ) to return to the initial state since[ ˆ P x + , ˆ P y + ] = 2 i ˆ P z + , and thus the HF phase possesses a ne-matic nature despite of the apparent spin (dipolar) order.Considering also the sublattice exchange, we concludethat the HF phase breaks (U(1) P z + / Z ) × Z . Thermal phase diagram. — Given the strong nematicnature of the zero-temperature phases, it is interestingto study the thermal phase transitions, especially asso-ciated with the (U(1)/ Z ) × Z symmetry breaking. Tothis end, we employ the sMC simulations [40] within thedirect-product approximation: | Ψ cl (cid:105) = ⊗ i | ψ i (cid:105) with lo-cal wave functions | ψ i (cid:105) = (cid:80) σ d i,σ | σ i (cid:105) ( | d i | = 1). Thestandard Metropolis updates are performed for the co-efficients d i,σ on L × L rhombic clusters under periodicboundary conditions, based on the Boltzmann distribu-tion p ∝ exp( − E cl /k B T ) with E cl ≡ (cid:104) Ψ cl | ˆ H| Ψ cl (cid:105) [40]. Wefurther develop the sMC method by applying a “relax-ation acceleration” (RA) with local unitary transforma-tions e ic ˆ H loc i | ψ i (cid:105) , where c are uniformly distributed ran-dom numbers and ˆ H loc i ≡ ( ⊗ j (cid:54) = i (cid:104) ψ j | ) ˆ H ( ⊗ j (cid:54) = i | ψ j (cid:105) ). Herewe choose, after some trials, | c | ≤ π (cid:107) ˆ H loc i (cid:107) F with (cid:107) · · · (cid:107) F being the Frobenius norm. The RA sweeps over latticesites are performed twice following each Metropolis up-date sweep. This method, applied to highly-symmetricsystems, is significantly more efficient in improving decor-relation and avoiding trapping in local minima. See theSM [43] for more technical details and the performanceexamination of the RA procedure.Figure 3 shows the thermal phase diagram obtained bythe sMC, which is reliable in the region away from thelow-temperature quantum regime, since it neglects theinter-site quantum correlations. It is seen that the samethree (LF, IF, and HF) phases are selected also by ther-mal fluctuations from the classical degenerate manifoldsat T = 0. The boundaries to the paramagnetic phase aredetermined from the divergence of the correlation lengthof the longitudinal ( S z , Q z ; indicated by square sym- FIG. 3: Thermal phase diagram obtained by the sMC simu-lations. We also mark the critical fields H c1 and H c2 obtainedby the CMF+S method at the quantum ( T = 0) limit. Thedashed lines are the sketches of the phase boundaries expectedfrom the combination of the sMC (valid at high temperatures)and CMF+S (valid at T = 0) results. bols) and transverse (all the others; indicated by trian-gles and hexagons) components of the three-sublatticeorder parameters. The LF-IF and IF-HF boundaries aredetermined by the scaling analyses of the susceptibilityin the corresponding transverse components.We show in Fig. 4(a) the stiffness ρ S z ( T ) for atwist of the spin and quadrupolar moments around S z in the transition region between LF and IF. It canbe clearly seen that ρ S z ( T ) at the transition point T = T c does not satisfy the standard universal rela-tion ρ S z ( T c ) = 2 k B T c /π for the Kosterlitz-Thouless (KT)transitions [44]. This is attributed to the nematic natureof ( Q x − y , Q xy ), which break U(1) / Z rotations around S z [shown in Fig. 2(b)]. Due to this, the ( Q x − y , Q xy )moments can form a topologically stable vortex with frac-tional vorticity ρ v = 1 / ρ v = 1. This half-vortex is analogous to the 180 degree disclination of ne-matic liquid crystals [45]. The transition from LF to IFis associated with the unbinding of pairs of half-vortexand half-antivortex, resulting in the modified universalrelation ρ S z ( T c ) = 2 k B T c /πρ v = 8 k B T c /π [46]. Thisparticular topological transition takes place also at theboundary of the HF and IF (or paramagnetic) phase[Fig. 4(c)], where it is related to the U(1) / Z rotationaround P z + mentioned above. This universal jump is as-sociated with the unbinding of half-vortex pairs in the( P x + , P y + ) plane. Note that we confirmed that the stiff-ness for a twist around P z − is always zero (within thenumerical error) in the HF phase.Let us comment briefly on the limit of zero mag-netic field. Since the classical ground state is given by d i = (1 , , , , , , 1) for sublattice A, B,and C, respectively, or SU(3) rotations thereof [25], thesymmetry is spontaneously broken down to U(1) × U(1). FIG. 4: (a) Stiffness ρ S z ( T ) along H/J = 2 . 0, which showsa universal jump ρ S z ( T c ) = 8 k B T c /π at the LF-IF transition,except for a slight finite-size effect. (b) Vortex and antivortexwith half vorticity ρ v = ± / Q x − y , Q xy )plane. The right panel is a schematic illustration of a topolog-ical half-vortex pair excitation on the background of a uniform quadrupolar order on, say, sublattice A. (c) Same as in (a) for ρ P z + ( T ) at the HF-paramagnetic transition along H/J = 3 . The fundamental group π (SU(3)/U(1) × U(1)) is triv-ial [47] and therefore there are no vortex-induced finitetemperature phase transitions [44, 48]. The tendencyof the IF-paramagnetic line towards the ( T, H ) = (0 , Experimental realization. — In optical-lattice experi-ments with alkaline-earth(-like) atoms, the field strength H corresponds to the difference in the chemical poten-tials of each component: µ − µ = µ − µ − = H . Thiscould be prepared via the population imbalance in dif-ferent colors [35]. More directly, if one could introducea state-dependent gradient potential, say, in the x di-rection, V ext σ ( x ) = σV x , the magnetization process in − H s < H < H s would be realized as the spatial depen-dence for local magnetic field H ( x ) = 2 V x (in the senseof the local density approximation [49]) with no over-all population imbalance. The formation of the three-sublattice orders can be detected by the time-of-flightimage of the momentum distribution [35, 50–54], and theIF state would appear as a spatial plateau. In addition,the extension of the quantum-gas microscope techniqueto SU( N ) systems could provide a wealth of detailed mea-surements of each phase found here. Conclusions. — We studied the quantum and thermalphase transition phenomena of the SU(3) Heisenbergmodel under magnetic fields by using the CMF+S andsMC methods. We demonstrated that pure quantum-fluctuation effects stabilize a magnetization plateau at2/3 of the saturation in the intermediate range of thefield strength. The uniform scalar nematic order param-eter also forms a plateau at zero value, which, more inter-estingly, appears already in the lower-field phase with nomagnetization plateau. This reflects the preservation ofadditional U(1) symmetry around Q z unique to the sys-tem with the SU(3) symmetry (reduced to U(1) × U(1)under magnetic fields). The high-field phase exhibitsan unexpected nematic nature steming from nontrivialpartial breaking of U(1) × U(1) symmetry. Moreover, thestrong nematic nature of the low- and high-field phasesgives rise to fractional vortices and antivortices, whosepair dissociation results in a topological phase transitionwith vorticity ρ v = 1 / ≈ . J/k B at the highest value), providea solid guideline for future experiments with alkaline-earth(-like) atoms. Besides, the physics we explored isrelevant to solid-state materials with nearly SU(3) sym-metric parameters and, more generally, to systems withmultipolar orders. Finally, it should be mentioned thatthe calculation methods used here can be applied to nu-merous types of classical and quantum simulations; es-pecially, the acceleration and decorrelation of the sMCupdates by the RA technique are expected to become in-creasingly important for models with higher symmetry,such as SU( N ≥ [1] S. Nakatsuji, Y. Nambu, H. Tonomura, O. Sakai, S. Jonas,C. Broholm, H. Tsunetsugu, Y. Qiu, and Y. Maeno, Sci-ence , 1697 (2005).[2] H. Tsunetsugu and M. Arikawa, J. Phys. Soc. Jpn. ,083701 (2006).[3] S. Bhattacharjee, V. B. Shenoy and T. Senthil, Phys. Rev.B , 092406 (2006).[4] K. I. Kugel and D. I. Khomskii, Sov. Phys. JETP , 725(1973).[5] D. P. Arovas and A. Auerbach, Phys. Rev. B , 10114(1995).[6] Y. Q. Li, M. Ma, D. N. Shi, and F. C. Zhang, Phys. Rev.Lett. , 3527 (1998).[7] Y. Tokura and N. Nagaosa, Science , 462 (2000).[8] A. J. Keller, S. Amasha, I. Weymann, C. P. Moca, I.G. Rau, J. A. Katine, H. Shtrikman, G. Zar´and, and D.Goldhaber-Gordon, Nature Physics , 145 (2014).[9] P. Chen, Z.-L. Xue, I. P. McCulloch, M.-C. Chung, C.-C. Huang, and S.-K. Yip, Phys. Rev. Lett. , 145301(2015).[10] V. L. Quito, Pedro L. S. Lopes, Jose A. Hoyos, E. Mi-randa, Eur. Phys. J. B , 17 (2020).[11] T. Fukuhara, Y. Takasu, M. Kumakura, and Y. Taka-hashi, Phys. Rev. Lett. , 030401 (2007). [12] M. A. Cazalilla, A. F. Ho, and M. Ueda, New J. Phys. ,103033 (2009).[13] H. Hara, Y. Takasu, Y. Yamaoka, J. M. Doyle, and Y.Takahashi, Phys. Rev. Lett. , 205304 (2011).[14] S. Taie, R. Yamazaki, S. Sugawa, and Y. Takahashi,Nat.Phys. , 825 (2012).[15] C. Hofrichter, L. Riegger, F. Scazza, M. H¨ofer, D. R. Fer-nandes, I. Bloch, and S. F¨olling, Phys. Rev. X , 021030(2016).[16] H. Ozawa, S. Taie, Y. Takasu, and Y. Takahashi, Phys.Rev. Lett. , 225303 (2018).[17] B. J. DeSalvo, M. Yan, P. G. Mickelson, Y. N. Mar-tinez de Escobar, and T. C. Killian, Phys. Rev. Lett. ,030402(2010).[18] M. K. Tey, S. Stellmer, R. Grimm, and F. Schreck, Phys.Rev. A , 011608(R) (2010).[19] C. Honerkamp and W. Hofstetter, Phys. Rev. Lett. ,170403 (2004).[20] 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. , 289 (2010).[21] P. Nataf and F. Mila, Phys. Rev. Lett. , 127204(2014)[22] J. S. Gardner, M. J. P. Gingras, and J. E. Greedan, Rev.Mod. Phys. , 53 (2010).[23] Introduction to Frustrated Magnetism , edited by C.Lacroix, P. Mendels, and F. Mila (Springer-Verlag, Berlin,2011).[24] A. L¨auchli, F. Mila, and K. Penc, Phys. Rev. Lett. ,087205 (2006).[25] A. Smerald and N. Shannon, Phys. Rev. B , 184430(2013).[26] B. Bauer, P. Corboz, A. M. L¨auchli, L. Messio, K. Penc,M. Troyer, F. Mila, Phys. Rev. B , 125116 (2012).[27] R. Moessner and A. P. Ramirez, Physics Today , 24(2006).[28] High Magnetic Fields: Applications in Condensed Mat-ter Physics and Spectroscopy , edited by C. Berthier, L.-P.Levy, and G.Martinez (Springer-Verlag, Berlin, 2002).[29] A. V. Chubukov and D. I. Golosov, J. Phys.: Condens.Matter , 69 (1991).[30] Y. Shirata, H. Tanaka, A. Matsuo, and K. Kindo, Phys.Rev. Lett. , 057205 (2012).[31] S. Nishimoto, N. Shibata, and C. Hotta, Nat. Commun. , 2287 (2013).[32] K. Nawa, M. Takigawa, M. Yoshida, and K. Yoshimura,J. Phys. Soc. Jpn. , 094709 (2013).[33] N. B¨uttgen, K. Nawa, T. Fujita, M. Hagiwara, P. Kuhns,A. Prokofiev, A. P. Reyes, L. E. Svistov, K. Yoshimura,and M. Takigawa, Phys. Rev. B , 134401 (2014).[34] S.-H. Baek, S.-H. Do, K.-Y. Choi, Y. S. Kwon, A. U. B.Wolter, S. Nishimoto, J. van den Brink, and B. B¨uchner,Phys. Rev. Lett. , 037201 (2017).[35] P. T. Brown, D. Mitra, E. Guardado-Sanchez, P. Schauß,S. S. Kondov, E. Khatami, T. Paiva, N. Trivedi, D. A.Huse, W. S. Bakr, Science , 1385 (2017).[36] D. Yamamoto, A. Masaki, and I. Danshita, Phys. Rev.B , 054516 (2012).[37] D. Yamamoto, G. Marmorini, and I. Danshita, Phys.Rev. Lett. , 127203 (2014); ibid , , 259901 (2014).[38] D. Yamamoto, H. Ueda, I. Danshita, G. Marmorini, T.Momoi, and T. Shimokawa, Phys. Rev. B , 014431(2017).[39] D. Yamamoto, G. Marmorini, M. Tabata, K. Sakakura, and I. Danshita, Phys. Rev. B , 140410(R) (2019).[40] E. M. Stoudenmire, S. Trebst, and L. Balents, Phys. Rev.B , 214436 (2009).[41] T. A. Toth, A. M. L¨auchli, F. Mila, K. Penc, Phys. Rev.B , 140403(R) (2012).[42] I. Niesen and P. Corboz, Phys. Rev. B , 180404(R)(2017); ibid . , 245146 (2018).[43] See Supplemental Material at http:[URL] for detailedstructure of the degenerate classical ground-state mani-fold and technical details of the sMC simulations.[44] J. M. Kosterlitz, Rep. Prog. Phys. , 026001 (2016).[45] N. D. Mermin, Rev. Mod. Phys. , 591 (1979).[46] S. E. Korshunov Phys. Rev. B , 054416 (2002).[47] H. T. Ueda, Y. Akagi, and N. Shannon, Phys. Rev. A ,021606(R) (2016).[48] H. Kawamura and S. Miyashita, J. Phys. Soc. Jpn. , 023701 (2010).[49] S. Bergkvist, P. Henelius, and A. Rosengren, Phys. Rev.A , 053601 (2004).[50] M. F. Parsons, A. Mazurenko, C. S. Chiu, G. Ji, D. Greif,and M. Greiner, Science , 1253 (2016).[51] M. Boll, T. A. Hilker, G. Salomon, A. Omran, J. Nespolo,L. Pollet, I. Bloch, and C. Gross, Science , 1257 (2016).[52] L. W. Cheuk, M. A. Nichols, K. R. Lawrence, M. Okan,H. Zhang, E. Khatami, N. Trivedi, T. Paiva, M. Rigol,and M. W. Zwierlein, Science , 1260 (2016).[53] A. Mazurenko, C. S. Chiu, G. Ji, M. F. Parsons, M.Kan´asz-Nagy, R. Schmidt, F. Grusdt, E. Demler, D. Greif,and M. Greiner, Nature (London) , 462 (2017).[54] T. A. Hilker, G. Salomon, F. Grusdt, A. Omran, M.Boll, E. Demler, I. Bloch, and C. Gross, Science , 484(2017). Supplementary Material for “Quantum and Thermal Phase Transitions of the Triangular SU(3) HeisenbergModel under Magnetic Fields”Classical degeneracy manifold of the triangular SU(3) Heisenberg model with magnetic fields Within the site-decoupling mean-field (MF) approximation, the ground state is assumed to be a direct product oflocal wave-functions ( | d i | = 1): | Ψ cl (cid:105) = ⊗ i | ψ i (cid:105) with | ψ i (cid:105) = (cid:88) σ d i,σ | σ i (cid:105) . (S1)The coefficient vector d i = ( d i, − , d i, , d i, ) normalized to unit length ( | d i | = 1) identifies the local state at site i as asuperposition of the three basis states ( σ = − , , µ = A , B , C) ansatz, the variationalenergy E cl ≡ (cid:104) Ψ cl | ˆ H| Ψ cl (cid:105) can be written as E cl N = J (cid:18) λ A + λ B + λ C − H J (cid:19) − h J − J, (S2)where N is the number of sites, λ µ ≡ (cid:104) Ψ cl | ˆ λ i µ | Ψ cl (cid:105) is an eight-component classical vector of length (cid:112) / 3, and H = (0 , , H, , , , , E cl is simply achieved when λ µ = 2 H / J for h ≤ J . The overlinemeans the average over µ = A , B , C. For H > J , the two conditions S zµ = 2 H/ J and Q z µ = 0, cannot besimultaneously satisfied because Q z µ / √ S zµ − / S zµ = H/ J + 3 / Q z µ = H/ J − / J < H < H s , with H s = 9 J being thesaturation field.From the above discussion, the magnetization M ≡ (cid:80) i (cid:104) ˆ S zi (cid:105) /N = S zµ is uniquely determined as shown in Fig. 1(dashed line). However, the specific spin and quadratic orders remain massively degenerate because the number ofconditions is smaller than that of variational parameters d i µ ,σ . Technical details of the semi-classical Monte Carlo analysis In the main text, we employ the semi-classical Monte Carlo (sMC) simulations [1] on L × L rhombic clusters underperiodic boundary conditions, since the fully-quantum Monte Carlo method suffers from the so-called sign problemfor frustrated quantum systems. First, we assume that the wave function of the entire system is described as a directproducts of local wave functions as in the MF approximation [Eq. (S1)], although the three-sublattice ansatz is notassumed. The total energy of the system is given by E cl ( { d i } ) = (cid:104) Ψ cl | ˆ H| Ψ cl (cid:105) = J (cid:88) (cid:104) i,j (cid:105) (cid:18) | d ∗ i · d j | − (cid:19) − H (cid:88) i (cid:0) | d i, | − | d i, − | (cid:1) (S3)within the direct-product approximation. We first set the initial values of d i on the entire lattice sites to L × L complexrandom vectors distributed homogeneously on the sphere of radius one in 3 (real) +3 (imaginary) dimensions. Startingwith the initial state, we perform the standard Metropolis local updates of d i to generate a sequence of states weightedby the probability proportional to the Boltzmann factor exp( − E cl ( { d i } ) /k B T ). Typical simulations contain 10 and2 × Monte Carlo steps for the thermalization of the state and the samplings of physical quantities, respectively.One Monte Carlo step consists of one Metropolis sweep over all sites followed by two “relaxation acceleration” (RA)sweeps (which will be explained in Sec. C).The quantum-mechanical expectation values of the local spin and quadrupolar moments can be calculated by λ Ai ≡ (cid:104) ψ i | ˆ λ Ai | ψ i (cid:105) = (cid:88) σ,σ (cid:48) (cid:104) σ i | ˆ λ Ai | σ (cid:48) i (cid:105) d ∗ i,σ d i,σ (cid:48) (S4)for a given site with vector d i . The eight components of the vector ˆ λ i correspond to the spin components ( ˆ S xi , ˆ S yi , ˆ S zi )for A = 1 , , Q x − y i , ˆ Q z i , ˆ Q xyi , ˆ Q yzi , ˆ Q xz ) for A = 4 , , · · · , 8, respectively, as in themain text. To discuss the spontaneous symmetry breaking, we calculate the correlation lengths of the diagonal andtransverse components: ξ (cid:107) (1 , = √ L π (cid:118)(cid:117)(cid:117)(cid:116) S (cid:107) (1 , ( Q K ) S (cid:107) (1 , ( Q K + (0 , π/ √ L )) − ξ ⊥ (1 , = √ L π (cid:118)(cid:117)(cid:117)(cid:116) S ⊥ (1 , ( Q K ) S ⊥ (1 , ( Q K + (0 , π/ √ L )) − S (cid:107) (1) ( k ) = 1 L (cid:88) i,j (cid:104)(cid:104) S zi S zj (cid:105)(cid:105) T e − i k · ( r i − r j ) , S (cid:107) (2) ( k ) = 1 L (cid:88) i,j (cid:104)(cid:104) Q z i Q z j (cid:105)(cid:105) T e − i k · ( r i − r j ) , S ⊥ (1) ( k ) = 1 L (cid:88) i,j (cid:104)(cid:104) Q x − y i Q x − y j + Q xyi Q xyj (cid:105)(cid:105) T e − i k · ( r i − r j ) , and S ⊥ (2) ( k ) = 1 L (cid:88) i,j (cid:104)(cid:104) S xi S xj + S yi S yj + Q yzi Q yzj + Q xzi Q xzj (cid:105)(cid:105) T e − i k · ( r i − r j ) . (S6)Here, (cid:104)(cid:104)· · · (cid:105)(cid:105) T means the thermal average in terms of Monte Carlo samplings and the ordering vector k = Q K ≡ (4 π/ , 0) corresponds to the three-sublattice order shown as the inset of Fig. 2(a).The stiffness ρ S z ( T ) for a twist generated by the unitary transformation, ˆ U S z ( q ) ≡ exp[ iq (cid:80) i x i ˆ S zi ], is defined inthe standard way as the second derivative of the free energy per unit area with respect to the twist angle q : ρ S z ( T ) = 1 L ∆ S (cid:42)(cid:42) ∂ (cid:104) ˆ U S z ( q ) ˆ H ˆ U † S z ( q ) (cid:105) ∂q (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) q =0 (cid:43)(cid:43) T − k B T (cid:42)(cid:42)(cid:32) ∂ (cid:104) ˆ U S z ( q ) ˆ H ˆ U † S z ( q ) (cid:105) ∂q (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) q =0 (cid:33) (cid:43)(cid:43) T = − √ L (cid:42)(cid:42) J (cid:88) (cid:104) i,j (cid:105) ( x i − x j ) (cid:16) S xi S xj + S yi S yj + 4( Q x − y i Q x − y j + Q xyi Q xyj ) + Q yzi Q yzj + Q xzi Q xzj (cid:17)(cid:43)(cid:43) T − √ L k B T (cid:42)(cid:42) J (cid:88) (cid:104) i,j (cid:105) ( x i − x j ) (cid:16) S xi S yj − S yi S xj + 2( Q x − y i Q xyj − Q xyi Q x − y j ) − ( Q yzi Q xzj − Q xzi Q yzj ) (cid:17)(cid:35) (cid:43)(cid:43) T , (S7)where ∆ S = √ / x -axis, althoughthe value of ρ S z ( T ) does not depend on this choice for L → ∞ . In a similar way, the stiffness ρ P z + ( T ) regardingˆ U P z + ( q ) ≡ exp[ iq (cid:80) i x i ˆ P z + ,i ] with ˆ P z + ,i ≡ ˆ S zi + √ ˆ Q z i is defined as ρ P z + ( T ) = − √ L (cid:42)(cid:42) J (cid:88) (cid:104) i,j (cid:105) ( x i − x j ) (cid:16) P x + ,i P x + ,j + P y + ,i P y + ,j ) + P x − ,i P x − ,j + P y − ,i P y − ,j + Q x − y i Q x − y j + Q xyi Q xyj (cid:17)(cid:43)(cid:43) T − √ L k B T (cid:42)(cid:42) J (cid:88) (cid:104) i,j (cid:105) ( x i − x j ) (cid:0) P x + ,i P y + ,j − P y + ,i P x + ,j ) − ( P x − ,i P y − ,j − P y − ,i P x − ,j )+ Q x − y i Q xyj − Q xyi Q x − y j (cid:17)(cid:35) (cid:43)(cid:43) T , (S8)where ˆ P x ± ≡ ( ˆ S x ± ˆ Q xz ) / √ P y ± ≡ ( ˆ S y ± ˆ Q yz ) / √ Numetical data of the sMC simulations Here we present some numerical data for the standard scaling analyses performed to determine the phase boundariesof Fig. 3 of the main text. All the three (LF, IF, and HF) phases possess a three-sublattice order in the diagonal FIG. S1: Typical examples of the numerical data for the scaling analyses that determine the finite-temperature phase diagramshown in Fig. 3 of the main text. The error bars assigned to each data point are estimated as the square root of variance ofabout 10 independent sMC simulations. The shaded bands indicate the estimated transition points with their error bar. components, S z and Q z . Therefore, the transition points to the paramagnetic phase can be identified by thedivergence of ξ (cid:107) (1 , . Figure S1 (a) shows a typical example of the transition from the IF to paramagnetic phase. Thecurves of the scaled correlation length ξ (cid:107) (1 , /L for different linear sizes L cross each other at a critical point, within theerror bar estimated from the square root of the variance over about 10 sMC simulations. We plot in Fig.3 the crossingpoints in ξ (cid:107) (2) /L , which are slightly less size-dependent, as the phase boundary from the ordered to the paramagneticstates.Note that, in the standard Kosterlitz-Thouless transition of the 2D XY model with rotational symmetry, the scaledcorrelation length of the transverse (XY) components does not exhibit an isolated critical (crossing) point but a finitecritical range from T = 0 with a constant value independent of L . Interestingly, the scaled correlation length of thetransverse components, ξ ⊥ (1) /L ( ξ ⊥ (2) /L ) in the present case shows a crossing behavior in the vicinity of the transitionsfrom the LF (HF) to parmagnetic transitions, in spite of the continuous nature of the rotational symmetries of thesystem around the S z and Q z axes [see an example for the LF-paramagnetic transition in Fig. S1(b)]. This maybe attributed to the combined effect of the simultaneous discrete (diagonal) and continuous (transverse) symmetrybreakings. A similar (apparent) crossing behavior of the scaled correlation length for the transverese componentshas been reported in previous studies on some related 2D models with combined discrete and continuous symmetrybreakings [2]. The crossing points in the scaled correlation lengths of the diagonal and transverse components arelocated at almost the same position (within the error bar) as seen in Fig. 3 of the main text.From the LF to the IF (HF to IF) phase, the topological transition associated with the unbinding of pairs of half-vortex and half-antivortex in the plane of Q x − y and Q xy ( P x + and P y + ) occurs, as explained in the main text. In thiscase, the corresponding scaled correlation length does not exibit an isolated critical point. Therefore, to locate thetopological transition points, we perform the scaling analysis on the susceptibilities of the corresponding quantities: χ ⊥ (1) = Jk B T S ⊥ (1) ( Q K ) and χ ⊥ (2) = Jk B T S ⊥ (2) ( Q K ) , (S9)which obey the following scaling relations: χ ⊥ (1 , = L − η ˜ χ ⊥ (1 , ( tL /ν ) (S10)with unknown universal functions ˜ χ ⊥ (1 , of t = ( T − T c ) /T c . Here, η and ν are the correlation function and correlationlength critical exponents, respectively. At the LF-to-IF (HF-to-IF) topological phase transition, χ ⊥ (1) ( χ ⊥ (2) ) is expectedto scale with the exact Kosterlitz-Thouless exponent η = 1 / L η − χ ⊥ (1 , become size-independent at the corresponding transition points t = 0 with η = 1 / 4. Figure S1(c) and S1(d) show0 FIG. S2: Comparison of the results with and without RA. The values of the scaled correlation length ξ (cid:107) (2) /L for H/J = 0 . L = 96 obtained by the sMC simulations (i) for 10 samples without RA (black), (ii) for 3 × samples without RA(red), and (iii) for 10 samples with RA (blue) are plotted. The error bars are estimated from the square root of the varianceover nine independent sMC simulations. The inset shows the energy of the system per site for the first 200 Monte Carlo stepsin the thermalization processes with (blue) and without (red) RA. typical examples of the scaling analysis performed to determine the LF-IF and HF-IF boundaries, respectively, plottedin Fig. 3. Relaxation acceleration techniques