Gapless quantum spin liquid in a honeycomb Γ magnet

A family of spin–orbit coupled honeycomb Mott insulators offers a playground to search for quantum spin liquids (QSLs) via bond-dependent interactions. In candidate materials, a symmetric off-diagonal Γ term, close cousin of Kitaev interaction, has emerged as another source of frustration that is essential for complete understanding of these systems. However, the ground state of honeycomb Γ model remains elusive, with a suggested zigzag magnetic order. Here we attempt to resolve the puzzle by perturbing the Γ region with a staggered Heisenberg interaction which favours the zigzag ordering. Despite such favour, we find a wide disordered region inclusive of the Γ limit in the phase diagram. Further, this phase exhibits a vanishing energy gap, a collapse of excitation spectrum, and a logarithmic entanglement entropy scaling on long cylinders, indicating a gapless QSL. Other quantities such as plaquette-plaquette correlation are also discussed.


INTRODUCTION
The ongoing search for exotic magnetic states in highly frustrated antiferromagnets [1][2][3][4][5][6] has been extended to a new class of correlated materials with a two-dimensional honeycomb structure 7-10 and its three-dimensional variants 11 . It is suggested that bond-dependent interactions could be realized in the spin-orbit coupled Mott insulators with the aforementioned lattice geometry 12 . In particular, the Kitaev honeycomb model exhibits a novel Kitaev quantum spin liquid (QSL), which hosts fractionalized Majorana fermions and flux excitations 13 . Realization of Kitaev interaction in real materials was first proposed in iridates [14][15][16] , and then turned toward α-RuCl 3 in which Ru 3+ ions are arranged in a honeycomb lattice and carry effective spin-1/2 particles 7,9 . Although α-RuCl 3 displays long-range zigzag magnetic order at low temperature [17][18][19][20][21][22] , it is argued to be proximate to the Kitaev QSL owing to the broad continuum of magnetic excitations identified in Raman scattering 23,24 and inelastic neutron scattering 9,[25][26][27] .
In spite of massive research efforts, it has been challenging to determine exchange parameters of the proposed spin Hamiltonian for α-RuCl 3 (see refs. 28,29 and references therein). However, there is a broad consensus on a sizable off-diagonal Γ interaction 30,31 , which is antiferromagnetic (AFM) and is potentially comparable to the celebrated Kitaev interaction 26,32,33 . Crucially, it is shown that the Γ interaction could help enhance the mass gap of Majorana fermions 34 and is responsible for the strongly anisotropic responses to the magnetic field observed in α-RuCl 3 provided that the Landé g-factor anisotropy is modest 28,30,35 . In contrast to the Kitaev model 13 , analytical solution of the honeycomb Γ model has not been found yet 36 . Previous classical studies have demonstrated that its ground state is a classical spin liquid 37 followed by a flux-ordered spin liquid, which is stabilized in a finite temperature window 38 . Given the infinite classical ground-state degeneracy 37 , determining the precise quantum nature of Γ model is nontrivial, and existing numerical works have already led to conflicting results. Parallel works by exact diagonalization 39 and density-matrix renormalization group (DMRG) study of a cylinder with a width of three unit cells 40 both claim that the ground state is a nonmagnetic phase. A variational Monte Carlo simulation, on the other hand, suggests that it is a zigzag order 41 . Furthermore, a recent study proposes that it is a nematic paramagnet that spontaneously breaks the lattice rotational symmetry 42 .
In this work, we study a model which consists of the Γ term and of a staggered Heisenberg (J) interaction along the bonds, dubbed the bond-modulatedJ-Γ model (see Eq. 1). Depending on the sign ofJ, it could either favor the zigzag order (J > 0) or stripy order (J < 0). If the ground state of Γ model is a zigzag ordered phase, then the zigzag order protruding from the pure Γ limit should compete with and survive up to a finite ferromagneticJ interaction. Otherwise, there will be an intermediate phase sandwiched between the two magnetically ordered states. Thus, this model works as a virtuous arena to clarify the debates by unfolding the competing states, although it is not a description of any particular material. By employing the DMRG method on both finite cylinders with circumferences of up to 10 sites and C 3symmetric hexagonal clusters 43,44 , we identify a disordered state in between. This phase manifests characters of a gapless QSL including a dense excitation spectrum, logarithmic entanglement entropy scaling, and short-range plaquette-plaquette correlation. The pure Γ limit belongs to this QSL and is separated from the zigzag order by a first-order transition.

Model
The Hamiltonian of the bond-modulatedJ-Γ model reads where S γ i (γ = x, y, z) is the γ-component of a spin-1/2 operator at site i, and α and β are the two other bonds on a honeycomb lattice. η γ = 1 for the bond 〈ij〉 γ along the horizontal direction and equals to −1 otherwise (see Fig. 1).J and Γ are parameterized using ϑ ∈ [0, π] so as toJ ¼ cos ϑ and Γ ¼ sin ϑ ð! 0Þ.
In what follows, we carry out a hierarchical study of Eq. (1) to provide multi-faceted evidences of the gapless QSL nature of Γ magnet. We start by mapping out the classical phase diagram via the parallel tempering Monte Carlo simulation 45,46 , and conclude that the ground state of Γ model sits exactly at the classical transition point on the verge of the zigzag phase. This makes sense because the zigzag ordering belongs to the macroscopic ground-state manifold of the classical Γ model 37 . Next, we show the energy reduction and sublattice magnetization within the linear spin-wave analysis (for a review, see ref. 47 ). Afterwards, we present a quantum phase diagram obtained by large-scale DMRG calculations on various distinct cluster geometries.
Classical phase diagram Figure 2 shows the classical phase diagram of Eq. (1) obtained by Monte Carlo simulations 45,46 , coincide exactly with the subsequent results of energy optimization method (see "Methods" section). Due to the bond-modulated η γ -term, the conventional zigzag and stripy orderings perpendicular to the Z bonds are induced when ϑ/π is 0 or 1, respectively. By introducing AFM Γ interaction, the ground state becomes more competitive, triggering the possibility of other magnetic orderings in the moderate interaction regime. The ground-state energy e g = E g /(NS 2 ) (E g is the total energy) is shown in Fig. 2a, while selected spin configurations of the corresponding phases are depicted in Fig. 2b-d. In the phase diagram, the leftmost is the zigzag order with e zz g ¼ Àð2Γ þ 3JÞ=2 and its magnetic moment direction is n ½111. The rightmost is occupied by the stripy order with e st g ¼ ÀðΓ À 3JÞ=2. Its spins are perpendicular to n ½111, but could vary freely in the plane spanned by e a ½112 and e b ½110, showing an emergent continuous symmetry. Further, an extensive intermediate region appears in between. It is dominated by a so-called mixed phase in which the AFM order and two twining zigzag orders are degenerate with energy e mixed g ¼ Àð2Γ ÀJÞ=2. Here, twining zigzag orders refer to the other two zigzag orders whose spin orientations are different from the one shown in Fig. 2b (for spin configurations, see Supplementary Note 1). The zigzag-mixed transition takes place exactly at ϑ cl t;l =π ¼ 0:5, reflecting the classical spin liquid of the Γ model 37 . There is no direct transition between the mixed phase and the stripy phase expected to occur at ϑ cl t;r =π = 1 À 1 π atan 2 ≈ 0.6476. Instead, a noncollinear phase (see Fig. 2d) with e g = À ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi f appears in a narrow window of ϑ/π that is less than 0.02, see inset of Fig. 2a.

Spin-wave theory
To understand the role played by quantum fluctuations and for the sake of comparison with the DMRG results later, we have performed the linear spin-wave calculation 47 based on the quadratic Hamiltonian . Therefore, the zigzag order could only survive for AFMJ (i.e., ϑ/π < 0.50), beyond which the magnon branch becomes imaginary and should be terminated by a transition. The magnon spectra for the representative stripy order with ϑ/ π = 0.75 are shown in Fig. 3a where the wave vector q is parameterized in units of (h, k) as q = ð 2π a1 h; 2π a2 kÞ 15 . The spectra are symmetric with the middle of the Γ-M line, so the Γ and M points are equivalent. Due to the emergent continuous symmetry of the classical stripy order, the magnon spectra are gapless. In the presence of quantum fluctuations, however, the degeneracy is lifted via order-by-disorder mechanism 48 , selecting two of them that are either parallel or antiparallel to e b axis (see inset of Fig. 3b). To illustrate it, we firstly define the quantum energy correction 2 ω qυ ðϕÞ, where ϕ is the angle in the e ae b plane 49 . For ϑ/π = 0.75, which is deep in the stripy order, we show ΔE(ϕ) vs. ϕ in the inset of Fig. 3b. The energy correction has its minima at ϕ = π/2 or 3π/2, corresponding to the two mostly favored configurations at the quantum level. The energy barrier δE, defined as the energy difference between ΔE(0) and ΔE(π/2), is approximately 0.0175. The main panel of Fig. 3b shows energy barrier at different ϑ/π in the stripy order. When ϑ/π = 1.00 the energy barrier is zero, consistent with the gapless Goldstone modes thereof. Beyond that, the energy barrier is finite, indicating that the stripy order should also be twofold degenerate in the quantum case.
Due to the magnon instabilities of zigzag and stripy orderings, they could only exist in their classically allowed regions. Illustration of an XC6 cylinder on a honeycomb lattice. η γ is +1 (−1) for horizontal (zigzag) bonds. The insets are (left) the unit cell for the zigzag/stripy order with a 1 = 3 and a 2 ¼ ffiffiffi 3 p , (middle) the hexagonal plaquette operatorŴ p with its six sites enumerated, and (right) the X (red), Y (green), and Z (blue) bonds. The corresponding phase transitions could be illuminated by their sublattice magnetization. It is observed that magnetization of the zigzag order is almost saturated when ϑ/π < 0.4. As ϑ/π approaches 0.5, it undergoes a considerable suppression and the lowest branch of the magnetization nearly vanishes at ϑ cl t;l =π ¼ 0:50 (see Supplementary Fig. 5). The stripy order is more stable and its magnetization only has a small reduction at ϑ cl t;r =π % 0:6476. However, spin-wave energy in the mixed phase, say AFM order, is overwhelmingly higher than its neighbors. It thus implies that the genuine phase in the intermediate region should be different from its classical counterpart, imposing restrictions on the applicability of the spin-wave analysis.
Intervening magnetically disordered state As discussed, the spin-wave calculation fails in the intermediate regime, hence the quantum study is necessary. We have performed the standard DMRG computation on three XC clusters of 12 × 6 (n = 6), 16 × 8 (n = 8), and 20 × 10 (n = 10) and compute the ground-state energy e g = E g /N which is shown in Fig. 4. The energy curves in the middle are very flat, while they have two sharp downwarping when away from the middle region, leading to two well-marked kinks that are signals of first-order transitions. These discontinuous phase transitions could also be advocated by the entanglement entropy, which has a trend to jump as the system size is increased (see Supplementary Figs. 9 and 11). Therefore, the DMRG result supports an intermediate region that impedes a direct transition between the zigzag and stripy phases 50,51 . For comparison, we also depict the spin-wave energy of the zigzag order (green belt) and stripy order (blue belt) in Fig.  4. It is clearly found that there is a further energy reduction of the zigzag order beyond the linear approximation, in accord with the dramatic suppression of magnetic order parameter which will be clarified later. Strikingly, as shown in the inset of Fig. 4, the energy e g of Γ model (ϑ/π = 0.5) exhibits a nonmonotonic scaling behavior 52 , indicative of a possible periodicity as revealed in the Kitaev model 13 . At each fixed circumference n, the energy is linearly decreasing with length L x of the cylinder. By varying the circumference n from 4 to 10, the extrapolated energy has an oscillation in a window of −0.357 < e g < −0.352. Therefore, we estimate that the energy of Γ model is e g = −0.354(3) in the thermodynamic limit.
In order to unveil the nature of the intermediate region and to pin down the precise phase boundaries, we resort to the magnetic order parameter, which is defined as M N ðQÞ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi S N ðQÞ=N p where S N ðQÞ is the magnetic structure factor with Q being the ordering wavevector (see "Methods" section for definition). The zigzag phase has a peak at M point in the Brillouin zone, while the stripy phase possesses both peaks at M and M 0 points. Crucially, magnetic structure factor of the intermediate region is diffuse with a soft peak. Figure 5a displays order parameters M N (Q) of the zigzag and stripy phases on four distinct XC clusters with a circumference ranging from 4 to 10. Akin to their spin-wave results, magnetization of the zigzag and stripy phases exhibits maxima at ϑ/π ≈ 0.25 and 0.75, respectively. This implies that these magnetic orderings are most stable, when Γ term is approximately of equal strength to the Heisenberg interaction. Away from these points, quantum fluctuations are enhanced so that the magnetic ordering in the intermediate region is dramatically suppressed, followed by an algebraically decay with the circumference n (see Fig. 5c, d). After a careful inspection of the finite-size effect, we conclude that the magnetization will disappear eventually as the best fitting gives M → 0.0 for Γ model (ϑ/π = 0.5).
The entire quantum phase diagram of Eq. (1) is presented in Fig.  5b. In addition to the conventional zigzag and stripy orderings, there is a disordered phase, which is later interpreted as a QSL, is stabilized in a large region between ϑ t,l and ϑ t,r with ϑ t,l /π ≃ 0.50 and ϑ t,r /π = 0.66 (1). It should be noticed that even though the bond-modulated Heisenberg interaction has a strong tendency to favor the zigzag ordering, the ground state of Γ model remains disordered albeit the transition point ϑ t,l is so close to π/2. A more elaborative study on hexagonal clusters of N = 24 and 32 suggests that ϑ t,l /π = 0.498(1) (see Supplementary Fig. 12). For other less  sensitive perturbations such as the third-nearest-neighbor Heisenberg (J 3 ) interaction, the zigzag order could only exist for J 3 /Γ > 0.075 (see Supplementary Fig. 20). While for another off-diagonal Γ 0 interaction, the zigzag order is generated at Γ 0 =Γ< À 0:015 53 . These, in turn, confirm the robust nonmagnetic character of Γ model notwithstanding the aggressive zigzag ordering.
Gapless excitation and entropy scaling Next, we investigate the nature of the disordered phase by calculating the excitation gap and entanglement entropy. For this purpose, we target the first sixteen energy states on a 24-site hexagonal cluster by the DMRG method (see "Methods" section), and the first fifteen low-lying excitation gaps, Δ υ = E υ − E g (υ = 1−15), is shown in Fig. 6a. It can be seen that Δ 1 is vanishing small while Δ 2 survives in the zigzag/stripy phases, indicative of the doubly degenerate ground states predicted by the semi-classical analysis. In the intermediate region, the ground state is unique and the density of state in the low-energy spectrum is higher than its neighbors. Such a collapse of excitation gaps could be interpreted as a sign of gapless spectrum 54,55 . To check the behavior of the lowest excitation gap as N is varied, we focus on four selected points at ϑ/π = 0.40, 0.50, 0.60, and 0.80, under hexagonal cluster of N = 18, 24, and 32. As can be seen from Table 1, the lowest excitation gap Δ 2 of the zigzag phase (ϑ/π = 0.40) and stripy phase (ϑ/π = 0.80) is considerably large and slightly grows with the increasing of system size. For the intermediate phase, however, the lowest excitation gap Δ 1 at ϑ/π = 0.50 and 0.60 declines quickly when N changes from 18 to 32, indicating that the excitation gap tends to close eventually.
We then turn to XC cylinders which enable us to calculate the excitation gaps on large system sizes. Above all, we calculate the first fifteen excitation gaps on a XC cluster of 8 × 4 and we find that there is unlikely a big ground-state degeneracy in the intermediate region since the gap increases gradually without abrupt change (see Supplementary Fig. 8). Therefore, we only present the two lowest excitation gaps on three larger XC clusters up to 200 sites (see Fig. 6b). The gaps in the middle are rather small, indicative of a gapless region. We also take a closer look of the gaps at ϑ/π = 0.50, where different YC clusters are also adopted. For either XC or YC cluster, Δ 1,2 show a decreasing trend with the increasing of circumference n. In spite of the oscillation in value, they appear to vanish within a reasonable round-off, showing the gaplessness of the excitation spectrum in Γ model (see Supplementary Fig. 13). As a final consistency check for the gapless nature, we calculate the excitation gap on cylinder geometry of 2 × L x × L y (for geometry, see inset of Fig. 7b) with N = 2L x L y sites in total. Although the three-leg cylinder (L y = 3) is gapped, excitation gap on cylinder of L y = 4 decreases quickly with L x and is expected to disappear as L x → ∞ (see Supplementary Fig. 14). Such a strong size-dependent behavior of the excitation gap is typical of gapless systems.
The entanglement entropy has appeared as a versatile tool in diagnosing quantum critical systems described by conformal field theory. In this regard, the von Neumann entanglement  entropy is introduced and it is defined as SðlÞ ¼ ÀTrðρln ρÞ where ρ is the reduced density matrix of a subregion with length l 56 . Figure 7a shows the representative behavior of SðlÞ in the Γ model on a 2 × 16 × 4 cylinder, which contains eight sites along each column. When l is a multiply of 8, it corresponds to a neat edge-cutting where the two halves have smooth margins. The entanglement entropy is minimized and forms a lower branch as marked by solid symbols. Otherwise the subsystems will be more entangled, gaining extra entropy over the lower bound. Therefore, we shall fix l = N/2 to extract the central charge upon a series of finite cylinders. For such a critical system, it is recognized that the entanglement entropy scaling takes the form of S vN SðN=2Þ ¼ c 6 ln ð 2Lx π Þ þ c 0 where c is the central charge and c 0 is a non-universal constant 56 . In Fig. 7b, we shows the logarithmic fitting of S vN for cylinders of length L x = 8, 16, 24, and 32. It is found that S vN obeys the formula well with the fitting constants ðc; c 0 Þ % ð2:92; 1:09Þ, showing that the central charge is close to 3. An alternative fitting of the lowest branch of the entropy on each cylinder also demonstrates that c ≃ 3 (see Supplementary Fig. 16). By contrast, for the three-leg cylinder the entropy is extremely insensitive to the length (see Fig. 7b), revealing a central charge of 0. The fact that the central charge depends highly on the width (L y ) of cylinders may imply the existence of spinon Fermi surface (SFS). In this scenario, the pockets of SFS might be detected by different cuts in the Brillouin zone, and thus the central charge could vary for different L y . We note in passing that the central charge argument has also been used to explore the possible SFS in the fieldinduced gapless QSL in the Kitaev model [57][58][59] . Another possibility is a Dirac QSL 6,60 with three Dirac Fermions around M points, which is potentially consistent with the central charge 3 on fourleg cylinders. Information of the central charge on wider cylinders should be helpful to distinguish between the two scenarios. Regardless of different QSL natures, the vanishing magnetization and excitation gap in the Γ model, together with the distinct central charges on cylinders of L y = 3 (c = 0) and L y = 4 (c = 3), manifest that its ground state is likely a gapless QSL.
Flux-like density and plaquette correlation So far, we have confirmed that there is no magnetic ordering in the honeycomb Γ model, yet little is known about the lattice symmetry breaking. Very recently, there is a proposal of plaquette ordering stemming from a broken translational symmetry in the classical Γ model 38 . It is thus of interest to examine whether there is a plaquette ordering in the quantum situation. To this end, we study the hexagonal plaquette operatorŴ p and its correlation. Actually, W p also has its own merit as it can capture the associated phase transitions 33 . The six-body plaquette operator is known as 13 which is the product of spin operators on out-going bonds around a plaquette (see Fig. 1). Figure 8 a shows the flux-like density hW p i ¼ P p hŴ p i=N p where N p = N/2 is the number of hexagonal plaquette on clusters of N = 24 and 32. Starting from ϑ/π = 0.0, hW p i is zero, followed by a continuous decrease before arriving at the transition point, ϑ/π ≃ 0.50. Afterwards, it begins to increase and then surpasses the critical line to enter into the stripy phase where hW p i > 0. Recalling the quantum phase diagram shown in Fig. 5b, our result corroborates that the flux-like density could signal phase transitions. For Γ model   The lowest excitation gap at ϑ = 0.40π, 0.50π, 0.60π, and 0.80π on hexagon clusters of N = 18, 24, and 32. For the zigzag/stripy phase the lowest excitation gap is Δ 2 , while for the intermediate region the lowest excitation gap is Δ 1 .
Q. Luo et al.
we have hW p i ¼ À0:25ð2Þ, which is about a quarter of that in the Kitaev model 13 . Based on the plaquette-plaquette correlation hŴ pŴ q i, we then introduce the plaquette order parameter P Np ðQÞ via the plaquette structure factor W Np ðQÞ 38 (see "Methods" for definition). In Fig. 8b, we show that the QSL phase has a vigorous peak at Γ point in the reciprocal space but a weaker intensity at K point, signifying a perceptible plaquette correlation. Nevertheless, The fact that the strength of P Np ðKÞ goes down rapidly suggests that there is unlikely a plaquette ordering in the region, further corroborating a QSL without a broken symmetry.

DISCUSSION
Ever since the seminal proposal of the Kitaev interaction in heavy 4d/ 5d transition metal oxides 12 , which triggers the thriving research direction of Kitaev materials, tremendous efforts have been devoted to realizing the Kitaev QSL in real materials, yet hampered by the ineluctable non-Kitaev terms such as the off-diagonal Γ interaction.
Whereas the honeycomb Γ model has drawn enormous attention, its quantum nature is still under debate. To this end, we introduce a bond-modulated Heisenberg interaction to check its tendency towards probable magnetic orderings. For the magnetic phase diagram of the proposed bond-modulatedJ-Γ model, we find an intermediate region which is intervened between the zigzag and stripy phases. Though exhibiting magnetic order at the classical level, quantum fluctuations suppress such ordering since it acquires a large energy according to the spin-wave result. In the quantum case, it turns out to be disordered and is separated from its two neighbors by first-order transitions. By taking massive numerical efforts on the Γ model, we are able to confirm the following three subtle physical issues. (i) The low-energy spectrum is rather dense on a 24-site hexagonal cluster, and the lowest excitation gap goes down gradually with the expansion of cluster size. The empirical extrapolation on large cylinders up to 200 sites gives a vanishing energy gap, in line with the logarithmic behaviors of entanglement entropy. (ii) The zigzag magnetic ordering vanishes eventually, consistent with the suppression of magnetization of the zigzag order by spin-wave analysis. (iii) In the plaquette structure factor, there is a perceptible short-range plaquette correlation because of a subleading peak at K point. These findings strongly corroborate the ground state of Γ model is a gapless QSL rather than a zigzag order, despite the latter being close in energy.
We would like to mention that, due to the gapless nature and for the lack of continuous spin symmetry, it is exceedingly challenging to capture the fractionalized excitation in the proposed QSL. The flux insertion method, which pumps fractional particles from one edge to the other, is usually a promising way to elucidate the topological characters of the ground state. It is performed by adiabatically twisting boundary conditions of the Hamiltonian so that the U(1) symmetry is required. The discrete symmetries of the Γ model thus inherently hinder this trick. Actually, topological degeneracy is in general not well-defined for a gapless QSL, as different gauge sectors are closely connected due to gapless excitations. Nonetheless, the Kitaev QSL is special because fluxŴ p is a conserved quantity, allowing for the identification of different flux sectors by the vison insertion 55 . Also of note is that a recent study suggests the existence of a nematicity due to the lattice rotational symmetry breaking 42 . We would like to stress that the symmetry of Γ model itself is discrete and the asymmetrical boundary condition could cause instability on the landscape of bond energy, making it hard to determine the nematicity in the thermodynamic limit. However, the pending lattice nematicity does not alter our proposal of the QSL, because it could be accompanied by a broken lattice symmetry as reported in other theoretical models 61,62 . Despite such challenges, our work emphasizes on the inspiring and intractable quantum nature of the Γ model. Notably, the dominating Γ region could be realized in α-RuCl 3 under compression where the magnitude of Kitaev interaction is small 53 .
In short, our results provide a significant guidance to further theoretical and experimental studies on honeycomb magnets.

Density matrix renormalization group
In order to check for finite-size effects, we have performed large-scale DMRG calculations 43,44 on three kinds of cluster geometries. Firstly, the frequently used geometry is a L x × L y XCn cluster under cylindrical boundary condition (see Fig. 1). Here, X indicates the orientation of the cylinder, while n is the circumference of the cylinder. We consider even circumferences n (=L y /a 0 ) ranging from 4 to 10 lattice spacing a 0 , and use fixed ratio L x /L y = 2 unless stated explicitly otherwise. N = L x L y is the total number of spins. Secondly, we consider the honeycomb cylinder of 2 × L x × L y where L x (L y ) is the number of unit cell along e 1 ¼ ð ffiffi ffi 3 p ; 0Þ (e 2 ¼ ð1=2; ffiffi ffi 3 p =2Þ) direction (see the inset of Fig. 7). Due to the limitation of modern computational capability, we focus primarily on four-leg cylinders (L y = 4). The maximal value of L x is 32, and the total number of spins N = 2L x L y . Lastly, we also use the C 3 symmetric hexagonal clusters with N = 24 and 32 sites under full periodic boundary conditions. In all cases, we keep up to m = 3000-5000 states and perform about 12 sweeps in the calculation, so as to ensure the truncation error is smaller than 10 −6 . When targeting the first few low-lying energy levels, we diagonalize a subspace of a sparse Hermitian matrix iteratively by Davidson algorithm, with the precision of each eigenvalue maintained at a desired standard. In addition, all the targeted states are used with an equal weight to construct the reduced density matrix.
The magnetic order parameter is defined by is the total static magnetic structure factor, with Here, R i is the position of spin and Q is the ordering wavevector. We also calculate the plaquette-plaquette correlation hŴ pŴq i, whereŴ p is the hexagon plaquette operator (see Eq. 2). Likewise, we define the static plaquette structure factor where R p is the central position of each plaquette, and N p = N/2 is the number of plaquette. To eliminate the strong finite-size effect due to the identity hðŴ p Þ 2 i = 1, we introduce the plaquette order parameter (see Supplementary Note 5)

Simulation and energy optimization
We use the parallel tempering Monte Carlo simulation with the heat-bath algorithm to prevent the possible metastable state at low temperatures 45,46 .
The simulation is carried out in a temperature range with a hundred of replicas. The heat-bath algorithm is performed at given temperature, followed by a so-called thermal replicas, where configuration swaps between different temperatures are allowed with a probability according to a detailed balance condition. The simulations are performed on three XC clusters of 16 × 16, 24 × 24, and 32 × 32, under toroidal boundary condition. The resulting energy-optimized spin configuration is then served as the benchmark for the analytical calculation. The classical spin can be written as where θ i ∈ [0, π) and ϕ i ∈ [0, 2π). Taking the zigzag (zz) order and stripy (st) order for instance, their classical energy are e zz g ¼ À and e st g ¼ where the explicit form of the auxiliary function is F ðθ; ϕÞ ¼ sin 2 θ sin 2ϕ À sin 2θðsin ϕ þ cos ϕÞ: Q. Luo et al.
Mathematically, the maximum of Eq. (9) is 2 with ðθ; ϕÞ ¼ ðπ À atanð ffiffi ffi 2 p Þ; π=4Þ or ðθ; ϕÞ ¼ ðatanð ffiffi ffi 2 p Þ; 5π=4Þ. This means that the classical energy of the zigzag phase is e zz g ¼ Àð2Γ þ 3JÞ=2 with the classical magnetic direction n ¼ ½111. The mimimum of Eq. (9) is −1, and the energy of the stripy phase is e st g ¼ ÀðΓ À 3JÞ=2. Its moment direction is free to vary in a plane that is perpendicular to n, indicating an emergent continuous symmetry in the classical stripy phase. The analytical energy and spin configurations of the other magnetic phases are shown in Supplementary Note 1.

Linear spin-wave theory
We summarize the derivation of the spin wave spectra of the zigzag phase, which is one of the degenerate ground states of the classical Γ model. In the framework of linear spin-wave theory, the local spin operator S i ¼ S x i ; S y i ; S z i À Á is represented by bosonic creation and annihilation operators a i and a y i . By virtue of the Holstein-Primakoff transformation, we havẽ S Here,S n i ðS Á nÞ is the spin component along the classical ordered moment n andS ± i ðS i Á eÞ ± {½S i Á ðn eÞ are the ladder operators consisting of the orthogonal spin components, with e being an (arbitrary) unit vector perpendicular to n and satisfying the right-hand rule 47 . The spin operator is thus where τ is introduced for classical spin which is either parallel (τ = +1) or antiparallel (τ = −1) to n. The γ-component of the spin S γ τ;i ¼ S τ;i Á e γ . For the zigzag phase, we choose the following orthogonal axis, e = [112], n e ¼ ½110, and n ¼ ½111. Because of the four-sublattice (n s = 4) nature, the magnetic unit cell is taken as the rectangle of the area a 1 × a 2 with a 1 = 3a 0 and a 2 ¼ ffiffi ffi 3 p a 0 , see Fig. 1. Within the magnetic unit cell, the wave vector q could be parameterized in units of (h, k) as q = ð 2π a1 h; 2π a2 kÞ 15 . In this notation, M and M 0 points in the Brillouin zone could be rewritten as (1, 0) and (0, 1), respectively. By introducing four flavors of Holstein-Primakoff bosons and using the Fourier transformation, we arrive at the following spin-wave Hamiltonian wherex y q ¼ ða y q ; b y q ; c y q ; d y q ; a Àq ; b Àq ; c Àq ; d Àq Þ is a vector of length 2n s and H q is a 2n s × 2n s matrix of the form and The parameters in Eqs. (14) and (15) are given by where ϱ = e ıπh/3 . Since B Àq ¼ B Ã q , C Àq ¼ C Ã q , E Àq ¼ E Ã q , and D Àq;τ ¼ D Ã q;Àτ , we hence deduce thatΔ y q ¼Δ q andΛ T Àq ¼Λ y q . The quadratic Hamiltonian Eq. (12) can be diagonalized via a bosonic Bogoliubov transformation T(q). To preserve the canonical commutation rules of the bosons, it should satisfy the orthogonality relations TΣT † = T † ΣT = Σ where Σ = diag(1, −1). The eigenvalues of ΣĤ q give the magnon spectrum ΩðqÞ ¼ diagðω q;1 ; ω q;2 ; Á Á Á ; ω q;ns Þ. The spin-wave dispersions of other magnetic orderings (including the mixed phase and the noncollinear phase) are shown in Supplementary Note 2.