Anisotropic-Exchange Magnets on a Triangular Lattice: Spin Waves, Accidental Degeneracies, and Dual Spin Liquids

We present an extensive overview of the phase diagram, spin-wave excitations, and finite-temperature transitions of the anisotropic-exchange magnets on an ideal nearest-neighbor triangular lattice. We investigate transitions between five principal classical phases of the corresponding model: ferromagnetic, N\'{e}el, its dual, and the two stripe phases. Transitions are identified by the spin-wave instabilities and by the Luttinger-Tisza approach. Some of the transitions are direct and others occur via intermediate phases with more complicated forms of ordering. In a portion of the N\'{e}el phase, we find spin-wave instabilities to a long-range spiral-like state. In the stripe phases, quantum fluctuations are mostly negligible, leaving the ordered moment nearly saturated even for the $S=1/2$ case. However, for a two-dimensional surface of the full 3D parameter space, the spin-wave spectrum in one of the stripe phases exhibits an enigmatic accidental degeneracy manifested by pseudo-Goldstone modes. As a result, despite the nearly classical ground state, the ordering transition temperature in a wide region of the phase diagram is significantly suppressed from the mean-field expectation. We identify this accidental degeneracy as due to an exact correspondence to an extended Kitaev-Heisenberg model with emergent symmetries that naturally lead to the pseudo-Goldstone modes. There are previously studied dualities within the Kitaev-Heisenberg model on the triangular lattice that are exposed here in a wider parameter space. One important implication of this correspondence for the $S=1/2$ case is the existence of a region of the spin-liquid phase that is dual to the spin-liquid phase discovered recently by us. We complement our studies by the density-matrix renormalization group of the $S=1/2$ model to confirm some of the duality relations and to verify the existence of the dual spin-liquid phase.


I. INTRODUCTION
Ever since the seminal works by Wannier [1] and Anderson [2], a motif of spins on a triangular-lattice network epitomizes the idea of geometric frustration that can give rise to non-magnetic spin-liquid states [3,4].A variety of materials and models with the triangular, kagomé, and pyrochlore lattices have provided a natural playground for geometric frustration and realized various exotic and quantum-disordered states [2][3][4][5][6][7].
More recently, magnetic materials with anisotropic spin-spin interactions, which arise from spin-orbit coupling in their magnetic ions, have offered a different path to achieve similar goals.A strong mixing of spin and orbital degrees of freedom leads to the bond-dependent anisotropic-exchange interactions providing an alternative mechanism for frustration [8,9].A particular case is that of the so-called compass model [10] on the tricoordinated honeycomb lattice with each spin component interacting selectively along only one of the bonds via an Ising-like interaction.In a celebrated work [11], Kitaev showed that it has a spin-liquid ground state with fractionalized excitations, a finding that set off a significant research effort.In real materials, however, desired terms occur along with the other diagonal and off-diagonal components of the anisotropic-exchange matrix that are allowed by the lattice symmetry [12][13][14][15][16].These terms have proven to be detrimental to the Kitaev spin liquid and so far have prevented its definite realization [12].
In our first study, Ref. [31], we argued that the experimental range of parameters of the model that should describe a disorder-free YMGO does not support a spinliquid state, in agreement with exact-diagonalization [37] and variational Monte Carlo studies [38].In our more recent work, Ref. [39], we provided a detailed exploration of the phase diagram of the most generic nearest-neighbor triangular-lattice model in order to find out whether anisotropic-exchange interactions on this lattice can po-tentially create much desired exotic states.We have used the density-matrix renormalization group (DMRG) for the S = 1/2 model and discovered a spin-liquid region in its 3D phase diagram [39].We have also established a close similarity and a direct isomorphism of this newly found spin-liquid phase to the much studied spin liquid of the fully isotropic Heisenberg J 1 -J 2 model on the triangular lattice [46][47][48][49][50][51].We have also pointed out that the spin liquid with open spinon Fermi surface [20,21,25,34] is not realized in the phase diagram of the model.This is also in accord with Ref. [38].
In the present study, we expand our previous work in several directions.First, we provide a quasiclassical description of the five principal magnetically-ordered single-Q phases that span the 3D phase diagram of the nearestneighbor anisotropic-exchange model on the triangular lattice.These phases are ferromagnetic, 120 • Néel, dual 120 • , and two different stripe states.For four of them, we find explicit expressions for their spin-wave spectra.To the best of our knowledge, the spin-wave spectrum of the 120 • Néel phase with anisotropic terms has not been discussed previously.We demonstrate that it is, generally, non-reciprocal in this case, ε α,k = ε α,−k .
Next, we analyze the transition boundaries between these principal phases as given by the instabilities of their magnon spectra.We find that such an approach closely and reliably reproduces phase boundaries obtained by a numerical optimization of classical energy in large clusters of spins [36], offering obvious advantages over this technique.In agreement with prior studies [32,35,36], we also find consistent discrepancies of the phase boundaries provided by the Luttinger-Tisza method [52] and discuss a potential reason for that.
For the 120 • phase close to the Heisenberg limit, we find an instability toward a long-range spiral state that is similar to the Z 2 vortex state found in the triangularlattice Kitaev-Heisenberg model [53,54].With the correspondence to that model discussed below in more detail, we note that the identified transition boundary is similar to the boundary of the spin-liquid phase advocated in our previous work [39] for the quantum S = 1/2 case, suggesting a possible relation between the two.
In the present study, we also explore the ferromagnetic and stripe parts of the phase diagram in order to check whether the regions dominated by anisotropic interactions can lead to strongly frustrated and highly degenerate states.The on-site magnetization is nearly classical and quantum fluctuations are negligible for most of these regions even in the quantum S = 1/2 limit.Enigmatically, however, the ordering Néel temperature calculated from the spin-wave spectrum in one of the stripe phases is suppressed in the vicinity of a surface of parameters in the 3D parameter space.This suppression originates from the gapless pseudo-Goldstone spin-wave modes, which occur due to an accidental degeneracy, with the Mermin-Wagner theorem dictating T N = 0 for a two-dimensional (2D) system.Although quantum fluctuations do induce a gap in the pseudo-Goldstone spectrum via an order-by-disorder effect [55], the ordering temperature remains suppressed in that region compared to the mean-field expectations.Thus, while the system is almost classical, large values of the factor f = T MF /T N , which is used to identify a proximity to a quantum-disordered state [23,35], can be highly misleading, questioning it as a useful measure in such cases.
Crucially, this surface of accidental degeneracy in the anisotropic-exchange model is identified as corresponding to an extended Kitaev-Heisenberg model.This latter model possesses emergent symmetries that naturally lead to the pseudo-Goldstone modes in the quasiclassical limit, thus explaining the enigmatic trends described above.There are also additional symmetry transformations within that model, known as Klein dualities [56,57], that allow one to make deeper connections between different parts of the parameter space.
Using these insights, we have performed DMRG studies of the quantum S = 1/2 anisotropic-exchange model in previously unexplored parts of the phase diagram [39].We have validated quasiclassical phase boundaries discussed above and verified previous studies of the Kitaev-Heisenberg model [53,54,[58][59][60] that are exposed here in a wider parameter space.We demonstrate that the so-called nematic phase [53,54,59,60] corresponds to the boundary between two stripe phases and does not represent a separate state in the quantum limit.
We emphasize that the region that we have previously identified as a spin-liquid phase in Ref. [39] includes a sector of the line that corresponds to the Kitaev-Heisenberg model.This implies that previous numerical studies of the S = 1/2 Kitaev-Heisenberg model [54,60] must have overlooked the spin-liquid phase, either due to smallness of their clusters [54] or due to periodic boundary conditions [60] that are unfavorable for DMRG.
However, the most important implication of the correspondence to the Kitaev-Heisenberg model is that it necessitates an existence of another spin liquid, the one that is Klein dual to the spin liquid found in Ref. [39].In our present DMRG study, we confirm the existence of this dual spin-liquid phase.Interestingly, for the exchange matrix written in crystallographic axes, the dual spin liquid occurs in the region dominated by anisotropic terms.We use the structure factor S(q) to argue that the dual spin liquid can be seen as a result of a "melting" of the dual 120 • phase, just as the spin liquid of Ref. [39] is a molten 120 • phase, with both phases maintaining the shapes of the structure factor similar to that of their parent ordered states.The confirmation of the dual spin liquid strengthens our case for both of them.
The paper is structured as follows.Sec.II presents the model and simplified classical phase diagram.Sec.III shows spin-wave spectra of the key phases.In Sec.IV phase boundaries are discussed.Sec.V discusses finitetemperature transitions.In Sec.VI a transformation to the extended Kitaev-Heisenberg model is given.Sec.VII presents the DMRG results.We conclude in Sec.VIII, and the Appendices contain further details.

II. MODEL AND CLASSICAL PHASES A. Model
In systems with spin-orbit coupling, magnetic degrees of freedom are entangled with the orbital orientations that are tied to the lattice due to crystal fields [8].Because of that, Hamiltonians of the low-energy effective pseudo-spins involve bond-dependent interactions that obey only discrete symmetries of the underlying lattice, thus explicitly breaking spin-rotational symmetries [35].
The most general nearest-neighbor spin-orbit-induced anisotropic-exchange Hamiltonian, applicable to a variety of systems, can be written as [16] where S T i = (S x i , S y i , S z i ), ij denotes nearest-neighbor sites, and Ĵij is a 3 × 3 exchange matrix that depends on the bond orientation.Since the spin-rotational symmetries are, generally, absent, constraints on the matrix elements of Ĵij come solely from the space group symmetry of the lattice.
The effect of these constraints on the Hamiltonian (1) for the triangular-lattice materials, such as YbMgGaO 4 and others, has been thoroughly discussed in Refs.[25,31,35,38,40].Here we would like to provide a brief and intuitive derivation of the main results.
Consider the Hamiltonian (1) on the δ 1 bond, see Fig. 1, with the x axis parallel to it As can be seen from Fig. 1, the symmetries of the lattice are the C 3 rotation around the z axis, C 2 rotation around each bond, site inversion symmetry I, and two translations, T 1 and T 2 along δ 1 and δ 2 , respectively [35].These symmetries eliminate most of the elements of the exchange matrix.First, the 180 • rotation around the δ 1 bond changes y → −y and z → −z, but should leave the two-site form (2) invariant, leaving us with Then, inversion with respect to the bond center, which is a combination of the site inversion and T 1 translation, and change i ↔ j should also leave (3) invariant, allowing only the symmetric off-diagonal term, J zy = J yz .
Renaming it as J zy = J z± , and rewriting the diagonal terms using XXZ-like parametrization J zz = ∆J, with J = (J xx + J yy )/2 and J ±± = (J xx − J yy )/4 yields the two-site Hamiltonian for δ 1 in a "spin-ice" form [61] A sketch of the triangular-lattice layer of magnetic ions (empty circles) embedded in the octahedra of ligands (black dots) with the primitive vectors.Thick (blue) bonds are between magnetic ions and ion-ligand bonds are the thin solid (dashed) lines for above (below) the plane.
For the other bonds, using the C 3 invariance with the z axis to transform (3) to the δ α bond in Fig. 1, changes the Ĵ1 matrix in Eq. ( 3) to Ĵα = R−1 α Ĵ1 Rα where is the rotation matrix, or, explicitly where the abbreviations are cα = cos φα and sα = sin φα .Altogether, the most general Hamiltonian (1) on the triangular lattice becomes where c(s) α = cos(sin) φα as above, the bond angles φα are that of the primitive vectors δ α with the x axis, φα = {0, 2π/3, −2π/3}, and the spin operators are in crystallographic axes that are tied to the lattice, see Fig. 1.The Hamiltonian (7) is naturally divided in the bond-independent XXZ part and the bond-dependent anisotropic J ±± and J z± terms, also referred to as the pseudo-dipolar terms [38], which generally break continuous spin-rotational symmetries down to discrete ones.

B. Classical phase diagram
Since there are four parameters in the model (7), its parameter space is three dimensional, with the fourth degree of freedom setting the energy scale.We apply one physical constraint on it by assuming the easy-plane type of the XXZ anisotropy, 0 ≤ ∆ ≤ 1, as this is natural for a variety of layered systems of interest [18,26].As was pointed out in Ref. [35], one needs to consider only positive J z± since the global π rotation around the z axis that should leave the Hamiltonian invariant is equivalent to changing J z± → −J z± .With these two constraints, one can map the entire 3D parameter space on a cylinder, with the vertical axis represented by the XXZ anisotropy ∆, J z± as the radial, and J ±± /J as the polar variables, so that each horizontal cut represents an entire 2D phase space of the model (7) for a fixed ∆.
In Figures 2(b) and (c) we use a parametrization such that J 2 + (2J ±± ) 2 + (5J z± ) 2 = 1, with the choice of numerical coefficients made to exaggerate the region where all parameters are of the same order, J ±± , J z± J.
The XY and the Heisenberg limits of the XXZ part of the model, ∆ = 0 and ∆ = 1, correspond to the bottom and the top of the cylinder, respectively.The 2D phase diagram of the latter is shown in more detail in Fig. 2(c).By the energy minimization for the commensurate single-Q states, there are five ordered phases shown in the classical phase diagram in Figs.2(b) and (c) with spin arrangements shown in Fig. 2(c) and the ordering vectors in the Brillouin zone in Fig. 2(a).Two of the phases are favored by the XXZ part of the model (7), the ferromagnetic phase with Q = Γ and the 120 • phase with Q 120 • = K, for J < 0 and J > 0, respectively.
Although not so obvious, the two stripe phases are favored by the J ±± and J z± bond-dependent terms that are selecting the states that satisfy them fully on one of the bonds and partially on the others [39].While the stripes have the same ordering vector, Q = M or equiv-alent, they differ by the mutual orientation of spins and bonds.In the stripe-x phase, favored only by the J ±± < 0 term, spins are in plane and along one of the bonds.In the stripe-yz phase, spins are perpendicular to one of the bonds and are also tilted out of plane, taking advantage of both J ±± and J z± terms, see Fig. 2(c).
The remaining small region is the dual 120 • phase with ordering vector Q d120 • = K/2 or equivalent.While the reason for this terminology and the logic behind this state will be made clear in Sec.VI, it is related to the conventional 120 • order via the so-called Klein duality transformation [56].The dual 120 • state is a 12-sublattice state, which is a combination of four counterrotating 120 • structures, shown in different colors in Fig. 2(d).
The classical per-site energies of these phases [in units of S(S + 1)] are as follows , where we abbreviated J c = [J(1 − ∆) + 4J ±± ] /2 and the (negative) out-of-plane tilt angle of the stripe-yz state is found by energy minimization as: tan 2θ = −2J z± / J c .Our discussion reproduces results of previous studies of the single-Q ordered ground states of the model (7) that have identified the stripe and 120 • states for J > 0 [32,35,38,39]; see, in particular, Ref. [37].We also extend the same approach to the entire available parameter space.However, as was first pointed out in Ref. [36] using numerical energy optimization in large clusters, more complicated multi-Q ordered structures become ground states of the classical model near the phase boundaries of the stripe and 120 • regions.We confirm these findings by studying instabilities in the spin-wave spectra,-see Sec.IV,-and also identify a different instability within the 120 • phase for a range of ∆ near the Heisenberg limit of the XXZ term in (7).This instability is toward a different multi-Q state with a long-range spiral-like distortion of the 120 • order that is similar to the Z 2 vortex state discussed previously for the triangular-lattice Kitaev-Heisenberg (or K-J) model [53,54].Our discussion of the correspondence of the model (7) to the K-J model in Sec.VI extends this earlier finding to a broader range of parameters.
Since the continuous spin-rotational symmetries in model ( 7) are broken, one generally expects gapped spin excitations in all ordered phases.This is indeed true for the stripe-x and for most of the stripe-yz parts of the phase diagram, where in the latter part, a peculiar accidental degeneracy exists along a 2D surface of parameters that is related to duality relations discussed in Sec.VI A. However, the ferromagnetic (FM) and the 120 • states of the classical model exhibit an accidental degeneracy everywhere in their 3D regions of stability, also referred to as an emergent symmetry in Ref. [38].This degeneracy means that their spectra are gapless and the orientation of their spin configurations is not fixed within the lattice plane for ∆ < 1, or at all for ∆ = 1, offering examples of the emergent U (1) and O(3) symmetries, respectively.However, since the model ( 7) breaks rotational symmetry, it means that quantum and thermal fluctuations will select a preferred direction and gap out the spectrum.We discuss the outcomes of such a quantum order-bydisorder effect in these phases in Sec.III.
It is worth noting that most of the phase diagram in Fig. 2 is occupied by the states with quantum fluctuations that remain insignificant even for the quantum S = 1/2 limit, such as stripes and FM states.Thus, by and large, strong anisotropic terms on the triangular lattice do not seem to result in a massive degeneracy of the classical states that would indicate possible exotic phases, contrary to some early expectations [18].

III. MAGNON SPECTRA
Although the linear spin-wave spectra and transverse dynamical structure factors for the ordered single-Q spin structures can be obtained numerically, see Ref. [62] and Supplemental Material of Ref. [55], their analytical forms can be tremendously useful and informative.In this Section, we present the linear spin-wave theory (LSWT) for four out of five ordered single-Q phases shown in Fig. 2 and discussed in Sec.II above: the stripe-x, stripe-yz, 120 • Néel, and ferromagnetic states.We do not consider the spectrum of the dual 120 • state because of its complicated 12-sublattice structure.Some of our results for the stripe phases have been been discussed previously either in limiting cases [31] or with minimal details [35].We would like to emphasize that the spectrum of the 120 • phase in the presence of the bond-dependent anisotropic terms has not been calculated previously.The same is true for the ferromagnetic phase, which is, however, much simpler.
The spin-wave expansion requires a rotation of the axes on each site from the laboratory reference frame {x, y, z}, in which the Hamiltonian is typically written, to a local reference frame {x, ỹ, z} with the z along the spin's quantization axis given by the classical energy minimization for a spin configuration where S i is the spin vector in the local reference frame at the site i and Ri is the rotation matrix for that site.Thus, the spin Hamiltonian in (1) can be rewritten as where the "rotated" exchange matrix is This procedure is followed by a standard bosonization of spin operators via the Holstein-Primakoff transformation, , and a subsequent diagonalization of the harmonic part of the Hamiltonian.
In the following, we will use a two-stage rotation from the crystallographic to local reference frame, R = Rϕ • Rθ , with the first rotation in the x-y plane by an angle ϕ where θ i are the out-of-plane canting angles of the spins.
While not unique, this approach is physically intuitive and is motivated by the in-plane orderings.
A particularly simple and useful example that is used below is that of the states that are coplanar with the lattice plane, such as stripe-x and 120 • Néel states, for which all θ i = 0 and the rotation matrix (15) After this rotation, the Hamiltonian (7) that yields the LSWT becomes where ϕ i(j) are the angles of the spins with the x axis in the laboratory frame and we have omitted the terms that contribute only to the anharmonic order of the SWT.

A. Stripe-x phase
For the stripe-x phase, the ordering vector is one of the M-points and spins are split into two sublattices with spins along one of the bonds, ϕ i = ϕ 0 + (Qr i ).Choosing the ordering vector Q = M = (0, 2π/ √ 3) for convenience, dictates ϕ A = 0 and ϕ B = π for the A and B sublattices, see Fig. 2.After some straightforward algebra with the Hamiltonian in (17), separating the bonds into intra-and inter-sublattice ones, yields the 4 × 4 LSWT matrix in (19) with the 2 × 2 matrices Âk and Bk and the elements of the matrices given by In this case, the eigenvalues of the Hamiltonian matrix (19) can be found analytically by diagonalizing ĝ Ĥk instead of ĝ Ĥk , giving two magnon branches Figure 3 gives an example of the two magnon branches for a representative point within the stripe-x phase along the two paths in the Brillouin zone shown in the inset with the ordering vector at M also indicated.The left inset shows a sketch of the J z± -J ±± phase diagram for ∆ = 1.0 in Cartesian coordinates with the point indicating the chosen set of parameters.As expected, the spectrum is fully gapped due to the symmetry-breaking terms.

B. Stripe-yz phase
In the stripe-yz phase, the ordering vector and the number of sublattices are the same as in the stripe-x case, but the spins are tilted away from the laboratory frame, see Fig. 2(c).Choosing again Q = M fixes all ϕ i = π/2 while the out-of-plane angles are θ i = θ + (Qr i ), leading to θ A = θ and θ B = θ + π, where θ < 0 according to the classical energy minimization in Sec.II B. Thus, the matrix of rotation to local reference frames (15) which yields a somewhat lengthly "rotated" Hamiltonian Eq. (11), shown in Appendix A. The subsequent spinwave expansion yields the LSWT matrix (19) with the structure identical to the stripe-x case, Eq. ( 20) above.We delegate explicit expressions for A k , B k , C k , and D k for the stripe-yz phase to Appendix A for brevity.Needless to say, the magnon energies are calculated using the same expressions (22).
Figure 4 shows magnon energies for a point within the stripe-yz phase, along the same two paths in the Brillouin zone as in Fig. 3, and for the same ordering vector Q = M .Although one expects the spectrum to be fully gapped due to the symmetry-breaking terms as in the stripe-x case, there is a gapless mode at the M point, which is not the ordering vector, and the spectrum is fully gapped at M .The gap vanishes because of the choice of parameters in Fig. 4, ∆ = 1.0,J ±± = 0.2J, and J z± = 2 √ 2J ±± , that belong to a 2D surface of parameters which yields an accidental degeneracy and a pseudo-Goldstone (gapless) mode.For ∆ = 1, Eq. ( 24) gives a line of points defined by J z± = 2 √ 2J ±± and the outof-plane tilt of spins in this case is given by tan θ = −1/ √ 2. At this stage, the condition (24) can be seen as a serendipitous discovery that is sufficient to ensure an accidental degeneracy of the magnon band.For the LSWT matrix elements in (20), Eq. ( 24) follows from requiring D M = 0, which simultaneously makes A M = |B M |, see Appendix A for their explicit expressions.
While this accidental degeneracy looks similar to the case of the J 1 -J 2 model on a triangular lattice [64], the symmetry associated with it is hidden.We identify this 2D surface of accidental degeneracies in the parameter space as corresponding to an extended Kitaev-Heisenberg model, which possesses emergent symmetries that naturally lead to the pseudo-Goldstone modes in the quasiclassical limit; see Sec.VI A. For a generic choice of parameters within the 3D region of the stripe-yz phase, the accidental degeneracy is not present, but the gap remains small in a large portion of the phase diagram.In Appendix A, we provide two additional plots of magnon energies for the stripe-yz phase to substantiate this point.
We also point out that the presence of this gapless or nearly gapless mode does not affect quantum fluctuations, which remain small even in the quantum S = 1/2 limit throughout the stripe-yz phase.However, the ordering Néel temperature is necessarily suppressed in a vicinity of the 2D surface in Eq. ( 24) due to the Mermin-Wagner theorem.We discuss this dichotomy in Sec.V.

C. 120 • phase
For the 120 • phase, the ordering vector is one of the K-points, see Fig. 2(a), and spins form a three-sublattice structure.For instance, choosing the ordering vector Q = K = (4π/3, 0) and fixing the angle on the A sublattices to ϕ A , defines the other angles via ϕ i = ϕ 0 + (Qr i ).Thus, as expected for the 120 • pattern.Using these phases in the Hamiltonian (17), after some tedious but straightforward algebra, one obtains the 3 × 3 Âk and Bk LSWT matrices in ( 19) Since the nearest-neighbor interactions couple only different sublattices, the resultant Mi matrices are all traceless and are built from the hopping amplitudes, which are either bond independent, as for M1 in the XXZ term or bond dependent, as for the ones originating from the anisotropic J ±± and J z± terms where and sums over α involve only three primitive vectors δ α .With the Âk and Bk matrices (25) written out explicitly, the 6×6 LSWT Hamiltonian (19) has to be diagonalized numerically.Figure 5(a) shows the resultant magnon spectrum for a representative point in the 120 • phase for J ±± < 0 indicated in the insets along the two paths in the Brillouin zone shown in Fig. 5(b), which also shows an intensity map of the lowest magnon branch.Note that due to the three-sublattice structure, the magnetic Brillouin zone is one-third of the full one and the spectrum in Fig. 5(a) is symmetric with the middle of the Γ − K line, so the K and Γ points are equivalent.Thus the spectrum possesses only one Goldstone mode for ∆ < 1.
While one expects the spectrum to be gapless because of the emergent U (1) symmetry of the classical model discussed in Sec.II B, the unusual feature in Figs.5(a) and 5(b) is the nonreciprocal character of the magnon dispersions, ε α,k = ε α,−k [65,66], which is due to the J z± term that breaks the inversion symmetry of the LSWT Hamiltonian in the 120 • state.While Ĥ−k = Ĥk does not automatically imply the nonreciprocity, the J z± term also makes Ĥ * −k = Ĥk , which together seem to be a sufficient condition.As one can see in Fig. 5(b), the 2π/3-rotation symmetry of the spectrum is preserved, but the inversion is not.An additional figure showing the same behavior for J ±± > 0 is given in Appendix A.
As discussed in Sec.II B, anisotropic terms do not contribute to classical energy of the 120 • state, yielding an accidental U (1) degeneracy for ∆ < 1.Since the magnon spectra depend on the angles of spins with bonds, this degeneracy will be lifted by zero-point fluctuations via the order-by-disorder mechanism [55,67,68] that should select a preferred spin direction and open a gap in the spectrum.In the 120 • structure, fixing an angle of one spin fixes the rest, so the quantum energy correction is given by δE(ϕ) = −3JS/2 + α,k ε α,k (ϕ)/6, where ε α,k (ϕ) are the magnon energies that depend on the angle ϕ of a spin in one of the sublattices with the x axis, and the sum is over the full Brillouin zone.
In Fig. 5(c) we show this quantum correction with its average value subtracted, ∆E(ϕ) = δE(ϕ)− δE(ϕ) vs ϕ.Thus, for the given parameters, fluctuations pin spins to the bond directions in the manner shown in the sketch of the 120 • configuration in Fig. 5(a).That is, each spin is perpendicular to one of the bonds and bisects an angle of the triangle, corresponding to the choice ϕ = π/6 + πn/3.This choice is for the parameters to the left of the dashed line in the inset in Fig. 5(a).To the right of it, a state with spins along the bonds is chosen with ϕ = πn/3, see Appendix A. Curiously, the fact that the energy minimum must switch implies that the U (1) symmetry is retained on a 2D surface within the 120 • phase.

D. Ferromagnetic phase
For the ferromagnetic phase, the ordering vector is at Q = Γ = (0, 0) and the spin state can be described with a single-sublattice picture, i.e., all ϕ i = ϕ and θ i = θ.As is mentioned in Sec.II, for ∆ < 1 the XXZ anisotropy reduces the symmetry to U (1) and makes θ = 0.For ∆ = 1, the classical energy is insensitive to the global direction of the ordered moment and the symmetry is O(3).Having this latter case in mind, the matrix of rotation to local reference frames (15) retains its general form and we list the relevant terms of the rotated exchange matrix in Appendix A.
The spin-wave expansion yields a simple Hamiltonian with somewhat involved expressions for A k and B k where, using notations c φ = cos φ and s φ = sin φ, where c α = cos kδ α and sums over α involve only three primitive vectors δ α as before.
The magnon spectrum is given by and expressions for A k and B k (32) simplify considerably for ∆ < 1 as it forces a coplanar state with θ = 0. Figure 6 shows the magnon spectrum along the two paths in the Brillouin zone for a representative point in the FM phase as indicated in the insets.Note that J < 0. Similar to the 120 • phase, the entire 3D region of the FM phase has an accidental U (1) degeneracy of the classical model for ∆ < 1 and O(3) for ∆ = 1 as discussed in Sec.II B, hence the gapless spectrum.One can show from Eqs. ( 32) and ( 33) that the Goldstone mode should be linear in k for all ∆ < 1 and quadratic for ∆ = 1; see Appendix A for another plot of the magnon spectrum demonstrating the latter result.
Similar to the stripe states, once the direction of the ordered moment is chosen, the C 3 symmetry is broken, while inversion remains; see Fig. 6.The direction of the ordered moment in the FM state is selected via order by disorder, as for the 120 • state.The same type of analysis for the parameters in Fig. 6 pins spins in the manner shown in the sketch of the ordered structure, i.e., perpendicular to one of the bonds with ϕ = π/2 + πn/3.We find that for ∆ < 1, the preferred orientation of the ordered moment in the FM phase of the phase diagram switches from the perpendicular (ϕ = π/2) to the parallel (ϕ = 0) orientation in a region roughly near J ±± ≈ 0, although the details of such a switch can be more complicated.In analogy with the 120 • phase, this behavior also implies a surface where U (1) symmetry is retained.
For ∆ = 1, the quantum selection of the direction of the ordered moment is more complicated because of the O(3) degeneracy.However, in case of parameters that fall on the line corresponding to the K-J model (see Sec. VI), order by disorder selects the so-called cubic axes as a set of preferred directions, in agreement with Ref. [53].

IV. INSTABILITIES OF MAGNON SPECTRA
In the preceding Sections, principal single-Q ordered phases of the model ( 7) have been identified and their spectra of excitations have been found.Here we advocate a straightforward and fruitful approach that provides further insights into the phase diagram of the model via an investigation of the stability boundaries of the magnon spectra.Generally, a magnon spectrum of a state is defined in a corresponding region of the phase diagram.As a function of the model parameters, a magnon branch may soften and become imaginary, or have ε 2 k < 0 at one or at a set of k points, indicating a transition to a different state, which is referred to as magnon instability.
The well-known examples of such a behavior are magnon-softening transitions in the J 1 -J 2 model on the square and triangular lattices and in some of its extensions [64,69,70], in which magnon instability occurs exactly at the classical phase boundary.There are also other examples, such as the XXZ version of the same model on the triangular lattice [71], in which the magnon stability region extends beyond the boundary of its classical phase.Such an overlap of the magnon stability regions of the neighboring phases suggests a first-order transition between them [72].Finally, magnon instability may occur before the boundaries of its classical phase are reached, indicating an intermediate phase.

A. Stripe phases
First, we explore magnon instabilities in the stripe phases.Ref. [36] used a numerical optimization of energy in large clusters of classical spins and was first to demonstrate that the single-Q classical phase diagram of the model (7) of the type shown in Fig. 2 may be incomplete.It was shown that the more complicated multi-Q states, which can be described as modulated stripe phases with incommensurate ordering vectors, create a layer of intermediate states between the stripe and 120 • phases.
We support this finding using magnon instabilities and explore it in a wider region of phase space.Our Figs. 7  and 8 show the 2D J z± -J ±± phase diagrams for several values of ∆ with magnon instability boundaries obtained from Eq. ( 22) for the stripe-x and stripe-yz phases.While not shown, magnon instability boundaries for the choice of ∆ ≈ 0.56 are virtually indistinguishable from the phase boundaries of the multi-Q regions identified numerically in Ref. [36].
Moreover, the k-points of the observed magnon instabilities also coincide with the ones identified in Ref. [36] as the new ordering vectors.For instance, for most of the boundary region with the 120 • phase, the magnon instability occurs at the incommensurate vector Q located along the line between the Γ and K points in the Brillouin zone, in agreement with Ref. [36].Needless to say, our method of identifying multi-Q transition boundaries offers obvious advantages over the numerical approach, as it requires only the knowledge of the magnon spectrum of the single-Q states.In Fig. 7, which displays magnon instability lines (solid lines) from the stripe phases for ∆ = 1 that occur before the classical boundaries (dotted lines) are reached, we also compare them to the results of the Luttinger-Tisza (LT) method [52], shown by the dashed lines.The LT method, with the so-called "weak" constraint, essentially amounts to finding a minimum of the lowest eigenvalue of the Fourier transform of the exchange matrix (6) J αβ (q) in the q-space to find a single-Q ground state.It has been used previously to study the phase diagram of the model (7) [32,35,36] and was instrumental in identifying stripe structures as the ground states of the model.However, it was noted that the LT method sometimes fails and finds an incommensurate state even when the true ground state is a commensurate single-Q state [32,35].In fact, these problems are not specific to the model ( 7) and have been known and understood as coming from the "weak" nature of the constraint [73].
Figure 7 illustrates these consistent discrepancies of the LT method with the results of the magnon instability approach for the phase boundaries of the stripe phases.Curiously, the two methods agree for one point {J ±± , J z± } ≈ {0.1, 0.3}, which is, actually, a point that belongs to the line corresponding to the K-J model, see Sec.VI, so one can suspect an emergent U (1)/O(3) symmetry along this line as the reason for restoring the validity of the LT approach.The other range where LT agrees with our method is the range of J z± 0.5 of the stripe-to-stripe boundary.The reason for that is not totally clear.We also note that the LT method breaks down within the 120 • phase for any ∆ at any non-zero 8. Same as in Fig. 7 for several ∆; dotted lines are transitions in Fig. 2, solid lines are magnon instabilities of the stripe and 120 • phases (half-ellipse shape).Filled areas are stable single-Q phases and blank regions are multi-Q states.
J ±± , predicting incommensurate states [32].Altogether, significant care must be exercised when using it.Lastly, we point out that the energy gain created by the modulation of the stripe states in the multi-Q structures was shown to be tiny, ∼ 10 −3 J [36].This smallness may explain why no multi-Q states were detected by DMRG [39] for the quantum S = 1/2 case.

B. 120 • phase
While the border regions of the stripe-x and stripe-yz phases with the 120 • phase get replaced by the multi-Q states, magnon instabilities within the 120 • phase also indicate intriguing behavior.First, the 120 • magnons defined by the algebra in Sec.III C are stable with respect to the anisotropic J ±± term far beyond the classical boundaries of the 120 • phase regardless of the value of the XXZ anisotropy ∆, with stability boundaries shown in Fig. 8 by the half-ellipse lines.This fact strongly suggests a first-order transition from the 120 • to the stripe phases if quantum fluctuations are included.
Second, the spectrum stability with respect to the anisotropic J z± term is more drastic.In fact, while the XXZ anisotropy provides a finite range of stability to magnons in the 120 • phase, magnons in the ∆ = 1 limit are unstable to any finite value of J z± , see Fig. 8.This instability is toward a different multi-Q long-range spiral state, which corresponds to the spectrum softening at three symmetric k points in the immediate vicinity of the Γ and K points.This new state is very similar to the so-called Z 2 vortex state that was discussed intensely in the triangular-lattice Kitaev-Heisenberg model [53,54].In Sec.VI, we discuss the correspondence of our model (7) to that model along the line in the ∆ = 1 plane.Our present consideration shows that a similar state exists in a significantly wider parameter space.
We note that both trends are commensurate with our previous DMRG results for the quantum S = 1/2 case in Ref. [39].Namely, DMRG observes only a direct transition from the 120 • to the stripe phases vs J ±± and the transition is likely first order [39].The behavior of the transition boundary of the spin-liquid phase found in Ref. [39] and the overall shape of the phase space occupied by it is similar to the boundary of the longrange spiral state discussed above.The main difference is that the footprint of the spin liquid in the ∆ = 1 plane is smaller and an ordered 120 • -like phase is stabilized for J z± 0.27.Both the spiral (Z 2 vortex) and the spin-liquid phases shrink and disappear upon reducing ∆.This similarity may argue for a possible relation between the two, suggesting the spin-liquid state of Ref. [39] to be a molten Z 2 vortex state rather than a molten 120 • state.However, to add a word of caution, we have found no indication of any sizable shift of the peaks in the spinliquid structure factor from the commensurate ordering vector of the 120 • state, nor have we seen any traces of the chirality in the spin-liquid wave-function that is non-zero in the Z 2 vortex state, see Ref. [39].

V. QUANTUM AND THERMAL FLUCTUATIONS IN THE STRIPE PHASES
In this Section, we verify whether strong anisotropies in the model (7) can result in strong quantum or thermal fluctuations in the stripe phases using a quasiclassical approach.

A. Quantum fluctuations
The on-site ordered moment for T = 0 within the LSWT takes the standard form where N is the total number of lattice sites, the k-sum is over the full Brillouin zone of the lattice, µ = 1, 2 numerates magnon branches in the stripe phase, v µk are the Bogolyubov parameters of the transformation that diagonalizes the Hamiltonian in Eq. ( 19), and we have used the symmetry of the sublattices in the stripe phases.
In the stripe phases, given the symmetry between sublattices, normalization of the Bogolyubov parameters is u 2 µk − v 2 µk = 1/2 and their squares can be found analytically following Ref.[74] where M 11 [λ] is the first minor of λ Î − ĝ Ĥk with the LSWT Hamiltonian matrix (19), the product in the denominator is over three out of four eigenvalues III A, and the explicit expression for M 11 [λ] given by with the matrix elements for the stripe-x phase given in Eq. ( 21) and for the stripe-yz phase in Appendix A.
The intensity map of S obtained from (35) for S = 1/2 and ∆ = 1 is shown in Fig. 9(a) throughout the stripe phases.One can see that the on-site magnetization is nearly classical and quantum fluctuations are really negligible for these regions all the way to the transition boundaries even in the quantum S = 1/2 limit.While this observation may seem natural for the strongly gapped stripex phase [31], it is somewhat less obvious in the stripeyz phase because its spectrum has low-lying pseudo-Goldstone modes owing to an accidental degeneracy, see Sec.III B. As will be argued in Sec.VI A, this is because the model in this region is related to a ferromagnetic state, which translates into a non-divergent contribution of the gapless region to the fluctuations in (35).
We also note that we do not show the results for the ordered moment S from the 120 • and ferromagnetic states.In the former state, the spectrum is unstable toward the spiral-like state for ∆ = 1, as discussed in Sec.IV B above.In the latter state, the order-by-disorder selection of the ordered moment direction is complicated for ∆ = 1 and has to be done numerically, see Sec.III D.
These problems are avoided for smaller ∆ and we present such calculations in our Fig.10(a) for ∆ = 0.5.Here the order-by-disorder selection in both 120 • and FM phases is straightforward, see Sec.III C and D, and for the 120 • state the formalism for calculating S is the same as in Eqs.(35) and (36).Note that the empty regions in Fig. 10(a) are from the dual 120 • state and the multi-Q regions discussed in Sec.IV A, for which it is simply challenging to perform similar calculations.There is, however, hardly any doubt that they are equally well ordered as the rest of the phase diagram.

B. Thermal fluctuations
It is common in studies of frustrated magnets and their models to characterize them with the "frustration ratio" f = θ CW /T N [35,75], the ratio of the Curie-Weiss temperature θ CW to the actual ordering temperature T N , as a measure of a proximity to potentially exotic states.
A more refined alternative to θ CW in anisotropic models is the mean-field transition temperature [76] where λ min (Q) is the lowest eigenvalue of the Fourier transform of the exchange matrix (6), Ĵij , at the ordering vector.For the stripe orders, Q = M and λ min (Q) is simply the classical energy per unit cell from Eq. ( 9): λ x = 2E stripe−x for the stripe-x phase and λ yz = 2E stripe−yz for the stripe-yz phase.
The transition temperature can be found from the vanishing point of the ordered moment within the LSWT where n (ε µk ) is the magnon occupation number.However, this approach is known to overestimate T N [77] as it leads to S T ∝ (T N − T ) for T → T N , not to a power law.Instead, we use the self-consistent RPA method of Refs.[74,78], in which one introduces the T -dependence in the magnon spectrum as ε k = 2 S T ε k .Then the spin Green's function is given by [78] and Eq. ( 39) is replaced by a self-consistent condition At T → T N , the ordered moment S T → 0 and (41) yields the ordering temperature Note that this approach is only valid for S = 1/2, with expressions being more complicated for larger spins [78].
The results of the calculations of f −1 = T N /T MF are shown in Figs.9(b) and 10(b) for ∆ = 1 and ∆ = 0.5, respectively.For the stripe-x phase, the ordering temperature is close to the mean-field temperature except near the transition boundaries as expected.The results for the stripe-yz phase are more intriguing.There is clearly a line of parameters where T N is exactly zero in both figures.This line is given by Eq. (24), which reduces to J z± = 2 √ 2|J ±± | for ∆ = 1 in Fig. 9(b).It corresponds to the condition for the accidental degeneracy and for the pseudo-Goldstone mode in the magnon spectrum, see Sec.III B. On a technical level, it is easy to understand why the transition temperature vanishes.In Eq. ( 42), u k , v k ∼ const and the magnon energy ε k ∼ k 2 near the accidental degeneracy point, leading to a logarithmic divergence.Physically, this divergence is a manifestation of the Mermin-Wagner theorem that forbids ordering in 2D in the presence of a continuous symmetry.
Several comments are in order.First, both ferromagnetic and 120 • magnon spectra are gapless and the states retain a continuous symmetry on the classical level.Therefore, by the same Mermin-Wagner theorem, the ordering transition temperature T N is zero throughout their regions, as is shown in Figs.9(b) and 10(b).
Second, we have provided a clear example that a large frustration ratio can be highly misleading as a guide to a quantum-disordered region in 2D.In our case, quantum fluctuations are negligible and the ordered moment is nearly classical, but T MF /T N can be large, which naively would imply a proximity to a spin-liquid state [35].Third, because the degeneracy leading to pseudo-Goldstone modes is accidental, quantum fluctuations will induce a gap in them via an order-by-disorder effect [55].We delegate the quantitative discussion of this effect to Appendix B which needs some further developments discussed in Sec.VI A. Although the degeneracy will be lifted, the ordering temperature will remain suppressed in that region compared to the mean-field expectations.
Lastly, while a detailed phenomenology of the pseudo-Goldstone modes and suppressed ordering temperature is achieved and demonstrated, it still leaves the question on the nature of the accidental degeneracy wide open.The system is in a well-ordered stripe phase, with the inplane and out-of-plane angles of spins seemingly pinned by the energy minimization, with no obvious combinations of ϕ and θ manifesting a continuous symmetry of the spin configuration in the crystallographic axes.In order to elucidate the nature of this symmetry, we need to consider a different set of axes used in the anisotropic bond-dependent models, which we discuss next.

VI. CUBIC AXES AND GENERALIZED KITAEV-HEISENBERG MODEL
The underlying crystal structure of the triangularlattice model considered in this work is that of the 2D arrangement of the edge-sharing octahedra of ligands, such as O 2− , surrounding magnetic ions, such as Yb 3+ in the case of YbMgGaO 4 , see Fig. 11.It is also a particular example of the 2D and 3D crystal structures built from the edge-sharing octahedra, such as the honeycomb and pyrochlore lattices, see Ref. [16] for an overview.
Because of the crystal field effect and since the superexchange processes between magnetic ions are mediated by the ligands, this geometry precipitates bond-dependent anisotropic-exchange interactions between magnetic moments of the ions with strong spin-orbit coupling [12,13,16,56].At the same time, this robust structure maintains high lattice symmetry of the resultant spin models, such as the one discussed in Sec.II for our model (7), which limits the number and the type of allowed terms.
This lattice arrangement also makes natural the choice of axes that is different from the crystallographic ones that have been used in this work so far, the so-called cubic axes.The cubic axes are directly tied to the edges of the cubes, that is, to the bonds of the magnetic ions with ligands, see Fig. 11, as opposed to the crystallographic ones that consider only magnetic ions.Crucially, aside from this physical justification, some of the hidden symmetries of the model become apparent in this language.
One of the choices for the cubic axes is illustrated in Fig. 11, where we also outline the octahedron of the ligand sites and the ligand-to-magnetic-ion bonds that are forming cubic shapes, with bonds and sites of the lattice being the same as in Fig. 1.Then, the transformation from the cubic to crystallographic reference frame, S cryst = Rc S cubic , is given by Next, the model ( 7) in the cubic axes can be rewritten as the extended Kitaev-Heisenberg (K-J) model where we use conventions and notations borrowed from the much-studied extended K-J model on the honeycomb lattice [12].The diagonals of the faces of the cubes in Fig. 11 that connect magnetic ions form the triangular lattice with three different bonds denoted as {X,Y,Z} ≡ {±δ 1 , ±δ 2 , ±δ 3 }.They are perpendicular to the corresponding cubic axes, i.e., the X-bond is perpendicular to the x axis, Y to y, and Z to z, respectively.In the model (44), these bonds are numerated as ij γ with the triads of {α, β, γ} being {y, z, x} on the X bond, {z, x, y} on the Y bond, and {x, y, z} on the Z bond, see Ref. [16] for the model in terms of exchange matrices.
The parameters of the extended K-J model ( 44) are related to the parameters of the original model (7) as This relation was previously discussed in Refs.[15,16] with slightly different factor and axes conventions.
In Sec.II B, we have discussed invariance of the model (7) to the simultaneous change of sign of the J z± -term and a global π-rotation of the crystallographic axes about the z axis as a justification to consider only J z± > 0. In order to access J z± < 0 in Eq. ( 45), rotation of the crystallographic axes S x(y) → −S x(y) , leads to the following transformation in the J 0 KΓΓ language [15]    (46) which returns the same extended Kitaev-Heisenberg model (44) and is, thus, self-dual.
r 5 A e 8 1 w 8 K T Z j 4 < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 2 G 0 + + B 4 U y x x F r K B q l A B 4 P p a F h w 8 r 5 A e 8 1 w 8 K T Z j 4 < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 2 G 0 + + B 4 U y x x F r K B q l A B 4 P p a F h w 8 r 5 A e 8 1 w 8 K T Z j 4 < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 2 G 0 + + B 4 U y x x F r K B q l A B 4 P p a F h w 8 s v H y 2 5 p z J z D r 6 A e f 1 A 9 X k k 7 E = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " F / / z h 9 s v H y 2 5 p z J z D r 6 A e f 1 A 9 X k k 7 E = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " F / / z h 9 s v H y 2 5 p z J z D r 6 A e f 1 A 9 X k k 7 E = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " F / / z h 9 s v H y 2 5 p z J z D r 6 A e f 1 A 9 X k k 7 E = < / l a t e x i t > e J 0 -only < l a t e i t s h a 1 _ b a s e 6 4 = " 2 G 0 + + B 4 U y x x F r K B q l A B 4 P p a F h w 8 r 5 A e 8 1 w 8 K T Z j 4 < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 2 G 0 + + B 4 U y x x F r K B q l A B 4 P p a F h w 8 r 5 A e 8 1 w 8 K T Z j 4 < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 2 G 0 + + B 4 U y x x F r K B q l A B 4 P p a F h w 8 r 5 A e 8 1 w 8 K T Z j 4 < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 2 G 0 + + B 4 U y x x F r K B q l A B 4 P p a F h w 8  7) for ∆ = 1 as in Fig. 2. The ellipse is the line of Jz± = 2 √ 2|J±±| for which model ( 7) corresponds to Kitaev-Heisenberg model (44).Special points with higher symmetries and some Klein duality connections are highlighted, see text.

A. Degeneracies and Klein duality for ∆ = 1
As was mentioned in Secs.III B and V B, in the ∆ = 1 plane of the phase diagram in Fig. 2 there is a special line defined by J z± = 2 √ 2J ±± , see Eq. ( 24), for which accidental degeneracy of the magnon spectrum in the stripeyz phase occurs.One can immediately see from Eq. ( 45) that along this line two parameters of the model vanish, Γ = Γ = 0, reducing the model in the cubic axes (44) to a simpler and more symmetric Kitaev-Heisenberg model with J 0 = J + 2J ±± and K = −6J ±± .This line of correspondence to the Kitaev-Heisenberg model is shown in the phase diagram in Fig. 12 as an oval, with the left half of it obtained via a relation in (46) for J z± < 0.Not only does this correspondence to the highersymmetry model hint at the source of the enigmatic degeneracies, but it also provides a deeper insight into connections between different parts of the phase diagram, thanks to the prior works on the model (47) on the honeycomb and triangular lattices [53,54,56,57].
Contrary to their explicit anisotropic character, compass models often exhibit continuous symmetries in the classical limit [10,79].For instance, it is easy to see that for a classical ferromagnetic state, the Kitaev term in ( 47) is invariant under the global spin rotation, thus, demonstrating an emergent O(3) symmetry [53,54].This consideration is directly relevant to the accidental degeneracy in the stripe-yz phase via the so-called Klein duality transformation [53,54,56,57].This is a foursublattice transformation, in which one spin, S(r), is left intact while spins connected to it via the X(Y,Z) bonds, S(r + δ γ ), are rotated around the x(y, z) cubic axes by π.Crucially, this transformation leaves the K-J Hamiltonian (44) invariant, with the parameters redefined as and in terms of J and J ±± as One can easily see in Fig. 11 that the described spin transformation converts the stripe-yz state into the ferromagnetic one.Since the FM state of the Kitaev-Heisenberg model ( 44) is invariant under global spin rotation, the stripe-yz state must also be invariant under a corresponding four-sublattice spin rotation, demonstrating accidental O(3) symmetry that naturally leads to the pseudo-Goldstone modes in the quasiclassical limit.
In the crystallographic notations, the out-of-plane spin angle in the stripe-yz phase along the K-J line is tan θ = −1/ √ 2, see Sec.III B, which is precisely along one of the cubic axes, see Fig. 11.Interestingly, as was shown for the K-J model on the honeycomb [80,81] and triangular lattices [54], precisely this orientation is chosen by quantum fluctuations from a classically degenerate O(3) manifold of states.In our case, this choice is made on the classical level by restricting ourselves with the two-sublattice stripe state of the surrounding phase in Fig. 12.
We underscore that the Klein duality is not restricted to the classical limit as the preceding discussion may seem to suggest, but it rather implies an exact similarity of the ground state and excitations, up to the described spin rotations [56,57].This means that the entire section of the K-J line in Fig. 12 belonging to the stripe-yz phase is dual to the entire section of the same line within the FM phase, as can be confirmed from Eq. ( 49), with the point marked as "K-only" being self-dual.For instance, one of the implications of the Klein duality is the existence of a point in the stripe-yz phase that is dual to the isotropic Heisenberg ferromagnetic point and thus must be free from any quantum fluctuations [14].Its coordinates, according to Eq. ( 49), are J = J ±± > 0 and J z± = 2 √ 2|J ±± |, and it is marked as the " J 0 -only" point in Fig. 12 with the dashed line emphasizing the duality relation.
In practice, Klein duality means that the calculations for quantities such as the Néel temperature in stripe-yz phase along the K-J line can be performed in the ferromagnetic phase instead.Since the K-J model is O(3) symmetric only on the classical level, of interest is the gap in the pseudo-Goldstone mode that is generated by the order-by-disorder effect.In Appendix B, we present calculations of the Hartree-Fock corrections to the magnon spectrum for the ferromagnetic Kitaev-Heisenberg model and use Klein duality to obtain the self-consistent transition temperature in the stripe-yz phase.Although the fluctuation-induced gap E g in the magnon spectrum is small, the transition temperature is very sensitive to it, T N ∝ T MF / ln (T MF /E g ) [77], leading to a finite T N except for the J 0 -only point.The resultant transition temperatures are still significantly suppressed compared to the mean-field expectations (38).
The implications for the rest of the phase diagram in Fig. 12 are the following.A section of the K-J line was previously identified as a "nematic" phase [53,54,59,60].In Fig. 12, this section coincides with the boundary between the stripe-x and stripe-yz phases and is unlikely to represent a separate phase on its own.
Another insight is provided by the Klein duality into the dual 120 • phase.The isotropic antiferromagnetic Heisenberg model has a Klein-dual point within the dual 120 • phase, marked as the second " J 0 -only" point in Fig. 12 with the coordinates J = J ±± < 0, see Eq. ( 49).This relation also elucidates the nature of the dual 120 • state.It is obtained from the three-sublattice 120 • state by the Klein rotations of spins around cubic axes resulting in the 12-sublattice state.Since the cubic axes are not collinear with the crystallographic ones, the resultant state is, generally, noncoplanar, see Ref. [56] for the projection of such a state onto the triangular-lattice plane.It is only when the plane of the initial 120 • state is chosen to be parallel to one of the principal planes of the cubic coordinate system, the Klein rotations will keep spins of the dual 120 • state in the same plane.This situation is sketched in Fig. 2.
The most important implication of the correspondence of the model (7) to the Kitaev-Heisenberg model (47) along the line J z± = 2 √ 2J ±± in Fig. 12 is that it necessitates the existence of a spin-liquid region in the S = 1/2 case that is Klein dual to the spin liquid found by us in Ref. [39].The confirmation of this is the subject of Sec.VII.

B. Gapless modes at ∆ < 1
As was discussed in Secs.III B and V B, the pseudo-Goldstone modes in the magnon spectrum of the stripeyz phase persist for ∆ < 1 along the lines defined by J z± = [4J ±± +J(1 − ∆)] / √ 2; see Eq. ( 24).For the extended K-J model (44) with parameters in Eq. ( 45), this condition means that only one of the off-diagonal terms remains zero along these lines (Γ = 0), while the other one does not (Γ = 0), leading to the K-J-Γ model with J 0 = J +2J ±± and ∆-dependent K and Γ Thus, while for ∆ = 1 the accidental degeneracies are associated with the high symmetry of the Kitaev-Heisenberg model, in this case their origin is more subtle.
For the K-J-Γ model, Klein duality transformation does not leave the Hamiltonian form invariant.While Heisenberg and Kitaev terms do preserve their structure with the change of the symmetric Γ -term becomes antisymmetric.For example, on the X-bond However, for ∆ < 1 and along the the K-J-Γ lines the spin orientation in the stripe-yz phase remains along one of the cubic axes.Therefore, the Γ -term does not contribute to the classical energy and leaves the accidental degeneracy of the Kitaev-Heisenberg model intact.

VII. QUANTUM REALM: DUAL SPIN LIQUIDS IN THE S = 1/2 MODEL
There are several new aspects of the anisotropicexchange model ( 7) that have been discussed in this work so far.These are the classification of all its single-Q classical phases, the identification of the instabilities of some of them to more complex multi-Q states, and various quantum and thermal effects in the magnetically ordered phases, see Secs.II-V.There is also a fruitful connection with an extended Kitaev-Heisenberg model ( 44) that uncovers hidden symmetries and relates different parts of the phase diagram to each other, see Sec.VI.Here, we build upon these insights using DMRG for the S = 1/2 model.
In our prior work, Ref. [39], we have discovered a spinliquid (SL) region of the 3D phase diagram of the model (7) using DMRG, with the sketch of its base shown in Fig. 13 as a triangle enveloping the tricritical point of the 120 • and the stripe phases.According to Ref. [39], the SL phase occupies a distorted cone shape with the base at ∆ = 1.0, the widest dimensions J z± [0.27, 0.45]J and J ±± [−0.17, 0.1]J, and the tip of the cone protruding along the XXZ axis down to ∆ 0.7.
It is important to note that for ∆ = 1.0, the SL area in Fig. 13 includes a segment of the line that corresponds to the pure Kitaev-Heisenberg model (47), discussed in Sec.VI A. Thus, in a retrospect, our discovery is also a discovery of a spin liquid in the K-J model on the triangular lattice, without the benefit of having an exact Kitaev-like solution in this geometry.We note that this finding is in disagreement with the previous numerical studies of this model [54,60], which, we believe, have missed the SL phase due to finite-size effects [54] or unfavorable boundary conditions [60].
Crucially, the Klein duality along the line of a correspondence to the Kitaev-Heisenberg model necessitates e J -only < l a t e x i t h a 1 _ b a s e 4 = " 2 G 0 + + B 4 U y x x F r K B q l A B 4 P p a F h w 8 r 5 A e 8 1 w 8 K T Z j 4 < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 2 G 0 + + B 4 U y x x F r K B q l A B 4 P p a F h w 8 r 5 A e 8 1 w 8 K T Z j 4 < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 2 G 0 + + B 4 U y x x F r K B q l A B 4 P p a F h w 8 R H g a y S N b J B f F I j e + S A 1 E m D C H J L 7 s k j e f L u v A f v 2 X v 5 b B 3 z R j M r 5 A e 8 1 w 8 K T Z j 4 < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 2 G 0 + + B 4 U y x x F r K B q l A B 4 P p a F h w 8 R H g a y S N b J B f F I j e + S A 1 E m D C H J L 7 s k j e f L u v A f v 2 X v 5 b B 3 z R j M r 5 A e 8 1 w 8 K T Z j 4 < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 2 G 0 + + B 4 U y x x F r K B q l A B 4 P p a F h w 8  12 with highlighted spin-liquid regions.The lower SL region is from Ref. [39].The magnified region emphasizes the vicinity of the dual spin-liquid region that is studied here by DMRG.
the existence of another spin liquid region, sketched in Fig. 13.Thus, in our present DMRG study, we investigate the previously unexplored parts of the phase diagram in order to confirm the existence of this dual SL phase.As in Ref. [39], we use several complementary approaches: the long-cylinder 1D "scans" with one parameter varied along the length of the 6 × 36 cylinder to explore different phases [31,46], the shorter cylinders with fixed parameters [6 × 20, "non-scans"], as well as the intensity maps of the structure factor, S(q), and correlation lengths; see also Appendix C. We use different boundary conditions and ranges of the varied parameter to exclude unwanted proximity effects.
For the DMRG calculations in the 6 × 36 cylinders, we typically perform 20 − 24 sweeps and keep up to m = 1600 − 2000 states depending on the complexity of the Hamiltonian with truncation error < 10 −5 .For the 6 × 20 cylinders, the protocol is 24 sweeps and up to m = 2000 states with truncation error < 10 −6 .In the real-space images of cylinders below, the size of the arrows represents the projection of local spin in the lattice plane and the color is used to indicate the sign of the out-of-plane tilt angle of the spins.The thickness of the nearest-neighbor bonds is proportional to the magnitude of the S i S j correlation, with the (solid) dashed lines representing (anti)ferromagnetic sign of it.In Fig. 14(a), the scan covers the entire stretch from the dual 120 • point, that is, the point that is Klein dual to the isotropic Heisenberg antiferromagnetic 120 • point, to the self-dual K-only point, K > 0, see also Fig. 12.In the notations of the Kitaev-Heisenberg model (47), this scan is from {J 0 , K} = {−1, 2}/ √ 5 to {J 0 , K} = {0, 1}, with the normalization of J 2 0 + K 2 = 1.Using duality transformation (48), one can see that the dual 120 • point indeed corresponds to the J 0 -only model.In the notations of the original anisotropic-exchange model (7), the scan is performed along the J z± = 2 √ 2|J ±± | line with ∆ = 1.0 by varying J/|J ±± | from −1 to 2, as is indicated in Figs.14(a) and 14(c).

A. Dual spin-liquid region
We note that no special boundary conditions are applied in any of the long-cylinder scans, yet the dual 120 • order appears naturally on the left end of the cylinder in Fig. 14(a) and is clearly discernible.Nearly half of this scan is occupied by what is clearly identifiable as a stripe state, which is intermediate between the stripe-x and stripe-yz states in that the spins are neither parallel nor perpendicular to any of the bonds.It is also tilted out of the lattice plane similar to the stripe-yz state.As was pointed out above, see Sec.VI A, this section of the K-J line, which is exposed here in a wider parameter space, corresponds to the boundary between two stripe phases and should not represent a separate phase in the quantum limit.In Refs.[54,60], this boundary is referred to as the nematic phase even in the quantum limit, while Ref. [82] called it a spin liquid.We believe that both terms must have been used in error, see also Appendix C for a DMRG scan across the right end of the scan in Fig. 14(a), the K-only point.
Despite a significant range of the scan in Fig. 14(a), there is a clear suppression of order between the ordered phases.The shaded regions in Figs.14(a) and (c) indicate the range that is obtained from the Klein duality of the extent of the "original" spin liquid along the K-J line in Fig. 13, J/|J ±± | ≈ [−0.22,0.23].It agrees with the minimum in the ordered moment and we refer to it as the dual SL region.Similar to the SL phase of Ref. [39], which comes from melting of the 120 • phase, the dual spin liquid appears to be born out of the dual 120 • phase.
The second scan through the dual SL region is shown in Fig. 14(b), with the vertical line within the shaded region marking its intersect with the first scan.It is also a tricritical point of the classical dual 120 • and the stripe phases.The scan in Fig. 14(b) is not on the K-J line, so the natural notations here are of the anisotropic-exchange model (7), in which this scan corresponds to the J = 0 line, the line where the XXZ part of the model is absent and the only two variables are J z± and J ±± .Keeping J z± = 2 √ 2 in order to  where magnetic order is suppressed.In Fig. 14(d), the dashed lines show the ordered moment as given by the 1/S calculations in the stripe phases, indicating that magnon instability boundaries are not unlike the boundaries of the putative dual SL phase.As is in the previous studies [39], the long-cylinder scans are only a part of the evidence for the spin liquid.In order to confirm the SL, we perform DMRG non-scans on 6 × 20 clusters with fixed parameters. Figure 15(a) shows one of them for the point that is at the core of the suggested dual SL region, at J = 0 and J z± = −2 √ 2J ±± , a tricritical point of the classical phases.Since it is on the K-J line, its coordinates in the K-J language correspond to {J 0 , K} = {−1, 3}/ √ 10.The cluster presented in Fig. 15(a) shows no discernible traces of any magnetic order and has weak nearestneighbor S i S j correlations of ferromagnetic sign.The boundary conditions are open with one site removed on each side of the cylinder to avoid spinon localization, common to the Z 2 spin-liquid states [39,46] and indicative of them.Without the sites removed, a weak pattern similar to the dual 120 • structure appears at the boundaries (not shown) and decays exponentially toward the center of the cluster, all supportive of the SL state.This behavior is in clear contrast with the results of a similar analysis of the dual 120 • ( J 0 -only) point offered in Appendix C. There, the long-cylinder scans also suggest a possible SL state, but the 6 × 20 non-scan cluster verification demonstrates a robust 12-sublattice order with a power-law decay from the boundary, as expected.
A more comprehensive insight into the spin-spin correlations and into the nature of the SL state is given by the static structure factor observable in experiments, We obtain it from the Fourier transform of the real-space spin-spin correlation function S α i S β j , determined from the DMRG ground-state wave function with all |i − j| distances that are available in the cylinder.
In Fig. 15(b), we present an intensity map of S(q) from the cluster in Fig. 15(a).It shows clear maxima at the Q d120 • = (2π/3, 0) and equivalent points, associated with the dual 120 • order, see Sec.II B. This is in accord with the "original" spin liquid of Ref. [39] having broad maxima in S(q) at the Q 120 • .These results also suggest that the dual spin liquid can be seen as a result of a "melting" of the dual 120 • phase in the same way that the spin liquid of Ref. [39] can be seen as a molten 120 • phase, both maintaining the shape of the structure factor similar to that of the parent ordered states.
Altogether, we have confirmed the existence of the second SL region in the phase diagram of the anisotropicexchange model, with properties that are in accord with the expectations based on the Klein duality relation along the K-J line.This confirmation of the dual spin liquid strengthens our case for both SL regions.
Two additional notes are in order.In Sec.IV B, we remark on a similarity of the spin-liquid region in the quantum S = 1/2 model from Ref. [39] with the quasiclassical region of instability of the 120 • phase to a multi-Q long-range spiral, or a Z 2 -vortex-like state, suggesting a possible connection between the two.This similarity may imply that the discussed spin liquids should originate from the melting of the Z 2 -vortex and dual Z 2vortex states rather than their more simple 120 • and dual 120 • counterparts.However, we reiterate here that while the long-range distortions of the ordered states found in DMRG clusters are possible, there are no traces of the expected shifts of the peaks from the commensurate ordering vectors in the SL structure factors.There is also no detectable residual chirality in the spin-liquid wave functions that may be expected to survive from the Z 2vortex states, see Ref. [39].Thus, the relation between the quasiclassical and quantum ground states of the problem in this region deserves further investigation.
Lastly, a superficial observation is that both the "original" and the dual spin liquids seem to be centered at the tricritical points of the classical single-Q ordered phases.In the K-J notations they correspond to the K = J 0 > 0 and K = −3J 0 > 0 points, respectively.Interestingly, in the anisotropic-exchange notations, the dual tricritical point corresponds to the model (7) with no XXZ terms and only anisotropic J z± and J ±± terms.While we are unable to draw any useful insight from this observation or identify any quasi-classical degeneracy that can be affiliated with this model, the chance that this point can have a special solution may be an intriguing possibility.Although the full parameter space of the nearestneighbor anisotropic-exchange model (7), or, equivalently, of the extended K-J-Γ-Γ model (44) on the triangular lattice is three dimensional, their relation to a simpler and more symmetric K-J model (47) along a one-dimensional line has proven to be very informative.Thus, we conclude this Section by summarizing our DMRG results for the quantum S = 1/2 K-J model (47) in the form of its 1D phase diagram, shown in Fig. 16, where we use the standard parametrization, J 0 = cos ϕ and K = sin ϕ, and the positive (counterclockwise) ϕ direction corresponds to the negative (clockwise) direction along the ellipse in our 2D phase diagram in Fig. 12, see Eq. (47).
We fully agree with the previous studies of the triangular-lattice K-J model [53,54,[58][59][60]82] on the location and boundaries of the FM and stripe-yz phases.Note that we deviate in the notations for what we refer to as the "120 • " and dual "120 • " phases.They were quasiclassically identified as the multi-Q distortions of such 120 • orders and were called the Z 2 vortex crystal and its dual phases [53,54].These regions are still designated as the Z 2 vortex and the dual Z 2 vortex states in the S = 1/2 case [54,60], although the evidence for the persistence of such distortions in this limit is slim.We note that some deviations from the precise spiral orders near the isotropic limit of the model may have been detected in our prior DMRG work, Ref. [39].While we do not conclusively confirm or rule them out, it seems that these deviations from the ideal orders, even if exist, play only a minor role in the energetics of the phases in the quantum limit, see Sec.IV.Hence, we simply use quotation marks for the "120 • " in referring to them.
We disagree with the previous works on the presence of a nematic phase in the quantum limit [54,60].While it may exist as a fully classical curiosity, quantum fluctuations select a stripe state that is intermediate between the stripe-x and stripe-yz phases, see Sec.VI A and Appendix C. Denoting this sector a spin liquid [82] must have been in error.In Fig. 16, we refer to it as a "stripe" phase to distinguish from the stripe-yz phase.
Our most important findings are the two regions of the spin-liquid phase.They occur in a proximity of, and, arguably, as a result of a melting of their respective parent "120 • " phases.This is also evidenced by the structure factor discussed in Sec.VII A and in Ref. [39].The connection of our "SL" region to a SL phase of the fully isotropic J 1 -J 2 Heisenberg model was established in Ref. [39], and was referred to as a SL isomorphism.We have found that the spin-spin correlations are very similar in them and that there is a path in the 4D phase diagram that provides a continuous link between the two.This connection further strengthens the argument for our "molten 120 • " phase scenario.
The respective SL regions in Fig. 16 are as follows.For the region marked "SL", K/J 0 ≈ [0.71, 1.40] and ϕ/π ≈ [0.20, 0.30].While this may or may not be a coincidence, the K = J 0 point is included in this range.By the Kleinduality transformation, the region referred to as "dual SL" occurs.Its boundaries agree well with the DMRG presented in Sec.VII A and are K/J 0 = [−3.4,−2.71] and ϕ/π ≈ [0.59, 0.61].As was discussed above, it is centered around K = −3J 0 point, which translates to the vanishing point for the XXZ part of anisotropic-exchange model (7) with only anisotropic J ±± and J z± terms present.
Lastly, the mean-field Schwinger-boson study [83] has discussed various spin liquids in the triangular-lattice K-J model.The SL region in Fig. 16 is affiliated with their "SL2" state, but that latter state has the structure factor that is very different from S(q) found by DMRG [39] and thus can be ruled out.We also reiterate that the spin liquids that we identify are not consistent with the "open spinon Fermi-surface" SL state proposed for YMGO [20].

VIII. SUMMARY
In this work, we have provided an extensive, if not exhaustive, overview of the phase diagram of the nearestneighbor triangular-lattice anisotropic-exchange model, which is relevant to a growing family of the rare-earthbased magnets and other materials with strong spin-orbit interactions.We have explored ordered phases, identified the principal ones that occupy a majority of the parameter space, and obtained explicit expressions for their non-trivial spin-wave excitation spectra.We have also identified and characterized transitions of some of the or-dered phases to more complex multi-Q states and demonstrated the effectiveness of the analysis of such transitions with the help of magnon instabilities.
In the studies of the quantum and finite-temperature effects in the well-ordered phases, a number of accidental degeneracies that lead to emergent continuous symmetries of the classical states have been discussed, and the effect of order-by-disorder on them has been analyzed.Another systematic and enigmatic accidental degeneracy has been found in the nearly classical stripe phase, leading to a strongly suppressed ordering temperature that can be falsely attributed to a proximity to an exotic state.This degeneracy has been connected to the correspondence of the original model to an extended Kitaev-Heisenberg model, in which the degeneracies have a more natural explanation.
This connection has been particularly fruitful for uncovering hidden symmetries and in relating different parts of the phase diagram to each other via the Kleinduality transformation.In a rather spectacular manifestation of the correspondence to the Kitaev-Heisenberg model, a new region of the spin-liquid state that is Klein dual to the spin liquid found by us in a prior work has been confirmed for the S = 1/2 model using an unbiased DMRG approach.Both the original and the dual spin liquids occur in proximity to their parent ordered phases that are also dual to each other.This finding strengthens our case for both spin-liquid regions in the phase diagram of the quantum model.As a corollary, we have also provided a one-dimensional phase diagram of the quantum Kitaev-Heisenberg model on the triangular lattice that updates and corrects previous results.
In conclusion, the present work, together with our prior works in the phase diagram of the anisotropic-exchange model on a triangular lattice, creates a foundation for the studies of the large group of materials with anisotropic exchanges, clears the path to a consistent interpretation of the current and future experiments, and gives important new insights into fundamental properties of quantum magnets with spin-orbit-generated low-energy Hamiltonians.observed static structure factor by a semiclassical simulation of the spin-spin correlations [23,24].Fig. 17 ε 1,2k from Eq. ( 22) vs k along the ΓM KΓ path and Fig. 18 along the ΓM K Γ path, respectively, where M is the ordering vector.The two models are away from the accidental degeneracy surface (24) discussed in Sec.III B, but still demonstrate a nearly gapless mode at the M point.This yields an ordering temperature that is suppressed compared to the mean-field value, which naively would suggest a proximity to a spin-liquid state.However, we have confirmed that no signs of strong quantum fluctuations are present in either of the cases.For S = 1/2, LSWT for the J z± -only model gives S = 0.4869, supported by the DMRG result S = 0.4795.Model B is also nearly classical with LSWT S = 0.4763 and DMRG giving S = 0.4694.120 • phase.-InFig. 19, we provide another exam- In Fig. 20, we demonstrate the magnon spectrum in the FM phase that is discussed in Sec.III D.Here it is shown for the ∆ = 1 plane of the phase diagram and also for the parameters that fall on the line corresponding to the K-J model, see Sec.VI, for which order-by-disorder fluctuations select the cubic axes as a preferred direction for the ordered moment.The Goldstone mode is clearly quadratic in k, while the rest of the features discussed in Sec.III D are preserved.Appendix B: Gap and TN on the K-J line Since the degeneracy leading to the pseudo-Goldstone modes along the K-J line in Fig. 12 is accidental, quantum fluctuations induce a gap and make the transition temperature T N (42) finite in 2D.Generally, calculations of such order-by-disorder gaps can be rather involved, see, e.g., Ref. [55].In the considered case of the stripeyz state in the K-J model, the problem is simplified by the absence of the cubic anharmonicities, because the spin orientation is along one of the cubic axes.The other simplification is the Klein duality to the ferromagnetic state, for which calculations are straightforward.
The LSWT spectrum in the ferromagnetic phase of the Kitaev-Heisenberg model (47) c α = cos kδ α , and γ k is from Eq. ( 33), see Ref. [58].The Hartree-Fock corrections from the quartic terms in the SWT Hamiltonian to A k and B k are with n = a † i a i , m α = a † i a i+δα , and ∆ α = a i a i+δα , where m 2 = m 3 , ∆ 2 = −∆ 3 , and ∆ 1 = 0. Ordering temperature for the FM phase of the K-J model calculated using Eq.(B7), lower curve, and the meanfield result, Eq. ( 38), upper curve.TN and K are in units of J 2 0 + K 2 1/2 and J0 < 0.
Within the 1/S-approximation, the spectrum becomes However, the most important effect of the quantum fluctuations is the gap at k = 0, which we have also verified by the recently developed method of Ref. [55].
Since the LSWT spectrum is parabolic, one can approximate the renormalized spectrum as which should suffice for the regularization of the divergences leading to the vanishing ordering temperature due to Bose factors in Eq. (41).Thus, our self-consistent transition temperature in the Kitaev-Heisenberg model is given by a simplified version of Eq. ( 42) with the renormalized spectrum (B6) The results of using (B7) for the FM phase of the K-J model are presented in Fig. 21, where we denote the transition temperature as T N having in mind Klein duality to the stripe-yz part of the K-J line in Fig. 12.Both T N and K are in units of J 2 0 + K 2 1/2 and J 0 < 0. Although the gap E g is small, the transition temperature depends singularly on it, T N ∝ T MF / ln (T MF /E g ) [77], leading to a finite and sizable T N with an exception to the phase boundaries and to the isotropic, O(3) symmetric K = 0 point ("J 0 -only" point), where T N remains zero.The resultant transition temperatures are still significantly suppressed compared to the mean-field results of Eq. ( 38).In Fig. 22(a), a scan is performed by varying J ±± /J from −0.7 to −0.3 for fixed ∆ = 1.0 and J z± = √ 2J, so that at the K-J line J z± = −2 √ 2J ±± as elsewhere.The scan shows a clear crossover from the stripe-x to stripeyz phase.The spins on the left edge of the cylinder in the stripe-x phase are slightly tilted off the lattice plane and continuously deform into the stripe-yz order on the right edge with only a small variation of the ordered moment near S ≈ 0.4, see Fig. 22(c).There is no indication of a nematic [54,60] or a spin-liquid state [82] at the K-J line, nor there is a sign of a proximity to any.
The dashed lines in Fig. 22(c) show the ordered moment as given by the 1/S calculations in the neighboring stripe phases, indicating a direct transition.We have also performed a DMRG non-scan 6 × 20 cluster calculation at the K-only point (not shown) with the results very similar to the central part of the scan in Fig. 22(a), showing a robust order in the form of a stripe phase with the same ordering vector at the M -point and spins oriented in between the stripe-x and stripe-yz orders.Given the smoothness of the crossover in Fig. 22(a) and that the ordering vector does not change, it is not clear to us whether the two stripe phases remain distinct in the presence of quantum fluctuations.
In Fig. 22(b), a scan is performed across the dual 120 • point by varying J ±± /|J| from −3.0 to 0.0 for fixed ∆ = 1.0 and J z± = 2 √ 2|J| (J < 0 here).There are wellordered stripe-x and stripe-yz phases in the two ends of the cylinder, separated by a region where the order is suppressed, see also Fig. 22(d).This may seem surprising as one expects the dual 120 • ordered phase in this region.One of the obvious reasons for a suppression of the order is a very large gradient of parameters in this scan, combined with a large unit cell of the 12-sublattice structure.In Fig. 22(d), the dashed lines show the ordered moment as given by the 1/S calculations in the stripe phases, showing that for S = 1/2 the dual 120 • state expands and claims some of the stripe regions, similar to the expansion of the "original" 120 • state, see Ref. [39].
To verify that the intermediate phase is not a spin liquid, we have performed a DMRG non-scan calculation in the 6 × 20 cluster for the J 0 -only, dual 120 • point (∆ = 1.0,J < 0, J ±± = J, and J z± = 2 √ 2|J|) with the results shown in Fig. 23(a).Here, the dual 120 • boundary conditions are applied at the left edge only, with a robust 12-sublattice order parameter showing a power-law decay with a finite asymptote, characteristic of the well-ordered phase, see Fig. 23(c).
The observed ordered moment of the dual 120 • state in this cluster is about S ≈ 0.09.Since the Klein duality implies that all observables between dual points should be the same, this value may seem to contradict the value of the ordered moment at the Heisenberg 120 • point, S ≈ 0.2 [84].As we verify in Fig. 23(b), the smaller value of the ordered magnetic moment is due to the mixed boundary conditions and the resultant aspect ratio of the cluster being far from being "optimal," according to Ref. [84].In Figs.23(b) and (c) we demonstrate that the 120 • order on the same cluster behaves very similarly and is in quantitative agreement with the Klein-duality expectations.

5 JFIG. 2 .
FIG. 2. (a) Brillouin zone of the triangular lattice with the ordering vectors for each phase.(b) The 3D classical phase diagram of the model (7) for the single-Q states, see text for their description.The vertical axis is 0 ≤ ∆ ≤ 1. (c) The 2D cut of (b) at ∆ = 1 with a sketch of spin structures and parametrization of the radial and angular coordinates.(d) A detailed sketch of the dual 120 • state.It is, generally, noncoplanar, see Sec.VI, and consists of 12 sublattices with the unit cell indicated.

5 J
FIG. 5. (a) Magnon energies for a point in the 120 • phase along two reciprocal k-contours (solid and dashed) in (b), with the ordering vector indicated.(b) Intensity plot of the lower branch.(c) The ϕ-dependent part of the zero-point energy vs ϕ that selects the structure sketched in (a).The spectrum retains C3, but not I symmetry, see text.

5 FIG. 6 .
FIG.6.Magnon energies along two contours shown in the inset for a representative point in the FM phase: J < 0, ∆ = 0.5, J±± = 0.4J, and Jz± = −0.1J.Inset shows the intensity plot of the magnon dispersion.A sketch depicts the orientation selected by the order-by-disorder mechanism.

FIG. 7 .
FIG.7.The Jz±-J±± phase diagram of the model(7) for ∆ = 1 with transition boundaries from Fig.2(dotted lines), from magnon instabilities in the stripe-x and stripe-yz phases (solid lines), and by the LT approach (dashed lines), see text.

FIG. 11 .
FIG.11.Same as Fig.1, side view.Octahedron of ligands is highlighted, cubic axes and bonds of the triangular lattice are shown.Spin ordering corresponds to the stripe-yz phase for the parameters along the gapless K-J line, see text.
e b a B N t I 0 8 V E F H 6 A S d o R p i 6 A b d o 0 f 0 5 N w 5 D 8 6 z 8 / L Z m n M m M + v o B 5 z X D 9 L a k 6 8 = < / l a t e x i t > K-only, > 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " F / / z h 9 FIG.12.Phase diagram of the model(7) for ∆ = 1 as in Fig.2.The ellipse is the line of Jz± = 2 √ 2|J±±| for which model (7) corresponds to Kitaev-Heisenberg model(44).Special points with higher symmetries and some Klein duality connections are highlighted, see text.
FIG.13.Phase diagram for ∆ = 1.0 from Fig.12with highlighted spin-liquid regions.The lower SL region is from Ref.[39].The magnified region emphasizes the vicinity of the dual spin-liquid region that is studied here by DMRG.

Figure 14
Figure 14 summarizes the results from the two longcylinder scans through the putative dual SL region, first along the K-J line and second normal to it at the tricritical point of the dual 120 • and the stripe phases.The portions of the phase diagram shown on the left of the cylinder images in Figs.14(a), (b) indicate the direction and extent of each scan, and Figs.14(c), (d) show the magnitude of the ordered moment S along the scans.

∆FIG. 14 .
FIG. 14. Long-cylinder DMRG scans (a) along the K-J line from the dual 120 • to the K-only point, and (b) normal to it at the classical tricritical point of the dual 120 • and the stripe phases, J = 0 and Jz± = −2 √ 2J±±.Sketches of the phases indicate the direction and extent of each scan, with (c) and (d) showing S along the scans.Shaded areas indicate putative dual spin-liquid regions.Dashed lines in (d) are the LSWT results for S , Eq. (35), in the stripe phases.

FIG. 16 .
FIG.16.The 1D phase diagram of the quantum S = 1/2 K-J model (47) using DMRG results from Sec. VII A and Ref.[39].Pure Heisenberg, pure Kitaev and their Klein-dual points are marked and the duality relation between them and SL and "120 • " phases is emphasized by dashed lines.

FIG. 17 .FIG. 18 .
FIG. 17. Magnon energies ε 1,2k from Eq. (22) vs k along the ΓM KΓ path for the model(7) with only Jz± term present (dashed lines) and for the Model B, J > 0, ∆ = 0.76, J±± = 0.26J, and Jz± = 0.45J (solid lines).Both are within the stripe-yz phase, and the ordering vector is at M as indicated.Left inset: the 2D phase diagram with the points indicating the chosen set of parameters.

5 J
FIG. 19.Same as Fig. 5. (a) Magnon energies along the two reciprocal k-contours (solid and dashed) in (b) for the parameters shown in the graph.(b) Intensity plot of the lower branch.(c) The ϕ-dependent part of the zero-point energy selecting the structure sketched in (a) with ϕ0 = πn/3.

FIG. 20 .
FIG.20.Magnon energies along the two contours shown in the upper inset for a point in the FM phase (right inset) that belongs to the line corresponding to the K-J model: J < 0, ∆ = 1.0,J±± = −0.1J,and Jz± = 2 √ 2J±±, see text.Upper inset shows the intensity plot of the magnon dispersion.A sketch depicts the orientation selected by the orderby-disorder mechanism with the chosen angles ϕ = π/2 and θ = arctan 1/ √ 2 .
FIG.21.Ordering temperature for the FM phase of the K-J model calculated using Eq.(B7), lower curve, and the meanfield result, Eq. (38), upper curve.TN and K are in units of J 2 0 + K 2 1/2 and J0 < 0.

FIG. 22 .
FIG. 22. Long-cylinder DMRG scans (a) across the K-only point (∆ = 1.0,J±± = −0.5J,and Jz± = 2J) vs J±±/J from −0.7 to −0.3, and (b) across the dual 120 • point (∆ = 1.0,J < 0, J±± = J, and Jz± = 2 √ 2|J|) vs J±±/|J| from −3.0 to 0. Both scans go from the stripe-x to the stripe-yz phase and are normal to the K-J line in Fig. 12, with the vertical lines in the clusters showing the intersect with it.Sketches of the phases indicate the direction and extent of each scan, with (c) and (d) showing S along the scans.Dashed lines in (c) and (d) are the LSWT results for S , Eq. (35), in the stripe phases.

FIG. 23 .
FIG. 23.DMRG non-scan calculation in the 6 × 20 cluster for the (a) dual 120 • point, and (b) 120 • point for comparison.In (a) the 12-sublattice structure is highlighted.(c) Ordered moment S along the length of the cluster.