Half-magnetization plateau in a Heisenberg antiferromagnet on a triangular lattice

TL;DR

This paper explores the phase diagram of a 2D Heisenberg antiferromagnet on a triangular lattice, revealing a half-magnetization plateau due to quantum fluctuations and competing ordered states with distinct symmetry breakings.

Contribution

It introduces a detailed phase diagram considering next-nearest-neighbor interactions, highlighting the emergence of a half-magnetization plateau from quantum order-by-disorder effects.

Findings

Identification of a half-magnetization plateau at 1/2 saturation.

Discovery of competing four-sublattice quantum states with different symmetry breakings.

Analysis of phase transitions between gapped and gapless states.

Abstract

We present the phase diagram in a magnetic field of a 2D isotropic Heisenberg antiferromagnet on a triangular lattice. We consider spin-$S$ model with nearest-neighbor ($J_1$) and next-nearest-neighbor ($J_2$) interactions. We focus on the range of $1/8<J_2/J_1<1$, where the ordered states are different from those in the model with only nearest neighbor exchange. A classical ground state in this range has four sublattices and is infinitely degenerate in any field. The actual order is then determined by quantum fluctuations via "order from disorder" phenomenon. We argue that the phase diagram is rich due to competition between competing four-sublattice quantum states which break either $\mathbb{Z}_3$ orientational symmetry or $\mathbb{Z}_4$ sublattice symmetry. At small and high fields, the ground state is a $\mathbb{Z}_3$-breaking canted stripe state, but at intermediate fields the ordered states break $\mathbb{Z}_4$ sublattice symmetry. The most noticeable of such states is "three up, one down" state in which spins in three sublattices are directed along the field and in one sublattice opposite to the field. Such a state breaks no continuous symmetry and has gapped excitations. As the consequence, magnetization has a plateau at exactly one half of the saturation value. We identify gapless states, which border the "three up, one down" state and discuss the transitions between these states and the canted stripe state.

Tables

Table 1. Table 1: Results for the boundaries of UUUD state for different J 2 / J 1 subscript 𝐽 2 subscript 𝐽 1 J_{2}/J_{1} (see SM for details of calculations). The UUUD state is stable in the range h l < h < h u subscript ℎ 𝑙 ℎ subscript ℎ 𝑢 h_{l}<h<h_{u} , where h l = h s ​ a ​ t / 2 − δ 1 subscript ℎ 𝑙 subscript ℎ 𝑠 𝑎 𝑡 2 subscript 𝛿 1 h_{l}=h_{sat}/2-\delta_{1} and h u = h s ​ a ​ t / 2 + δ 2 subscript ℎ 𝑢 subscript ℎ 𝑠 𝑎 𝑡 2 subscript 𝛿 2 h_{u}=h_{sat}/2+\delta_{2} .

$J_2/J_1$	$1/8$	$1/4$	$1/3$	$1/2$	$1$
${\delta }_1\left(1/S\right)$	0.46	0.15	0.11	0.11	0.28
${\delta }_2\left(1/S\right)$	1.2	0.80	0.75	0.75	1.09

Equations

$$\mathcal{H}$$

$$E_{\Phi }/N=-\mu \sum _{i=1,2,3}|{\Phi }_i{|}^2+\frac{1}{2}{\Gamma }_1\sum _{i=1,2,3}|{\Phi }_i{|}^4$$

$$+{\Gamma }_2(|{\Phi }_1{|}^2|{\Phi }_2{|}^2+|{\Phi }_1{|}^2|{\Phi }_3{|}^2+|{\Phi }_2{|}^2|{\Phi }_3{|}^2)$$

$$+{\Gamma }_3({\Phi }_1^2{\Phi }_2^2+{\Phi }_2^2{\Phi }_3^2+{\Phi }_3^2{\Phi }_1^2+h.c.)$$

$${\Gamma }_2-{\Gamma }_1$$

$${\Gamma }_3$$

$$E_{\Delta }/N=$$

$$+\frac{1}{2}K|{\Delta }_1^2+{\Delta }_2^2{|}^2$$

$$K=$$

$$-\frac{{\beta }_K}{S}.$$

$${\mathcal{H}}_{uuud}=$$

$$+(-2h_0+h)(a_{\mathbf{k}}^{†}{a}_{\mathbf{k}}+b_{\mathbf{k}}^{†}{b}_{\mathbf{k}}+c_{\mathbf{k}}^{†}{c}_{\mathbf{k}})+(6h_0-h)d_{\mathbf{k}}^{†}{d}_{\mathbf{k}}\}$$

$$H_{\mathbf{k}}^{(2)}=S\left(\begin{array}{cccc}-2h_0+h& 2{\xi }_{\mathbf{k}}^c& 2{\xi }_{\mathbf{k}}^b& 2{\xi }_{\mathbf{k}}^a\\ 2{\xi }_{\mathbf{k}}^c& -2h_0+h& 2{\xi }_{\mathbf{k}}^a& 2{\xi }_{\mathbf{k}}^b\\ 2{\xi }_{\mathbf{k}}^b& 2{\xi }_{\mathbf{k}}^a& -2h_0+h& 2{\xi }_{\mathbf{k}}^c\\ 2{\xi }_{\mathbf{k}}^a& 2{\xi }_{\mathbf{k}}^b& 2{\xi }_{\mathbf{k}}^c& 6h_0-h\end{array}\right),$$

$$\left(\begin{array}{c}{a}_{\mathbf{k}}\\ b_{\mathbf{k}}\\ c_{\mathbf{k}}\\ d_{-\mathbf{k}}^{†}\end{array}\right)=\left(\begin{array}{cccc}\frac{1}{\sqrt{2}}& \frac{-1}{\sqrt{6}}& \frac{1}{\sqrt{2}}& \frac{-1}{\sqrt{6}}\\ \frac{-1}{\sqrt{2}}& \frac{-1}{\sqrt{6}}& \frac{1}{\sqrt{2}}& \frac{-1}{\sqrt{6}}\\ 0& \frac{2}{\sqrt{6}}& \frac{1}{\sqrt{2}}& \frac{-1}{\sqrt{6}}\\ 0& 0& \frac{-1}{\sqrt{2}}& \frac{3}{\sqrt{6}}\end{array}\right)\left(\begin{array}{c}{e}_{\mathbf{k}}\\ f_{\mathbf{k}}\\ {\bar{c}}_{\mathbf{k}}\\ {\varphi }_{-\mathbf{k}}^{†}\end{array}\right),$$

$${\mathcal{H}}_0=$$

$${\mathcal{H}}^{(4)}=$$

$$\displaystyle+a^{\dagger}_{\gamma,1}d^{\dagger}_{2}d^{\dagger}_{3}d_{1+2+3})+h.c.\big{]}+2\xi^{\gamma}_{1-2}\big{(}a^{\dagger}_{\alpha,1}a_{\alpha,2}a^{\dagger}_{\beta,3}a_{\beta,1+3-2}-a^{\dagger}_{\gamma,1}a_{\gamma,2}d^{\dagger}_{3}d_{1+3-2}\big{)}\}.$$

$$⟨a_{\alpha }^{†}{a}_{\beta }{⟩}_{0,\mathbf{k}}=⟨(T_{{\alpha }^{\prime }\alpha }{\tilde{a}}_{{\alpha }^{\prime }}^{†}+T_{4\alpha }\tilde{d})(T_{{\beta }^{\prime }\beta }{\tilde{a}}_{{\beta }^{\prime }}+T_{4\beta }{\tilde{d}}^{†}){⟩}_{0,\mathbf{k}}=\frac{1}{N}\sum _{\mathbf{k}\in B.Z.}{T}_{4\alpha }{T}_{4\beta }\equiv ⟨T_{4\alpha }{T}_{4\beta }{⟩}_{\mathbf{k}}$$

$$\delta {\mathcal{H}}_0=S\{{\delta }_1(e_0^{†}{e}_0+f_0^{†}{f}_0)+{\delta }_2{\varphi }_0^{†}{\varphi }_0\}.$$

$$⟨a_0⟩=\frac{-1}{\sqrt{6}}⟨{\varphi }_0^{†}⟩,⟨b_0⟩=\frac{-1}{\sqrt{6}}⟨{\varphi }_0^{†}⟩,⟨c_0⟩=\frac{-1}{\sqrt{6}}⟨{\varphi }_0^{†}⟩,⟨d_0^{†}⟩=\frac{3}{\sqrt{6}}⟨{\varphi }_0^{†}⟩$$

$$\displaystyle\Delta_{\alpha}=\frac{1}{2\sqrt{3}}\big{(}\Delta_{+}e^{i\boldsymbol{Q}\cdot\boldsymbol{R}_{\alpha}}+\Delta_{-}e^{-i\boldsymbol{Q}\cdot\boldsymbol{R}_{\alpha}}\big{)}=\Delta_{1}\boldsymbol{v}_{1\alpha}+\Delta_{2}\boldsymbol{v}_{2\alpha}$$

$$e_{\mathbf{k}}=\sqrt{N}{\Delta }_1{\delta }_{\mathbf{k},0}+{\tilde{e}}_{\mathbf{k}},f_{\mathbf{k}}=\sqrt{N}{\Delta }_2{\delta }_{\mathbf{k},0}+{\tilde{f}}_{\mathbf{k}}$$

$$a_{\mathbf{k}}=$$

$$c_{\mathbf{k}}=$$

$$E_{\Delta }/N=$$

$${\mathcal{H}}^{(4)}=\frac{1}{N}\sum _{i,j,l,k}{V}_{i,j,l,k}({\mathbf{k}}_1+\mathbf{q},{\mathbf{k}}_2-\mathbf{q},{\mathbf{k}}_1,{\mathbf{k}}_2){\psi }_{i,{\mathbf{k}}_1+\mathbf{q}}^{†}{\psi }_{j,{\mathbf{k}}_2-\mathbf{q}}^{†}{\psi }_{l,{\mathbf{k}}_1}{\psi }_{k,{\mathbf{k}}_2}.$$

$$S_{\mathbf{r}}^z(\mathbf{z})=S-a_{\mathbf{r}}^{†}{a}_{\mathbf{r}},S_{\mathbf{r}}^{+}(\mathbf{z})=\sqrt{2S}\sqrt{1-\frac{a_{\mathbf{r}}^{†}{a}_{\mathbf{r}}}{2S}}{a}_{\mathbf{r}},S_{\mathbf{r}}^{-}(\mathbf{z})=\sqrt{2S}{a}_{\mathbf{r}}^{†}\sqrt{1-\frac{a_{\mathbf{r}}^{†}{a}_{\mathbf{r}}}{2S}}$$

$$S_{\mathbf{r}}^{+}=\sqrt{2S}(1-\frac{1}{4S}(1+\frac{1}{8S}+\frac{1}{32S^2}+...)a_{\mathbf{r}}^{†}{a}_{\mathbf{r}})a_{\mathbf{r}}+\mathcal{O}(a^5)$$

$$\delta {\mathcal{H}}^{(4)}=$$

$$\displaystyle+d^{\dagger}_{1}a^{\dagger}_{\gamma,2}a^{\dagger}_{\gamma,3}a_{\gamma,1+2+3}+a^{\dagger}_{\gamma,1}d^{\dagger}_{2}d^{\dagger}_{3}d_{1+2+3})+h.c.\big{]}\}$$

$$-iV({\mathbf{k}}_1+\mathbf{q},{\mathbf{k}}_2-\mathbf{q},{\mathbf{k}}_1,{\mathbf{k}}_2)$$

$$\displaystyle V^{(0)}(\boldsymbol{k}_{1}+\boldsymbol{k},\boldsymbol{k}_{2}-\boldsymbol{k},\boldsymbol{k}_{1},\boldsymbol{k}_{2})V^{(0)}(\boldsymbol{k}_{1}+\boldsymbol{q},\boldsymbol{k}_{2}-\boldsymbol{q},\boldsymbol{k}_{1}+\boldsymbol{k},\boldsymbol{k}_{2}-\boldsymbol{k})\big{)}$$

$$V^{(1)}=$$

$$=$$

$${\int }_{\omega }\frac{\mathrm{d}\omega }{2\pi }\frac{1}{(\omega -{\omega }_{{\mathbf{k}}_1})(-\omega -{\omega }_{{\mathbf{k}}_2})}=\frac{-i}{{\omega }_{{\mathbf{k}}_1}+{\omega }_{{\mathbf{k}}_2}}.$$

$${\mathcal{H}}_0=S\sum _{\mathbf{k}}\{k^2(e_{\mathbf{k}}^{†}{e}_{\mathbf{k}}+f_{\mathbf{k}}^{†}{f}_{\mathbf{k}})+(2k^2+\delta h){\varphi }_{\mathbf{k}}^{†}{\varphi }_{\mathbf{k}}\},$$

Full text

**Introduction ** Recent experimental and theoretical advances renewed the interest in the physics of frustrated spin systems. In many of these systems the classical ground state is infinitely degenerate, and the actual ground state spin configuration is selected by quantum fluctuations (the “order from disorder” phenomenon). The resulting ground state is often rather unconventional and in several cases displays a non-monotonic behavior of magnetization in an applied field, with kinks, jumps, and plateaus Coldea et al. (2002); Tokiwa et al. (2006); Takigawa and Mila (2011); Starykh (2015); Wosnitza et al. (2016); Starykh et al. (2014); Yamamoto et al. (2017). The most known example of such behavior is in the case of a two-dimensional (2D) quantum antiferromagnet on a triangular lattice with nearest-neighbor exchange $J_{1}$ Chubukov and Golosov (1991); Coletta et al. (2016). Classically, all spin configurations, which satisfy $\boldsymbol{S}_{r}+\boldsymbol{S}_{r+\boldsymbol{\delta}_{1}}+\boldsymbol{S}_{r+\boldsymbol{\delta}_{2}}=\boldsymbol{h}S/(3J_{1})$ for each triad of neighboring spins, have the same ground state energy. Quantum fluctuations lift the degeneracy and select a set of three coplanar configurations, between which the systems transforms upon increasing field. The middle configuration, which exists at $h$ around $1/3$ of the saturation field $h_{sat}=9J_{1}$ ($h$ scaled), is a collinear state with two spins up (U) and one spin down (D) in every elementary triangle (an UUD state). In such a state only a discrete $\mathbb{Z}_{3}$ symmetry is broken (one spin in a triad is selected to be antiparallel to a field), and, as a result, all excitations are gapped and the magnetization has a plateau at exactly one-third of the saturation value Chubukov and Golosov (1991); Nikuni and Shiba (1993); Griset et al. (2011); Chen et al. (2013); Coletta et al. (2016). This plateau has been observed in Cs2CuBr4 Ono et al. (2003); Honecker et al. (2004); Ono et al. (2005); Tsujii et al. (2007); Fortune et al. (2009) and in Ba3CoSb2O9 Susuki et al. (2013). An UUD state survives in a finite range of perturbations, like the spatial anisotropy of the exchange interaction Nikuni and Shiba (1993); Alicea et al. (2005, 2009); Tay and Motrunich (2010); Griset et al. (2011); Chen et al. (2013); Chubukov and Starykh (2013), multiple-spin ring exchange Kubo and Momoi (1997), and the next nearest neighbor exchange Ye and Chubukov (2017), as long as the perturbations are not strong enough to close the minimal excitation gap in the UUD state. What replaces the UUD state at larger perturbations has been the subject of intensive research over the last several years Alicea et al. (2009); Tay and Motrunich (2010); Chen et al. (2013); Chubukov and Starykh (2013); Coletta et al. (2013); Veillette et al. (2005); Veillette and Chalker (2006); Nikuni and Shiba (1993); Ueda and Totsuka (2009); Griset et al. (2011)

In this communication we study spin-$S$ Heisenberg antiferronmagnet on a triangular lattice with nearest ($J_{1}$) and second-nearest ($J_{2}$) exchange interaction. Previous studies have found that the UUD phase and other two three-sublattice coplanar ground states in a field are immune to $J_{2}$ up to $J_{2}/J_{1}<1/8$. At larger $J_{2}$, however, the set of classical ground states changes discontinuously from three-sublattice configurations to four sublattice ones, in which four spins on two neighboring triads satisfy $\boldsymbol{S}_{r}+\boldsymbol{S}_{r+\boldsymbol{\delta}_{1}}+\boldsymbol{S}_{r+\boldsymbol{\delta}_{2}}+\boldsymbol{S}_{r+\boldsymbol{\delta}_{3}}=\boldsymbol{h}S/(2(J_{1}+J_{2}))$ (see Fig. 3). This condition does not uniquely specify spin order, even at zero field. The selection of the order by quantum fluctuations at $h=0$ has been analyzed by various means Jolicoeur et al. (1990); Li et al. (2015); Hu et al. (2015); Zhu and White (2015), and the consensus is that for $1/8<J_{2}/J_{1}<1$ the winner is the stripe order with ferromagnetic alinement of spins along one of three principle axes on a triangular lattice and antiferromagnetic along the other two. The same order (the canted stripe state, see Fig. 2(a)) is selected by quantum fluctuations near the saturation field, and semiclassical (large $S$) spin-wave analysis shows Ye and Chubukov (2017) that this state remains stable at all fields. It would seem natural to conjecture that this state, with monotonic magnetization $M(h)$, is the true quantum ground state for $1/8<J_{2}/J_{1}<1$ in all fields.

We argue that the phase diagram of $J_{1}-J_{2}$ model in a field is actually rather complex, with multiple phases (see Fig. 1), and the stripe order is the ground state configuration only in some range of fields and of $J_{2}/J_{1}$. For other values of $h$ and $J_{2}/J_{1}$ the ground state configurations are the co-planar states, similar to those at small $J_{2}$. In particular, around $h=h_{sat}/2$, the ground state is the UUUD state, in which spins in three sublattices are aligned along the field and in the forth sublattice opposite to the field. This spin order breaks $\mathbb{Z}_{4}$ sublattice symmetry, but doesn’t break any continuous symmetry. As a result, spin-wave excitations are gapped, and the magnetization has a plateau at exactly $1/2$ of the saturation value. We argue that the UUUD state exists for all $J_{2}$ in the interval $1/8<J_{2}/J_{1}<1$, i.e., the magnetization plateau exists for all $J_{1}-J_{2}$ systems, either at $1/3$ of the saturation value, at $J_{2}/J_{1}<1/8$, or at $1/2$ of the saturation value, at $1/8<J_{2}/J_{1}<1$. We also analyze the proximate states to the UUUD state. Above the upper critical field $h_{u}$, the UUUD state becomes unstable towards a state in which three up-spins rotate in one direction from the direction of $\boldsymbol{h}$, and the down-spin rotates in the opposite direction (see Fig. 2(b)). Below the lower critical field $h_{l}$, we found, at large $S$, a particular coplanar state, in which down-spin does not move, while three up-spins again rotate, but now one of these three spins splits from the other two (see Fig. 2(c)). A non-coplanar, chiral umbrella state Korshunov (1993); Kubo and Momoi (1997) is close in energy and may be the ground state near $h_{l}$ at smaller $S$ (see Fig. 4).

A cascade of field-induced magnetic transitions at fields below $h_{sat}/2$ has been observed in 2$H-$AgNiO2 Coldea et al. (2014); Wawrzyńska et al. (2008). It has been argued Wheeler et al. (2009) that in this material Ni2+ ions are localized and form a $S=1$ triangular lattice antiferromagnet with $J_{2}=0.15J_{1}$, single-ion easy axis anisotropy $D$, weak ferromagnetic exchange between layers. And Classical Monte-Carlo calculations for this model have found Seabra and Shannon (2010) the region of UUUD phase, whose width at $T=0$ scales with $D$. We show that in a quantum model the UUUD phase is stable in a finite range of $h$ already at $D=0$. We expect that future measurements of the magnetization in 2$H-$AgNiO2 at higher fields will be able to detect the UUUD phase and also the cascade of phases above $h_{sat}/2$. The analysis of the high-field phases will allow one to distinguish whether UUUD order is stabilized predominantly by quantum fluctuations or by single-ion anisotropy 111If the UUUD order is dominated by quantum fluctuations, one should expect to see both ${\bar{V}}$ phase and canted stripe phase at higher fields, like in Fig. 1. If UUUD order is mostly due to single-ion anisotropy, only ${\bar{V}}$ phase is present, see Ref. Seabra and Shannon (2011).

**Model and the high field phase diagram ** The $J_{1}-J_{2}$ Heisenberg antiferromagnet on a triangular lattice is described by

[TABLE]

where $\langle i,j\rangle$ and $\langle\langle i,j\rangle\rangle$ run over all the nearest and next nearest neighbor bonds. Due to the global spin-rotational symmetry, the direction of $\boldsymbol{h}$ does not matter. We choose $\boldsymbol{h}=h\,\hat{e}_{z}$, and consider the range of $1/8<J_{2}/J_{1}<1$.

The first indication that the stripe phase is not the only ground state in a field comes from the Ginzburg-Landau analysis of the order immediately below the saturation field. We said in the Introduction that this analysis yields the stripe order. This is true for all $J_{2}$ in the interval of interest, however, with one exception – $J_{2}=J_{1}/3$. To see why this $J_{2}$ is exceptional, we note that spin-wave excitations soften at $h=h_{sat}$ at three points in the Brillouin zone ($M_{1}$, $M_{2}$, $M_{3}$ in Fig. 3). To understand the order below $h_{sat}$ one then needs to introduce three condensates $\Phi_{1},\,\Phi_{2},\,\Phi_{3}$. The ground state energy in terms of $\Phi$ is:

[TABLE]

where $\mu\sim S(h_{\text{sat}}-h)$. The type of spin order that minimizes $E_{\Phi}$ depends on the interplay between the quartic coefficients $\Gamma_{i}$. In the classical limit, $\Gamma_{1}=\Gamma_{2}=8(J_{1}+J_{2}),~{}\Gamma_{3}=0$, i.e., any state from the manifold $|\Phi|_{1}^{2}+|\Phi|_{2}^{2}+|\Phi|_{3}^{2}\equiv\mu/\Gamma_{1}$ is the ground state. Quantum fluctuations lift the degeneracy. To leading order in $1/S$ we found Ye and Chubukov (2017), near $J_{2}=J_{1}/3$,

[TABLE]

where $\beta_{1,2}>0$ are numbers of order one. The logarithm $|\log(h_{\text{sat}}-h)|$ is present because of quadratic dispersion near $M$-points in Fig. 3: e.g., near ${\bf M}_{1}$, $\omega_{\boldsymbol{k}}=SJ_{1}((1+\frac{9}{2}\alpha)k_{x}^{2}+(1-\frac{3}{2}\alpha)k_{y}^{2})-\mu$, where $\boldsymbol{q}=\boldsymbol{k}+\boldsymbol{M}_{1}$ and $\alpha=J_{2}/J_{1}-1/3$. Because of the logarithm, $\Gamma_{2}>\Gamma_{1}$, A straightforward analysis then shows that only one $\Phi_{i}$ is non-zero because it costs extra energy to develop simultaneously condensates from different valleys. The resulting order is the stripe state. A selection of $\Phi_{i}$ breaks $\mathbb{Z}_{3}$ symmetry, which for the stripe state can be understood as an orientational symmetry (spins align ferromagnetically along one of the three spatial directions). However, the prefactor for the logarithm in $\Gamma_{2}-\Gamma_{1}$ in Eq. 3 is non-zero only when the dispersion is anisotropic, and it vanishes at $J_{2}=J_{1}/3$, when $\omega_{\boldsymbol{k}}$ becomes isotropic ($\alpha=0$). For this $J_{2}/J_{1}$, the sign of $\Gamma_{2}-\Gamma_{1}$ is determined by regular $1/S$ terms, along with the sign of $\Gamma_{3}$. We computed these terms and found $\Gamma_{2}-\Gamma_{1}<0$, $\Gamma_{3}<0$. As a result, at $J_{2}/J_{1}=1/3$, all three condensates emerge with equal amplitudes and relative phases 0 or $\pi$ (because $\Gamma_{3}<0$). The four choices for $(\Phi_{1},\Phi_{2},\Phi_{3})$ are $(\Phi,\Phi,\Phi),~{}(\Phi,-\Phi,-\Phi),~{}(-\Phi,\Phi,-\Phi),~{}(-\Phi,-\Phi,\Phi)$. In each of these states spins in three sublattices tilt to one direction from the field, and in one sublattice tilt to the opposite (see Fig. 1). We label such a state $\bar{\text{V}}$ by analogy with the corresponding $V$ state 222The three-sublattice $V$ state has spins in two sublattices tilt in one direction from the field, and in another sublattice to the opposite direction. at $J_{2}<J_{1}/8$ Nikuni and Shiba (1993); Griset et al. (2011); Chen et al. (2013); Coletta et al. (2016). The $\bar{\text{V}}$ state breaks $U(1)$ spin-rotational symmetry in the plane perpendicular to the field, and also breaks a $\mathbb{Z}_{4}$ sublattice symmetry by selecting a sublattice in which spin direction is different from that in other three sublattices.

Immediately below $h_{sat}$, the $\bar{\text{V}}$ state is stable in the infinitesimally small range around $J_{2}=J_{1}/3$, at $(J_{2}/J_{1}-1/3)^{2}<1/|\log(h_{sat}-h)|$. As $h$ decreases, the width grows and becomes $\mathcal{O}(1)$ at $h_{\text{sat}}-h=\mathcal{O}(1)$. The $\bar{\text{V}}$ and the stripe state break different discrete symmetries ($\mathbb{Z}_{4}$ and $\mathbb{Z}_{3}$, respectively), hence the transition between the two states is likely first order. The increase of the width of the $\bar{\text{V}}$ state with decreasing field can be understood as a generic consequence of the fact that this state is favored by regular $1/S$ terms, i.e., by quantum fluctuations at short length scales, while the stripe phase is favored by $|\log(h_{\text{sat}}-h)|$, which comes from long-wavelength fluctuations. As the magnitude of the transverse order increases with decreasing field, long wavelength fluctuations are suppressed, and $\bar{\text{V}}$ state becomes more favorable.

**Half-magnetization plateau ** As the field decreases towards $h_{sat}/2$, the $\bar{\text{V}}$ state evolves: the spin in one sublattice continuously rotates away from the field direction towards the direction antiparallel to $\boldsymbol{h}$. The spins in three other sublattices remain parallel to each other and first rotate away from the field, and then rotate back. Eventually, near $h=h_{sat}/2$, spins in the three sublattices become parallel to $\boldsymbol{h}$ and spins in the fourth sublattice become antiparallel to $\boldsymbol{h}$ (see Fig. 2(b)). Once this happens, the system enters into the new, UUUD phase. In this phase, $U(1)$ symmetry is restored (there is no sublattice spin component transverse to the field), but $Z_{4}$ symmetry is still broken. To obtain the boundaries of the UUUD phase, we compute its excitation spectrum. For this, we introduce four sets of Holstein-Primakoff (H-P) bosons and do spin-wave calculations to order $1/S$. In the classical, $S\to\infty$ limit, the spin-wave excitations are stable only at $h=h_{sat}/2$, where the spectrum consists of one gapped spin wave branch (in-phase precession of all spins around the field), and three gapless branches, with zero modes at $\Gamma$ point of the four-sublattice Brillouin zone (see Fig. 3). Quantum $1/S$ correction to spectrum, however, make it stable in a finite range of $h$ around $h_{sat}/2$. Namely, all spin-wave branches become gapped (and positive) in a range $h_{l}<h<h_{u}$, where $h_{l}=h_{sat}/2-\delta_{1}$ and $h_{u}=h_{sat}/2+\delta_{2}$. We show the details of the calculations in the Supplementary Material (SM) and present the results for $\delta_{1}$ and $\delta_{2}$ in Table 1. We found, somewhat unexpected, that the stability width of the UUUD phase is finite for all $J_{2}$ in the interval $1/8<J_{2}/J_{1}<1$. We further computed the ground state energy of the UUUD phase to order $1/S$ (classical energy plus $1/S$ corrections from zero point fluctuations), and compared with that of the stripe phase. We found that for all $J_{2}$ the energy of the UUUD state is lower. Because of this and because the UUUD state naturally emerges from the $\bar{\text{V}}$ state, we argue that the UUUD state is the true ground state near $h=h_{sat}/2$ for all $1/8<J_{2}/J_{1}<1$. As all excitations in the UUUD state are gapped, this state has magnetization fixed at exactly $1/2$ of the saturation value.

We also verified that at the upper critical field of the UUUD state, it becomes unstable towards $\bar{\text{V}}$ state. Namely, at $h=h_{u}$ one of the spin-wave branches condenses, and the condensate leads to $\langle S_{x}\rangle=a$ for spins on three up-spin sublattices, and $-3a$ for the spins on the down-spin sublattice. This result in turn implies that the $\bar{\text{V}}$ state, which started at a point $J_{2}=J_{1}/3$ at $h=h_{sat}$, extends over the whole range of $J_{2}$ near $h_{sat}/2$ (see Fig. 1).

At the lower boundary of the UUUD phase, two other spin-wave modes become unstable at the $\Gamma$ point. To determine the the spin order below $h_{l}$, we again perform Landau Free energy analysis in terms of the corresponding two complex order parameters $\Delta_{1}$ and $\Delta_{2}$. We present the details in the SM. The Free energy has the form Nandkishore et al. (2012); Venderbos et al. (2016):

[TABLE]

Classically, $\Gamma=h_{sat}/4$, $K=0$. Then $|\Delta_{1}|^{2}+|\Delta_{2}|^{2}\equiv\mu/\Gamma$, i.e. different ordered states are degenerate. Quantum fluctuations lift the degeneracy, and the result depends on the sign of $K$. If $K>0$, $\Delta_{1}=\pm\,i\Delta_{2}$. It can be checked that this gives rise to a non-coplanar umbrella state, in which the down-spin remains intact, and three up-spins split out and form a cone. Such a state breaks $U(1)\times\mathbb{Z}_{4}\times\mathbb{Z}_{2}$ symmetry. If $K<0$, the relative phase between $\Delta_{1}$ and $\Delta_{2}$ is either [math] or $\pi$, and the order is coplanar (see Fig. 4).

We computed $K$ to accuracy $1/S$. The details of calculations are presented in SM, and here we quote the result: $K$ is the sum of logarithmical, $|\log(h_{l}-h)|/S$, $\log S/S$, and non-logarithmical, $\mathcal{O}(1/S)$ terms, much like Eq. 3. The logarithmical term yields $K<0$, however the prefactor for the logarithm vanishes at $J_{2}=J_{1}/3$, and at this value of $J_{2}$ non-logarithmical terms become relevant. Near $J_{2}=J_{1}/3$, we have

[TABLE]

Where the $|\log(h_{l}-h)|/S$ term is a contribution from spin wave modes which go as $k^{2}$ at $h=h_{l}$, and $\log S/S$ term comes from another spin wave mode that softens at $h=h_{u}=h_{l}+\mathcal{O}(1/S)$. In distinction to the situation near $h_{sat}$, here we found that $K$ remains negative, even for $J_{2}/J_{1}=1/3$. This implies that the state below $h_{l}$ is a co-planar state. An umbrella state is not ruled out, however, for smaller $S$ as we computed $\beta_{K}$ in Eq. 30 at $S\gg 1$.

To determine the structure of the coplanar state below $h_{l}$ more work is actually required because for $K<0$, the Free energy to order $\Delta^{4}$ is $E_{\Delta}/N=-\mu(|\Delta_{1}|^{2}+|\Delta_{2}|^{2})+\frac{1}{2}(\Gamma-|K|)(|\Delta_{1}|^{2}+|\Delta_{2}|^{2})^{2}$, i.e., the degeneracy is not fully lifted. To select the order, one has to compute $\mathcal{O}(\Delta^{6})$ terms in the Free energy. We found (see SM for detail) that sixth-order terms select the order in which of the three up-spins two are tilting in one direction and another in the opposite direction, while the down spin remains intact (see Fig. 4). This state breaks $U(1)\times\mathbb{Z}_{4}\times\mathbb{Z}_{3}$ symmetry. It can potentially transform gradually into the stripe state, which breaks $U(1)\times\mathbb{Z}_{3}$, if the down spin begins rotating at higher deviations from $h_{l}$ and match the spin from up-triad, which is separated from the other two. Or, the transition can be first order. Either way, at small fields, the order becomes a stripe. A more complex phase diagram at low fields is expected in the presence of a single-ion anisotropy Seabra and Shannon (2011); Coldea et al. (2014).

**Conclusions ** In this work we analyzed the phase diagram of a Heisenberg antiferromagnet on a triangular lattice, with nearest and second nearest neighbor interactions ($J_{1}-J_{2}$ model), in a magnetic field. We focused on the case $1/8<J_{1}/J_{2}<1$, when semiclassical description involves four sublattice representation. We argued that the phase diagram is quite rich and contains several phases, besides the stripe state, which breaks $\mathbb{Z}_{3}$ orientational symmetry and has been detected at $h=0$ and near a saturation field. The most substantial result of our study is the identification of the UUUD phase, in which spins in the three sublattices are directed along the field, and spins in the fourth sublattice are directed opposite to the field. Such a state breaks a discrete $Z_{4}$ sublattice symmetry, but no orientational and continuous symmetry. As a result, all excitations in the UUUD phase are gapped, and the magnetization is fixed at exactly $1/2$ of the saturation value. We demonstrated that this phase is stable in a finite range of fields near $h_{sat}/2$ and is likely the true ground state of the model at all $J_{2}$ from the interval $1/8<J_{1}/J_{2}<1$. We identified gapless planar states around the UUUD phase. The one at higher fields is the $\bar{\text{V}}$ state. It breaks $U(1)\times\mathbb{Z}_{4}$ symmetry. The one at lower fields breaks $U(1)\times\mathbb{Z}_{4}\times\mathbb{Z}_{3}$ symmetry. A close competitor to this last state is a non-coplanar umbrella state. Such a state may potentially develop at a smaller $S$. We call for magnetization measurements in quasi-2D triangular-lattice antiferromagnets with $J_{2}>J_{1}/8$, best matetials with $S=1/2$, as they normally have no single-ion anisotropy, but also $S=1$ materials, like 2$H-$AgNiO2 Coldea et al. (2014); Wawrzyńska et al. (2008), to verify the existence of the plateau at a half of the saturation value of magnetization and quantum phases above this field.

**Acknowledgement ** We acknowledge with thanks useful conversations with C. Batista, A. Coldea, R. Coldea, S-W Cheong, J. Kang, N. Perkins, and O. Starykh. We are particularly thankful to R. Coldea for careful reading of the manuscript and useful comments. The work was supported by the NSF DMR-1523036.

Supplementary Material

In the Supplementary Material we present technical details of calculations, which we reported in the Manuscript. We focus on the analysis of the half-magnetization plateau state (UUUD). Calculations right below the saturation field have been presented in our earlier work Ye and Chubukov (2017).

In the formulas below, $N$ is defined as the number of sites in a given sublattice, i.e. $N=\frac{N_{tot}}{n_{subl}}$. For example, for the four sublattice states, $n_{subl}=4$, and $N=\frac{1}{4}N_{tot}$.

Appendix A A: Low energy spectrum of the UUUD state

The excitation spectrum of the UUUD state can be straightforwardly obtained by using a four-sublattice representation for three spin-up and one spin-down sublattices and introducing four sets of Holstein-Primakoff (H-P) bosons $a,~{}b,~{}c$ for the spin-up sublattices and $d$ for the spin-down sublattice. The linear spin wave Hamiltonian reads:

[TABLE]

where $h_{0}\equiv(J_{1}+J_{2})$. $\xi_{\boldsymbol{k}}^{\alpha_{i}}~{}(\alpha_{i}=a,b,c)$ are the structure factors due to the exchange interactions between sublattices. They take the form of $\xi_{\boldsymbol{k}}^{\alpha_{i}}=J_{1}\cos\boldsymbol{k}\cdot\boldsymbol{\delta}_{i}+J_{2}\cos\boldsymbol{k}\cdot\boldsymbol{l}_{i}$, where $\boldsymbol{\delta}_{2}=(1,0),\,\boldsymbol{\delta}_{1,3}=(\frac{1}{2},\mp\frac{\sqrt{3}}{2})$, $\boldsymbol{l}_{2}=(0,\sqrt{3}),\,\boldsymbol{l}_{1,3}=(\pm\frac{3}{2},\frac{\sqrt{3}}{2})$ are the nearest and next-nearest neighbor bonds respectively in unit of the lattice constant. To obtain the classical spin wave spectrum, i.e. to diagonalize the classical quadratic Hamiltonian of Eq. 6, we write $\mathcal{H}^{(2)}$ in the matrix form as: $\mathcal{H}^{(2)}=\sum_{\boldsymbol{k}}\Psi_{\boldsymbol{k}}^{\dagger}H_{\boldsymbol{k}}^{(2)}\Psi_{\boldsymbol{k}}$.

[TABLE]

where $\Psi_{\boldsymbol{k}}=\{a_{\boldsymbol{k}},b_{\boldsymbol{k}},c_{\boldsymbol{k}},d^{\dagger}_{-\boldsymbol{k}}\}^{T}$. To preserve the commutation relation of bosons, the canonical transformation $\Psi_{\boldsymbol{k}}=T_{\boldsymbol{k}}\Phi_{\boldsymbol{k}}$ satisfies $g=TgT^{\dagger}$, where $g=\text{diag}~{}(1,1,1,-1)$Colpa (1978). As a result, $T_{\boldsymbol{k}}$ that diagonalizes the Hamiltonian satisfies $T_{\boldsymbol{k}}^{-1}g\,H^{(2)}T_{\boldsymbol{k}}=g\,\Lambda_{\boldsymbol{k}}$ Ye and Chubukov (2017), where $\Lambda_{\boldsymbol{k}}=\text{diag}\,(\omega_{\tilde{a},\boldsymbol{k}},\omega_{\tilde{b},\boldsymbol{k}},\omega_{\tilde{c},\boldsymbol{k}},\omega_{\tilde{d},\boldsymbol{k}})$. $\omega_{\boldsymbol{k}}$ are the spin-wave spectrum, and $\Phi_{\boldsymbol{k}}=\{\tilde{a}_{\boldsymbol{k}},\tilde{b}_{\boldsymbol{k}},\tilde{c}_{\boldsymbol{k}},\tilde{d}^{\dagger}_{-\boldsymbol{k}}\}^{T}$ are canonical eigenmodes.

Diagonalizing the linear spin wave Hamiltonian yields three gapless and one gapped spin wave branches. The latter one, with gap $h_{\text{sat}}/2$, describes the in phase precession of UUUD state around the magnetic field. All three quadratic dispersing gapless modes soften at $\boldsymbol{k}=0$, i.e. $\Gamma$ point of the four-sublattice Brillouin zone. In the following, we focus on the spectrum at $\Gamma$ point, and analyze the stability of the UUUD state and its proximate states. $H_{\boldsymbol{k}}^{(2)}$ at the $\Gamma$ point can be diagonalized through a global rotation of basis,

[TABLE]

where $\{e_{\boldsymbol{k}},f_{\boldsymbol{k}},\phi^{\dagger}_{\boldsymbol{k}}\}$ are the low energy modes. The spectrum at $\Gamma$ reads,

[TABLE]

which indicates that classically only at $h=4h_{0}=h_{\text{sat}}/2$, the spectrum of all spin-wave branches is non-negative and the UUUD state is among the generate ground state manifold. To see if quantum fluctuations stabilize the UUUD state in a range of field near $h_{\text{sat}}/2$, corrections to the spectrum at $1/S$ order should be calculated. As there is no 3-boson interactions for a collinear state in the isotropic Heisenberg model, the $1/S$ corrections to the spectrum only come from the 4-boson interactions through the 1/S Hartree-type self-energy. The 4-boson interaction term of the UUUD state reads:

[TABLE]

$N$ is the total number of sites. For brevity, we denote $\xi_{1}\equiv\xi_{\boldsymbol{k}_{1}}$, $d_{1}\equiv d_{\boldsymbol{k}_{1}}$, and so forth. Note that $\xi_{\boldsymbol{k}}=\xi_{-\boldsymbol{k}}$. The set of $\{a_{\alpha},a_{\beta},a_{\gamma}\}$ runs over all cyclic permutations of $\{a,b,c\}$ as $\{a,b,c\}$, $\{c,a,b\}$, $\{b,c,a\}$.

The corrections to quadratic terms from $\mathcal{H}^{(4)}$ can be obtained by contracting two magnons, i.e. calculating magnon densities such as $\langle a^{\dagger}_{\alpha}a_{\beta}\rangle_{0,\boldsymbol{k}}$, $\langle\xi a^{\dagger}_{\alpha}a_{\beta}\rangle_{0,\boldsymbol{k}}$. The subscript of $\langle...\rangle_{0,\boldsymbol{k}}$ labels the quantities averaged over, e.g. [math] labels averaging over the quadratic Hamiltonian, $\boldsymbol{k}$ averages over the crystal momentum in the sublattice Brillouin zone. The density averages can be obtained by the eigenvectors (columns of matrix $T$). For example,

[TABLE]

$\langle a^{\dagger}_{\alpha}d^{\dagger}\rangle_{0,\boldsymbol{k}}=\langle T_{4\alpha}T_{44}\rangle_{\boldsymbol{k}}$, $\langle a_{\alpha}d\rangle_{0,\boldsymbol{k}}=\langle T_{1\alpha}T_{14}+T_{2\alpha}T_{24}+T_{3\alpha}T_{34}\rangle_{\boldsymbol{k}}=\langle T_{4\alpha}T_{4\beta}\rangle_{\boldsymbol{k}}$, $\langle d^{\dagger}d\rangle_{0,\boldsymbol{k}}=\langle T_{14}^{2}+T_{24}^{2}+T_{34}^{2}\rangle_{\boldsymbol{k}}=-1+\langle T_{44}^{2}\rangle_{\boldsymbol{k}}$, and so forth. $g=TgT^{\dagger}$ has been applied to simplify the above expressions. As the matrix $T_{\boldsymbol{k}}$ is a regular function of $\boldsymbol{k}$, the self-energy should also be regular near $\Gamma$. Thus it is enough to calculate the Hartree terms at $\Gamma$, where classical spin-wave modes soften. Expressing them in terms of the canonical eigenmodes $\{e_{\boldsymbol{0}},f_{\boldsymbol{0}},\phi^{\dagger}_{\boldsymbol{0}}\}$, we have

[TABLE]

where $\delta_{1},\,\delta_{2}$ are linear combinations of the averages. They sets the boundary of the UUUD state by requiring all the three modes are gapped at $\Gamma$, i.e. $h-4h_{0}+\delta_{1}>0$, $4h_{0}-h+\delta_{2}>0$. From $h>4h_{0}-\delta_{1}$, we define the lower critical field as $h_{l}=4h_{0}-\delta_{1}$, below which $\{e,\,f\}$ modes soften. From $h<4h_{0}+\delta_{2}$, we obtain the upper critical field above which $\phi$ mode softens. The numerical values of $\delta_{1},\,\delta_{2}$ at different $J_{2}/J_{1}$ are listed in Table I. As $\delta_{1},\,\delta_{2}$ are positive for all $J_{2}/J_{1}$ in the range $1/8<J_{2}/J_{1}<1$, the UUUD state is stable in this full range of $J_{2}/J_{1}$.

The pattern of transverse magnetic order above $h_{u}$ and below $h_{l}$ can be identified straightforwardly from Eq. 8. Above $h_{u}$, the $\phi$ mode develops condensate. From Eq. 8,

[TABLE]

Due to the rotation symmetry along the field, $\langle\phi_{0}\rangle=|\langle\phi_{0}\rangle|e^{i\varphi}$, $\varphi\in(0,2\pi)$. We relate the spin order with the magnon condensate through the H-P transformation, i.e. $\langle S_{\alpha,x}\rangle=\sqrt{2S}\frac{\langle a_{\alpha}\rangle+\langle a^{\dagger}_{\alpha}\rangle}{2}=-\sqrt{\rho}\cos\varphi$, $\langle S_{\alpha,y}\rangle=\sqrt{2S}\frac{\langle a_{\alpha}\rangle-\langle a^{\dagger}_{\alpha}\rangle}{2i}=-\sqrt{\rho}\sin\varphi$ for $\{a,b,c\}$ and $\langle S_{d,x}\rangle=\sqrt{2S}\frac{\langle d_{\alpha}\rangle+\langle d^{\dagger}_{\alpha}\rangle}{2}=3\sqrt{\rho}\cos\varphi$, $\langle S_{d,y}\rangle=-\sqrt{2S}\frac{\langle d_{\alpha}\rangle-\langle d^{\dagger}_{\alpha}\rangle}{2i}=3\sqrt{\rho}\sin\varphi$ for $d$ sublattice, where $\rho\equiv\frac{S}{3}|\langle\phi_{0}\rangle|^{2}$. Thus the transverse spin order has three up-spins point to one direction and the down-spin point to the opposite direction, which matches the spin order of the $\bar{\text{V}}$ state near $h_{\text{sat}}/2$. Fixing the phase of the condensate $\varphi$ specifies the plane of the $\bar{\text{V}}$ state and breaks $U(1)$ rotation symmetry around the field.

The pattern of the transverse magnetic order below $h_{l}$ is more complex due to the degeneracy of the zero modes $\{e_{\boldsymbol{0}},\,f_{\boldsymbol{0}}\}$ at $h_{l}$. And the order can be non-coplanar if the relative phase between condensates $\Delta_{1}\equiv\langle e_{0}\rangle$ and $\Delta_{2}\equiv\langle f_{0}\rangle$ is non-zero mod $\pi$. In particular, if $\Delta_{2}=\pm\,i\Delta_{1}$, the transverse orders associated with $\Delta_{1}$ and $\Delta_{2}$ takes a relative angle $\pm\pi/2$, thus the transverse order of $|\Delta_{1}|\pm i|\Delta_{2}|$ form a equilateral triangle as shown in Fig. 5 and (b), where $\pm$ interchanges $b,\,c$ sublattice labels. And the total magnetic order (transverse plus longitudinal) has non-zero chirality, i.e. an umbrella state.

Appendix B B: Symmetry constraint on the Landau free energy

We show details of obtaining the form of the Landau free energy in powers of magnon condensates $\Delta$ below $h_{l}$, the lower critical field at which the UUUD state becomes unstable. As explained at the end of Sec. A, the order below $h_{l}$ should be described by a two-component order parameter, two degenerate The Landau free energy of $\Delta$ should respect a $C_{3v}$ symmetry, i.e. threefold rotation around the center of a triangle formed by three nearest neighbor up-spins ($C_{3}$) and three reflections in the symmetry lines of the equilateral triangle ($\sigma_{v}$). Upon lowering the field, the only classical ground state configuration that respect the $C_{3v}$ symmetry is umbrella state, which breaks $\mathbb{Z}_{2}$ chiral symmetry. The transverse magnetic order of spin $\alpha$ ($\alpha=a,b,c$) can be expressed as $\Delta_{+}e^{i\boldsymbol{Q}\cdot\boldsymbol{R}_{\alpha}}$ or $\Delta_{-}e^{-i\boldsymbol{Q}\cdot\boldsymbol{R}_{\alpha}}$, where $\Delta_{\pm}$ denotes the spiral order of different chirality, $\boldsymbol{R}_{a}=(0,0)$, $\boldsymbol{R}_{b}=(-\frac{1}{2},\frac{\sqrt{3}}{2})$, $\boldsymbol{R}_{c}=(-\frac{1}{2},-\frac{\sqrt{3}}{2})$. $\{e^{i\boldsymbol{Q}\cdot\boldsymbol{R}_{\alpha}},e^{-i\boldsymbol{Q}\cdot\boldsymbol{R}_{\alpha}}\}$ forms the basis of the $\Gamma_{3}$ irreducible representation of $C_{3v}$ symmetry group. Thus a generic order parameter in $\Gamma_{3}$ representation is:

[TABLE]

where $\boldsymbol{v}_{1}=\{0,\frac{-1}{\sqrt{2}},\frac{1}{\sqrt{2}}\}$, $\boldsymbol{v}_{2}=\{\frac{2}{\sqrt{6}},\frac{-1}{\sqrt{6}},\frac{-1}{\sqrt{6}}\}$. $\{\boldsymbol{v}_{1},\boldsymbol{v}_{2}\}$ forms another basis of the $\Gamma_{3}$ representation, and it is defined in such a way that $\Delta_{1},\Delta_{2}$ is the same as defined in previously $\Delta_{1}\equiv\langle e_{0}\rangle$ and $\Delta_{2}\equiv\langle f_{0}\rangle$ up to exchanging the sublattice labels $a,\,c$. $\Delta_{\pm}$ and $\Delta_{1,2}$ are related as $\Delta_{1}=\frac{-i}{\sqrt{2}}\,(\Delta_{+}-\Delta_{-})$, $\Delta_{2}=\frac{1}{\sqrt{2}}(\Delta_{+}+\Delta_{-})$. The Landau free energy at order $\Delta_{\pm}^{2n}$ can be obtained by finding all channels that contribute to $(\Gamma_{3}\otimes\Gamma_{3})^{n}\rightarrow\mathbb{I}$.

At quadratic order ($n=1$), the only term is $|\Delta_{+}|^{2}+|\Delta_{-}|^{2}\equiv|\Delta_{1}|^{2}+|\Delta_{2}|^{2}$.

At quartic order, as $(\Gamma_{3}\otimes\Gamma_{3})^{2}=(\Gamma_{1}\oplus\Gamma_{2}\oplus\Gamma_{3})^{2}$, there are three channels contributing to identity. From $\Gamma_{1}\otimes\Gamma_{1}$, one can get $(|\Delta_{+}|^{2}+|\Delta_{-}|^{2})^{2}\equiv(|\Delta_{1}|^{2}+|\Delta_{2}|^{2})^{2}$; from $\Gamma_{2}\otimes\Gamma_{2}$, one can get $(|\Delta_{+}|^{2}-|\Delta_{-}|^{2})^{2}\equiv|\Delta_{1}^{*}\Delta_{2}-\Delta_{1}\Delta_{2}^{*}|^{2}$; from $\Gamma_{3}\otimes\Gamma_{3}$, one get $(|\Delta_{+}^{*}\Delta_{-}|^{2}+|\Delta_{-}^{*}\Delta_{+}|^{2})\equiv|\Delta_{1}^{2}+\Delta_{2}^{2}|^{2}$. Combining all three contributions, the quartic term takes the form of Eq. 4, which can break the $\mathbb{Z}_{2}$ ($\sigma_{v}$) chiral symmetry and select $\Delta_{+}$ or $\Delta_{-}$, but cannot break $\mathbb{Z}_{3}$ ($C_{3}$) symmetry.

We show that $\mathbb{Z}_{3}$ symmetry can be spontaneously broken by the six-order term. From the channel $(\Gamma_{3}\otimes\Gamma_{3})^{3}\rightarrow(\Gamma_{3})^{3}\rightarrow\Gamma_{3}\otimes\Gamma_{3}\rightarrow\mathbb{I}$, one can obtain terms of $\Delta_{+}^{3}\Delta_{-}^{*3}+\Delta_{+}^{*3}\Delta_{-}^{3}$ and $\Delta_{+}^{6}+\Delta_{-}^{6}$. In terms of $\Delta_{1}$ and $\Delta_{2}$, when their relative phase is [math] or $\pi$, the above two terms take the same form of $\Delta_{1}^{6}-15\Delta_{1}^{4}\Delta_{2}^{2}+15\Delta_{1}^{2}\Delta_{2}^{4}-\Delta_{2}^{6}$, which select either $\Delta_{1}$ or $\Delta_{2}$ and thus spontaneously breaks the $\mathbb{Z}_{3}$ symmetry.

It turns out that for a system with certain trigonal or hexagonal symmetry, i.e. $C_{6v}$, $D_{6h}$, etc., the Landau free energy of a two-component order parameter takes the same form Nandkishore et al. (2012); Venderbos et al. (2016). Technically, it’s because the multiplication table Koster (1963) for the two-dimensional irreducible representation of these symmetry groups are the same.

Appendix C C: Calculation of quartic coefficient $K$

Below $h_{l}$, both $e$ and $f$ magnon potentially softens. We split the magnon operators $e_{\boldsymbol{k}}$ and $f_{\boldsymbol{k}}$ into condensate and non-condensate as

[TABLE]

or in terms of $\{a_{\boldsymbol{k}},b_{\boldsymbol{k}},c_{\boldsymbol{k}}\}$ as

[TABLE]

From (B), the ground state energy density in powers of magnon condensates up to quartic order is:

[TABLE]

As explained in the manuscript, the condensates manifold that minimized the energy depends on the sign of $K$. In the following, we show how the sign of $K$ is obtained up to order $1/S$.

For convenience, Eq. 10 can be expressed as

[TABLE]

where $\{\psi_{i,\boldsymbol{k}}\}=\{a_{\boldsymbol{k}},b_{\boldsymbol{k}},c_{\boldsymbol{k}},d^{\dagger}_{-\boldsymbol{k}}\}$ or $\{e_{\boldsymbol{k}},f_{\boldsymbol{k}},\tilde{c}_{\boldsymbol{k}},\phi^{\dagger}_{-\boldsymbol{k}}\}$, depending on the basis of canonical modes in the context. $V_{i,j,l,k}$ is chosen such that double counting has been avoided.

The classical value for $K$ is obtained by expressing all magnon operators by the condensates $\Delta_{1},\,\Delta_{2}$ following Eq. 15 or Eq. 16. We find $K^{(0)}\equiv\frac{1}{2}V_{ffee}(\boldsymbol{0},\boldsymbol{0},\boldsymbol{0},\boldsymbol{0})=0$.

We now calculate $K^{(1)}$ from quantum corrections at order $1/S$. Only the sign of $K$ is significant, which generally doesn’t change across the critical field $h_{l}$ below which the quadratic term becomes negative. We calculate the 4-boson interaction at order $1/S$ at $h_{l}$. There are two sources of quantum corrections, one from the normal ordering of the Holstein-Primakoff bosons, another from quantum fluctuations at second order perturbation.

C.1 Corrections from normal ordering

The Holstein-Primakoff (H-P) transformation

[TABLE]

contains $\sqrt{1-a_{\boldsymbol{r}}^{\dagger}a_{\boldsymbol{r}}/2S}$. To express the Hamiltonian in terms of the H-P bosons, we expand $\sqrt{1-a_{\boldsymbol{r}}^{\dagger}a_{\boldsymbol{r}}/2S}$ in powers of the bosons. Due to the normal ordering of the bosons in the expansion, e.g. $(a^{\dagger}_{r}a_{r})^{2}=a^{\dagger}_{r}a^{\dagger}_{r}a_{r}a_{r}+a^{\dagger}_{r}a_{r}$, $S_{\boldsymbol{r}}^{+}$ can be written as:

[TABLE]

Eq. 10 is obtained keeping the leading term in $1/S$, i.e. $S_{\boldsymbol{r}}^{+}\approx\sqrt{2S}(1-\frac{1}{4S}a_{\boldsymbol{r}}^{\dagger}a_{\boldsymbol{r}})a_{\boldsymbol{r}}+\mathcal{O}(a^{5})$, etc.. To obtain the 4-boson interaction to the second order in $1/S$, we keep to order $1/S^{2}$ in Eq. 20, i.e. $S_{\boldsymbol{r}}^{+}\approx\sqrt{2S}(1-\frac{1}{4S}(1+\frac{1}{8S})a_{\boldsymbol{r}}^{\dagger}a_{\boldsymbol{r}})a_{\boldsymbol{r}}+\mathcal{O}(a^{5})$.

The normal ordering contribution to the quartic term reads:

[TABLE]

Plug Eq. 16 into Eq. 21, the normal ordering contribution to $K$ is $K^{(1)}_{a}=\frac{J_{1}+J_{2}}{24S}$.

C.2 Corrections from quantum fluctuations

Corrections to the 4-boson interaction at order $1/S$ can be expressed diagrammatically as shown in Fig. 6. Following the Feynman rules that associate each 4-boson vertex with $-iV$ and each boson propagator with $iG_{k}=\frac{i}{\omega-\omega_{\boldsymbol{k}}}$, the 4-boson interaction up to $1/S$ is

[TABLE]

As each ladder contributes a $1/S$ factor due to the magnon dispersion $\omega\sim S$, $1/S$ corrections to the 4-boson interaction and magnon dispersion is not relevant. The subscript of $G_{k}$, $\int_{k}$ is short for $k=(\omega,\boldsymbol{k})$. The indices for multiple magnon branches and the corresponding symmetry factor are suppressed for brevity. $V^{(1)}(\boldsymbol{k}_{1}+\boldsymbol{q},\boldsymbol{k}_{2}-\boldsymbol{q},\boldsymbol{k}_{1},\boldsymbol{k}_{2})$ can be expressed as

[TABLE]

To get the expression in the second line, we integrate the frequency using

[TABLE]

From the expression of Eq. 22, one can see that when $\omega_{\boldsymbol{k}_{1}+\boldsymbol{k}}+\omega_{\boldsymbol{k}_{2}-\boldsymbol{k}}\sim Sk^{2}$, the integral is logarithmical divergent if $V^{(0)}\sim\mathcal{O}(1)$. In the following, we first calculate the logarithmical part, which can be done analytically, and then calculate the regular $1/S$ part if necessary.

$K^{(1)}$ at logarithmical accuracy come from quadratic dispersive soft modes in the low energy $\{e,~{}f,~{}\phi\}$ sector (see Eq. 8). First, we diagonalize the low energy spectrum at $h=h_{l}$ as follows. In terms of $\{e,f,\phi^{\dagger}\}$, the low energy quadratic Hamiltonian can be written in two parts as $\mathcal{H}_{uuud}=\mathcal{H}_{0}+(J_{2}/J_{1}-1/3)\mathcal{H}^{\prime}$. $\mathcal{H}_{0}$ is diagonal in $\{e,f,\phi^{\dagger}\}$, while $\mathcal{H}^{\prime}$ mixes $\{e,f,\phi^{\dagger}\}$ modes.

[TABLE]

where $\delta h=\delta_{1}+\delta_{2}=\mathcal{O}(1/S)$ is the width of the UUUD state. In the matrix form, $\mathcal{H}^{\prime}=S\sum_{\boldsymbol{k}}\Psi^{\dagger}_{\boldsymbol{k}}H^{\prime}\Psi_{\boldsymbol{k}}$, $\Psi_{\boldsymbol{k}}=\{e_{\boldsymbol{k}},f_{\boldsymbol{k}},\phi^{\dagger}_{-\boldsymbol{k}}\}^{T}$,

[TABLE]

where $\varphi_{\boldsymbol{k}}=2\theta_{\boldsymbol{k}}+\pi/3$, and $\theta_{\boldsymbol{k}}$ is defined as the angle between $\boldsymbol{k}$ and the positive x-axis. When $k^{2}\ll 1/S$, the $\phi$ mode decouples from $e,~{}f$ modes that are gapless at $h_{l}$. Below $h_{l}$, the condensates develop, the quadratic dispersing gapless modes either acquire a gap that scales as $h_{l}-h$ or become linear dispersing Goldstone mode, both of which introduce a natural cutoff to the logarithm – $|\log(h_{l}-h)|/S$ as a primary contribution at logarithmic accuracy. While $\phi$ mode starts to contribute to $K^{(1)}$ significantly when $k^{2}\gtrsim 1/S$ in the integral, at order $\log S/S<<|\log(h_{l}-h)|/S$. To obtain $K$ at order $|\log(h_{l}-h)|/S$, it is enough to use the 4-boson interaction in terms of $e,~{}f$ near $\boldsymbol{k}=0$.

[TABLE]

where $h_{0}\equiv J_{1}+J_{2}$. Though there is no $e^{\dagger}e^{\dagger}ff+h.c.$ term at classical level, as $\{e,f\}$ are not the eigenmodes at $\boldsymbol{k}\neq 0$ when $J_{2}/J_{1}\neq 1/3$, $e^{\dagger}e^{\dagger}ff$ at order $|\log(h_{l}-h)|/S$ is expected to be non-zero. First, condense two incoming or outgoing magnon modes and the non-condensed ones will be served as internal propagators, $\mathcal{H}_{l}$ yields

[TABLE]

The $\boldsymbol{k}$-dependence of the interaction, not relevant at logarithmic accuracy, has been suppressed. Following Eq. 22, the leading order correction to the energy in proportion to $\Delta_{1}^{2}\bar{\Delta}_{2}^{2}+h.c.$ is:

[TABLE]

where $\langle\hat{O}\rangle_{0}$ is the average of operator $\hat{O}$ over the quadratic Hamiltonian, which is equivalent to the boson scattering problem presented diagrammatically in Fig. 6. Though different from the fully polarized state Starykh et al. (2014); Ye and Chubukov (2017), the quadratic Hamiltonian of the UUUD state gets corrections from quantum fluctuations, the calculation can be self-consistently carried out in powers of $1/S$. For our purpose, it is enough to calculate $K$ up to order $1/S$. Comparing Eq. 17 and Eq. 26, $K=-2h_{0}^{2}\frac{4}{N}\sum_{\boldsymbol{k},\boldsymbol{q}}\langle e^{\dagger}_{\boldsymbol{k}}e^{\dagger}_{-\boldsymbol{k}}f_{\boldsymbol{q}}f_{-\boldsymbol{q}}\rangle_{0}$. To calculate $\langle e^{\dagger}_{\boldsymbol{k}}e^{\dagger}_{-\boldsymbol{k}}f_{\boldsymbol{q}}f_{-\boldsymbol{q}}\rangle_{0}$, we obtain the canonical eigenmodes $u,~{}v$ at low energy through a rotation of basis:

[TABLE]

where $\tilde{\theta}_{\boldsymbol{k}}=\theta_{\boldsymbol{k}}-\pi/3$. The low energy spectrum of $u,~{}v$ modes are $\omega_{u,\boldsymbol{k}}=SJ_{1}k^{2}$, $\omega_{v,\boldsymbol{k}}=SJ_{1}(1+3\alpha)k^{2}$, where $\alpha\equiv J_{2}/J_{1}-1/3$. $\langle e^{\dagger}_{\boldsymbol{k}}e^{\dagger}_{-\boldsymbol{k}}f_{\boldsymbol{q}}f_{-\boldsymbol{q}}\rangle_{0}$ in terms of canonical eigenmodes $u$ and $v$ is:

[TABLE]

Following Eq. 22,

[TABLE]

where the factor $2$ at the front of the RHS of Eq. 29 comes from two choices to contract $u^{\dagger}u$. Similarly, we obtain other contributions in Eq. 28 and $K$ at $|\log(h_{l}-h)|/S$ accuracy reads

[TABLE]

Again $\alpha\equiv J_{2}/J_{1}-1/3$. A plot of $K^{(1)}_{\log}$ at primary logarithmical order is shown in Fig. 7.

Following a similar procedure, we obtain $K^{(1)}$ at the secondary logarithmical order $\log S/S$, and reach a similar result, i.e. $K^{(1)}<0$ when $J_{2}/J_{1}\neq 1/3$. It implies that upon lowering the field even to $h_{l}-h\gtrsim 1/S$, $K^{(1)}$ remains negative when $J_{2}/J_{1}\neq 1/3$, i.e. long wave length fluctuations select planar configuration at $J_{2}/J_{1}\neq 1/3$.

At $J_{2}/J_{1}=1/3$, the mixing part $(J_{2}/J_{1}-1/3)\mathcal{H}^{\prime}$ in $\mathcal{H}_{uuud}$ vanishes, and $e,~{}f,~{}\phi$ are canonical eigenmodes at $k^{2}$ order. Thus $K$ at both primary and secondary logarithmical order vanishes. We calculate regular $1/S$ part of $K$ from quantum fluctuations, and define it as $K^{(1)}_{b}$ (remind that $K^{(1)}_{a}$ comes from normal ordering). The calculation is straightforward but tedious. First, similar to how Eq. 25 is obtained, a complete $\mathcal{H}_{i,\boldsymbol{k}}$ is obtained by condensing two incoming or two outgoing bosons while keeping all other bosons as non-condensate operators at all $\boldsymbol{k}$. The complete expression of $\mathcal{H}_{i,\boldsymbol{k}}$ reads

[TABLE]

Second, we express the condensates $\langle a_{\alpha,0}\rangle$ in terms of $\Delta_{1},\Delta_{2}$ as defined in Eq. 16, express the non-condensed boson in terms of canonical modes $\{\tilde{a},\tilde{b},\tilde{c},\tilde{d}^{\dagger}\}$ defined below Eq. 7, and calculate the one-ladder contribution to $K$ from the four magnon branches at all $\boldsymbol{k}$. Terms of the form $\langle V^{2}\,\tilde{a}^{\dagger}_{\alpha}\,\tilde{a}^{\dagger}_{\beta}\,\tilde{a}_{\alpha}\,\tilde{a}_{\beta}\rangle_{0},~{}\langle V^{2}\,\tilde{d}^{\dagger}\,\tilde{d}^{\dagger}\,\tilde{d}\,\tilde{d}\rangle_{0}$ contribute to $K^{(1)}_{b}$, where $\tilde{a}_{\alpha},\,\tilde{a}_{\beta}$ run over $\{\tilde{a},\tilde{b},\tilde{c}\}$. We obtain $K^{(1)}_{b}=-0.37J_{1}/S$. Adding quantum corrections from normal ordering and quantum fluctuations, we have $K^{(1)}=K^{(1)}_{a}+K^{(1)}_{b}=\frac{J_{1}}{S}(\frac{(1+J_{2}/J_{1})}{24}-0.37)=\frac{-0.32J_{1}}{S}<0$.

Appendix D D: Evolution of the UUUD state below $h_{l}$

From $K<0$ for all $1/8<J_{2}/J_{1}<1$ in the limit $S\gg 1$, the relative phase between $\Delta_{1},\Delta_{2}$ is fixed to be zero $\mod\pi$, and Eq. 17 can be expressed as:

[TABLE]

However, the degeneracy in the ground state manifold is not fully lifted, and $\Delta=\cos\phi\,|\Delta_{1}|+\sin\phi\,|\Delta_{2}|$, $\phi\in(0,2\pi)$. As shown in App. B, the free energy at order $\Delta^{6}$ breaks the degeneracy and either $\Delta_{1}$ or $\Delta_{2}$ is selected depending on the sign in the sixth order term:

[TABLE]

If $C>0$, $\Delta_{2}$ is selected, and if $C<0$, $\Delta_{1}$ is selected. Obtaining the sign of $C$ requires calculating the six-boson interaction with $1/S$ correction, which requires much heavier numerical calculations than simply comparing the energy between a sequence of coplanar states that can be described by $\Delta=\cos\phi\,|\Delta_{1}|+\sin\phi\,|\Delta_{2}|$ close to $h_{\text{sat}}/2$. We find that the difference in their energy can indeed be fitted by the expression of $E^{(6)}_{C}$, and $C>0$ in this case. Thus $\Delta_{2}\equiv\langle f_{0}\rangle$ is selected. Similar to how the $\bar{\text{V}}$ state is identified above $h_{u}$ in App. A, From Eq. 8,

[TABLE]

it can be checked that the magnetic order selected by the sixth-order term has of the three up-spins two are tilting in one direction and another in the opposite direction, while the down spin remains intact (see Fig. 4). We also checked that condensates up to $\Delta^{3}$ order doesn’t couple to the $d$-magnon linearly, suggesting that down-spin remains pointing down below $h_{\text{sat}}/2$. If it is the case at even higher deviations from $h_{l}$, the transition from the coplanar state to the canted stripe state has to be first order, which is also consistent with the kink in the energy vs. $h$ plot (Fig. 7).