Gapless quantum spin liquid in a honeycomb $\Gamma$ 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 $\Gamma$ 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 $\Gamma$ model remains elusive, with a suggested zigzag magnetic order. Here we attempt to resolve the puzzle by perturbing the $\Gamma$ region with a staggered Heisenberg interaction which favours the zigzag ordering. Despite such favour, we find a wide disordered region inclusive of the $\Gamma$ 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.

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 ferromagnetic J 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 3 -symmetric 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 plaquetteplaquette correlation. The pure Γ limit belongs to this QSL and is separated from the zigzag order by a firstorder 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 groundstate 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 bondmodulated η γ -term, the conventional zigzag and stripy , and 32×32 (blue circle). The solid black line stands for the exact solution given by the energy optimization method. The classical phase diagram is then determined by kinks in energy curves. Inset: Zoom in of energy curves near ϑ cl t,r /π ≈ 0.6476. A noncollinear (NCL) phase appears in a narrow window of 0.6368 < ϑ/π < 0.6543. b and c depict the zigzag order and stripy order, respectively. d A NCL phase with a unit cell of 4 × 2.
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 /(N S 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ã [112] andb [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 FIG. 3: Spin-wave analysis of the stripy phase. a Four branches of the magnon spectra ωqυ in the stripy phase where ϑ/π = 0.75. The path along the symmetry directions in the momentum space is depicted in the inset. b Energy barrier δE between the stripy phases of different orientations in the classically allowed zone. Inset: Quantum energy correction ∆E(φ) vs angle φ suited at theã-b plane. The parameter is fixed to ϑ/π = 0.75, which is marked as a black hexagram in the main panel.

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 spinwave calculation 47 based on the quadratic Hamiltonian · · · is the Nambu spinor, and M q is a 2 × 2 block matrix (see "Methods" and Supplementary Note 2). There are four spin-wave dispersion branches ω qυ (υ = 1-4) for the four-sublattice (n s = 4) zigzag and stripy orderings. In the zigzag order, there exists a magnon gap ∆ at M point in the Brillouin zone (see inset of Fig. 3a). When approaching Γ limit, J/Γ 1, the lowest magnon branch is softened and the gap vanishes as ∆/Γ √ 30 3 J /Γ. 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 par- allel or antiparallel tob axis (see inset of Fig. 3b). To illustrate it, we firstly define the quantum energy cor- where φ is the angle in theã-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(π/2) and ∆E(0), 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. 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. S5). 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 groundstate 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 Fig. S9 and Fig. S11). Therefore, the DMRG result supports an intermediate region that impedes the 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) = S N (Q)/N where S N (Q) is the magnetic structure factor with Q being the ordering wavevector (see "Methods" 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 points. Crucially, magnetic structure factor of the intermediate region is diffuse with a soft peak. Fig. 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 and Fig. 5d). 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 bondmodulated 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. S12). 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. S20). While for another off-diagonal Γ interaction, the zigzag order is generated at Γ /Γ < −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 lowlying 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 Tab. I, 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. S8). 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 roundoff, showing the gaplessness of the excitation spectrum in Γ model (see Supplementary Fig. S13). 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. S14). 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 where c is the central charge and c 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 ) ≈ (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. S16). 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 four-leg 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,Ŵ 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).

FIG. 8:
Flux-like density and plaquette correlation. a Flux-like density W p and b plaquette order parameter P(Q) for N = 24 (red triangular) and 32 (blue circle). The inset exhibits the plaquette structure factor of Γ model, which has a relatively weak peak at K point (corner of the Brillouin zone) in the reciprocal space. Figure 8a shows the flux-like density W p = p Ŵ p /N p where N p = N/2 is the number of hexagonal plaquette on clusters of N = 24 and 32. Starting from ϑ/π = 0.0, W p 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 W p > 0. Recalling the quantum phase diagram shown in Fig. 5b, our result corroborates that the fluxlike density could signal phase transitions. For Γ model we have W p = −0.25 (2), which is about a quarter of that in the Kitaev model 13 . Based on the plaquetteplaquette correlation Ŵ pŴq , 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 lowenergy 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 ther-modynamic 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 Fig. 7). Due 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 cluster with N = 24 or 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 M N (Q) = 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 Ŵ pŴq 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 (Ŵ p ) 2 = 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 and where the explicit form of the auxiliary function is F(θ, φ) = sin 2 θ sin 2φ − sin 2θ(sin φ + cos φ). (9) Mathematically, the maximum of Eq. (9) is 2 with (θ, φ) = (π−atan( √ 2), π/4) or (θ, φ) = (atan( √ 2), 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 for 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 † i . By virtue of the Holstein-Primakoff transformation, we havẽ 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 = √ 3a 0 , see Fig. 1a. 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 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 spinwave Hamiltonian is a vector of length 2n s andĤ q is a 2n s × 2n s matrix of the form and∆ The parameters in Eq. (14) and Eq. (15) are given by 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 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.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding authors upon reasonable request. and X.W. checked the calculations. All authors together discussed the numerical details and drafted the article.

ADDITIONAL INFORMATION
Supplementary information is available in the online version of the paper.
Correspondence and requests for materials should be addressed to J.Z. or X.W.
Supplementary Note 1: Classical energy and spin configurations

The mixed antiferromagnetic-twining zigzag phase
Due to the bond-modulated η γ (±1) term, theJ-Γ model shown in Eq. (1) in the main text does not posses C 6 rotational symmetry. This leads to a discrepancy among the three conventional configurations of the zigzag ordering with different orientions. We find that the other two twining zigzag orders (see Fig. S1(b) and (c)) have higher energy than the standard zigzag order (see Fig. 2(b) of the main text) when ϑ/π ∈ [0, 1/2). When ϑ/π is slightly larger than 1/2, the twining zigzag orders overcome the latter and become the ground state. Interestingly, the antiferromagnetic (AFM) order (see Fig. S1(a)) has the same energy and contributes to the degenerate manifolds. Moreover, due to the two nonequivalent sites per unit cell, each configuration is two-fold degenerate. Consequently, we conclude that the mixed phase has six-fold degenerate ground states. For the AFM ordering out of the mixed phase, its energy is given by where F(θ, φ) = sin 2 θ sin 2φ − sin 2θ(sin φ + cos φ). The function F(θ, φ) is plotted in Fig. S2, which has a maximum of 2. Since the translation of φ by π does not change the magnitude of the function F(θ, φ), we then obtain e mixed g = −(2Γ −J)/2 with (θ, φ) = (atan( √ 2), π/4) or (θ, φ) = (π − atan( √ 2), 5π/4). In this case the classical magnetic direction of the AFM ordering is n = [111].

The noncollinear phase
When θ/π ≈ 0.65 there is a noncollinear (NCL) phase consisting of two kinds of spins (or four if we consider that two of them are anti-parallel to their partners) which are neither (anti-)parallel nor perpendicular in the classical Supplementary Figure S2: Illustration of function F(θ, φ) in the (θ, φ) parameter space.
For example, if θ = 3π/4 − ψ 0 /2 and φ = π/4 where ψ 0 = atan Γ 4J , we have the classical energy We also note that the angles between the two kinds of spins are −ψ 0 or its supplementary angle π + ψ 0 . Since ψ 0 is ϑ-dependent, the polar angle θ also varies with ϑ. The fascinating character of the noncollinear phase is that it may also possess other ground states with larger unit cell. For example, we find such a ground state whose (enlarged) unit cell has 64 lattice sites, see Fig. S3(b).
Supplementary Figure S3: Configurations of the noncollinear phase in the real space. The size of the unit cell marked by the yellow shadow is (a) 4 × 2 and (b) 8 × 8 (enlarged).

Supplementary Note 2: Quantum fluctuations in the Linear spin-wave theory
Here, we go beyond the classical level by considering the quantum fluctuations to find out where the nonmagnetic state may appear in the phase diagram. To this end, we utilize the linear spin-wave theory where each 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 † i . We adopt some of the notations used by Janssen and Vojta 1 . By virtue of the Holstein-Primakoff transformation, 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. The spin operator is thus where τ is introduced for classical spin which is either parallel (τ = +1) or antiparallel (τ = −1) to n. For each γ-component S γ τ,i = S i · e γ , we have where ξ γ = e · e γ , η γ = (n × e) · e γ , and ζ γ = n · e γ . Before going into detail, we shall define an auxiliary function

four-sublattice stripy order
Like the zigzag order, the stripy order also has four sublattices. Following a very similar procedure shown in "Method: Linear spin-wave theory" in the main text, we get the spin-wave Hamiltonian for the stripy order as is a vector of length 2n s andĤ q is a 2n s × 2n s matrix of the form Those matrix entries are given by where δ x = − a1 6 , a2 2 , δ y = − a1 6 , − a2 2 , and δ z = a1 3 , 0 . Here, a 1 × a 2 is the unit cell of the zigzag/stripy order with a 1 = 3 and a 2 = √ 3, see Fig. 1 of the main text. For simplicity, the wave vector q is parameterized in units of (h, k) as q = 2π a1 h, 2π a2 k 2 . Actually, there is an emergent continuous U (1) symmetry for the stripy order. The spins are perpendicular tõ c [111], but could vary freely in the plane spanned byã [112] andb [110]. As demonstrated in Fig. 3(b) in the main text, these degeneracy is lifted by quantum fluctuations, leading to a two-fold degenerate ground state which is either parallel or antiparallel tob. Thus, we choose n =b as the magnetically ordered moment direction and e =c as the arbitrary unit vector, then we find that the explicit formula in Eq. (S12) is given by (S13) where s = 1+2 √ 2ı 3 and = e ıπh/3 . The quadratic Hamiltonian Eq. (S9) can be diagonalized via a bosonic Bogoliubov transformation 1 , where Ω(q) = diag (ω q,1 , ω q,2 , · · · , ω q,ns ). The transformation matrix satisfies the orthogonality relations T ΣT † = T † ΣT = Σ where Σ = diag (1, −1). The spectrum of the Hamiltonian can be obtained by the eigenvalue equation is normalized with respect to the inner product involving the matrix Σ, i.e., n(q)|Σ|n(q) = 1 with n(q)| ≡ |n(q) † , then the columns of the matrix T (q) are given by the two vectors |n(q) . In addition to the spin-wave dispersion relations which are usually of prime interest, there are also two other important quantities which can be easily obtained using the spin-wave calculation. Namely, (i) the value of the total ordered moment M per site, and (ii) the total energy per site ε. For the classical moment M , it is straightforwardly to get 1 where v * (υ) −q denotes the lower half of the normalized υ-th eigenvector occurring in Eq. (S15), with positive energy ω q,υ . The momentum integral is over all wavevectors q = (q x , q y ) in the Brillouin zone. Likewise, the spin wave energy ε is given by 1 The AFM order has two sublattices in the honeycomb lattice. Without loss of generality we can assume that τ = +1 for the A-sublattice and τ = −1 for the B-sublattice. Quite directly, we can obtain (S19) For the AFM ordering, the spins are found to be along the [111] direction. So we choose the following crystalline axis, e = a [112], n × e = b [110], and n = c * [111]. In this case we have where ω = e 2πı/3 . In light of above equations we find that γ 0,q = e ıqδx + e ıqδy − e ıqδz γ 0,q = 1 3 e ıqδx + e ıqδy + e ıqδz γ 1,q = 1 3 ω −1 e ıqδx + ωe ıqδy + e ıqδz .
3. Spin-wave dispersion ωqυ, energy ε, and sublattice magnetization M In this section, we begin by calculating the spin-wave dispersions for the selected points in the zigzag phase (0.00 ≤ ϑ/π < 0.50), the mixed phase of AFM ordering and twining zigzag ordering (0.50 ≤ ϑ/π < 0.6368), and the noncolinear (NCL) phase (0.6368 ≤ ϑ/π < 0.6543). The spin-wave dispersion of the stripy phase, which shows an order-by-disorder mechanism, is shown in the main text. The spin-wave dispersions at ϑ/π = 0.25 (zigzag), ϑ/π = 0.60 (AFM), ϑ/π = 0.60 (twining zigzag), and ϑ/π = 0.6450 (NCL) are shown in Fig. S4, and the path in the reciprocal space is depicted in the inset for each subplot. For the zigzag ordering (see Fig. S4(a)), it is gapped with a magnon gap ∆ of 2.12. Generally, we have ∆ = 3(Γ 2 + 2ΓJ) when ϑ 0.40π. As can be seen from Fig. S4(b) and (c), the magnon gap is relatively small in the mixed phase of the AFM ordering and twining zigzag ordering. In the NCL phase, its unit cell contains 8 sites (cf. Fig. S3(a)). Therefore, the dispersion shown in Fig. S4(d) contains eight branches and the minimal gap is around 0.25.
We then turn to the spin-wave energy ε and sublattice magnetization M of the magnetically ordered states. In Fig. S5(a), the top black solid lines is the classical energy for all the phases. The transitions between the neighboring phases are of first order because of the kinks in the energy curve. The bottom yellow star line represents the quantum energy on a 24-site hexagonal cluster. The spin-wave energy of the zigzag phase (green triangle), AFM phase (red circle), and stripy phase (blue square) are calculated in their classically allowed region. In the zigzag phase, the spin-wave energy correction ∆E (when compared to the classical energy) is the largest at ϑ/π = 0.50. Here, the classical energy is -0.25 while the spin-wave energy is -0.33, showing an energy correlation of 0.08. We note that the quantum ground-state energy at the Γ limit is estimated to be -0.354(3) (see Fig. 4 of the main text), and there is a large quantum fluctuation due the classical ground-state degeneracy. In contrast, in the intermediate region, the energy correction for the AFM phase is very small, indicating that the AFM order is unlikely the true ground state at the quantum level. We emphasize here that this phenomenon is directly related to the QSL phase by large-scale DMRG calculation. In the stripy phase, the spin-wave energy is very close to the quantum result.
Supplementary Figure S5: (a) The spin-wave energy for zigzag order (green), AFM order (red), and stripy order (blue). The black line is the classical energy (S is set to be 1/2) illustrated in the main text, and the yellow star line is the quantum energy on a 24-site hexagonal cluster. (b) The four branches of the sublattice magnetization Mυ (υ = 1 − 4) for the zigzag order and stripy order. The averaged magnetizationM (thick black line) is also shown for comparison. Figure S5(b) shows the sublattice magnetization M υ (υ = 1 − 4) of the zigzag and stripy phases. Since both phases contain four sites in their unit cells, there are four different branches of M υ . Deep into the zigzag phase, the four branches have very close values, but the difference becomes pronounced at the boundaries. In the Γ limit where ϑ/π → 0.50, the lowest branch M 1 is dramatically suppressed and it tends to vanish. Since the quantum fluctuation has a strong impact on the lowest branch, it thus suggests that the ground state of Γ model is likely a QSL in the quantum situation. However, we note that the averaged magnetizationM = (M 1 + M 2 + M 3 + M 4 )/4 is nonzero, with a value ofM ≈ 0.28. We will go back to this point later. In the stripy phase, the lowest branch M 1 is relatively smaller than the remaining branches and it is also insensitive to the interaction except near the phase boundary. This phenomenon may relate to the fact the lowest branch is gapless due to the emergent continuous symmetry. In contrast, the second branch M 2 is more amenable to reveal the magnetization, and it share the same behavior as the quantum case.
Before ending this section, we wish to get further spin-wave signatures of the QSL in the pure Γ model. The classical Γ model is known to have an infinite number of degeneracy, with the two-sublattice AFM ordering, four-sublattice zigzag ordering, six-sublattice vortex-like 120 • ordering, and many large-unit-cell (LUC) ordering such as 18-site ordering. In the quantum level, some of them are not favored (e.g., the AFM order, which has a higher spin-wave energy), but a subset of them (whose number should be large) are competitive. For example, the spin configurations of the zigzag ordering (peaks at M point) and 18-site ordering (peaks at 2M/3 point) are shown in Fig. S6. For both of the two, all the spins have the same weight of three components, i,e., |S x i | = |S y i | = |S z i | = S/ √ 3. As a result, the zigzag ordering and the 18-site ordering have the same spin-wave energy and the same magnetizationM 0.28. Since the zigzag ordering and the 18-site ordering have different peaks in the reciprocal space, the linear superposition of the two leads to a reduced magnetization ofM =M / √ 2 ≈ 0.20. In addition, by an alternative rearrangement of the spins, larger LUC orderings with 36-site and 48-site unit cells could be created, and the superposition of them will further weaken the magnetization, giving rise to a nonmagnetic phase eventually. In this sense, the spin-wave calculation is consistent with our large-scale DMRG calculation.
Supplementary Figure S6: (a) Configuration of the zigzag ordering in the real space. The spins are parameterized bŷ S = S (sin θ cos φ, sin θ sin φ, cos θ). Here, θ is represented by the color (see colormap) while φ is represented by the orientation of the arrow in the plane. The shape of the unit cell (with 4 sites in total) is marked by the yellow shadow. The right panel is the static structure factor in the momentum space with peaks at M points. (b) Configuration of one of the 18-site orderings in the real space. The shape of the unit cell (with 18 sites in total) is marked by the yellow shadow. The right panel is the static structure factor in the momentum space with peaks at 2M/3 points.
Supplementary Note 3: The bond-modulatedJ-Γ model: cylinder vs hexagonal cluster 1. XC cylinder: precision, low-lying excitation, and entanglement entropy In the current computational capability, the number of block states m in the two-dimensional (2D) DMRG calculation is limited to a few thousands for the Hamiltonian without U (1) symmetry. Therefore, it is crucial to do a proper finite-size scaling of the measured quantities, such as energy and order parameter, with respect to m. Taking the 12 × 6 XC cylinder as an example, we start by keeping m = 500 states in the warming-up process. Then we continue to sweep by increasing the block states to 1000, 2000, and 3000, respectively. For each m kept we perform 4 sweeps and up to 12 sweeps are applied. The numbers of sweeps are extended to 16∼20 times in the gapless region and the block state m is increased to 4000 occasionally. Figure S7 shows the extrapolation of the energy for the gapped zigzag phase (ϑ/π= 0.40) and the gapless QSL phase with ϑ/π= 0.50 (Γ limit) and ϑ/π = 0.60. The energy in Fig. S7(a) is almost m independent, showing that the energy converges very quickly in the gapped phase. In the QSL phase, the energy at m = 1000, 2000, and 3000 shows a linear scaling of 1/m. Typically, the energy difference between the E g (m = 3000) and the extrapolated one E g (m → ∞) is not very big. For example, when ϑ/π= 0.50, the two are -25.40831921 and -25.40866029, respectively. The energy per-site is -0.35289332 and -0.35289806, with an even smaller difference. On the other hand, it is essential to check the ground-state degeneracy on XC cylinders, the number of which could be large due to the possible gapless edge excitation. To get the first few dozens of target states, we diagonalize a sparse Hermitian matrix by Davidson algorithm. Suppose that we already have the first k eigenvectors, we then extend the subspace and perform the iteration until some criterions are met, and the (k + 1)-th eigenvalue could be obtained with the default precision. We emphasize that in the DMRG calculation, we need to use all the targeted states to construct the reduced density matrix. In our implementation, all the target states are used with equal weight. In principle, this method is general and can be used to target a large number of low-energy states. But due to the truncation error, there should be a balance between the system sizes and the number of target states to guarantee the numerical accuracy. For a cylinder with ∼50 sites, the first 20-30 states could be targeted precisely with an error bar of O 10 −6 or less. In our calculation, we calculate the first sixteen energy levels on a small system size of 8 × 4 cylinder, which enables us to get the low-lying energy accurately. The excitation gap ∆ υ = E υ − E 0 with υ = 0 − 15 is shown in Fig. S8. For the zigzag order (ϑ = 0.4π), the ground state is indeed doubly degenerate with a large energy gap of ∼ 0.35 (Note: we note that this value is somewhat small due to small system size and cylinder boundary condition. As shown in Fig. 6a in the main text, excitation gap for larger cylinder is ∼ 0.65, which is consistent with the PBC case shown in Tab. 1 in the main text). For the intermediate phase at ϑ = 0.5π and ϑ = 0.6π, the excitation gap is successive increasing without large abrupt change. Thus, there is unlikely a large number of ground-state degeneracy.
In Fig. 4 of the main text, we show the ground-state energy on XC clusters of 12 × 6, 16 × 8, and 20 × 10. Here, Fig. S9 show the von Neumann entropy S on the same cylinders. It is observed that the entropy in the gapped zigzag and stripy phases are lower than that in the intermediate phase. With the increasing of the total sites, the entropy in the intermediate phase increases, consistent with the conclusion that the intermediate phase is gapless. In addition, the entropy has drops at ϑ l /π 0.50 and ϑ r /π ≈ 0.66, indicating that both of the transitions are of first order.

Hexagonal cylinder: Zigzag-QSL transition
Due to the bond-modulated Heisenberg (J) interaction which favors the zigzag ordering, the Zigzag-QSL transition near ϑ l /π 0.50 is more intricate. To determine the transition type and the transition point accurately, we have done extra calculation on hexagonal clusters of N = 18, 24, and 32 under full periodic conditions (for geometries, see Fig. S10). We note that although these clusters have less sites than the cylinder cases, they do not have the boundary effect and thus show clearer tendencies of physical quantities as N is increased. To begin with, we calculate the energy and von Neumann entropy on the 18-site cluster. It is found that there is a kink in the energy e g and a sharp jump in the entropy S (not shown), supporting the first-order transition. Nevertheless, since the 18-site cluster does not match with the zigzag ordering whose unit cell is 4, this result is more or less unjust.
Supplementary Figure S10: The lighted regions with yellow frames mark hexagonal clusters of (a) N = 18, (b) N = 24, and (c) N = 32. The black solid dots represent real sites whose number should be N . The types of bonds are X-bond (red), Y-bond (green), and Z-bond (blue). In (a) and (b), the dotted vertical lines represent the partition of the systems into two equal halves with N/2 sites in each part. In (c), the blue and cyan dotted lines represent [l, r]-partitions of [11,21] and [19,13], respectively, with four Z-bonds been cut.
In what follows, we resort to the 24-site and 32-site clusters, which are accommodative with the zigzag phase. First of all, since the entropy is amenable to signify the transition type, we thus calculate the von Neumann entropy S by cutting the same number of bonds. As can be seen from Fig. S10(b) and S10(c), the red cut in the 24-site case and the blue and cyan cuts in the 32-site case meet this requirement where four Z-Type bonds are cut, and the entropy are shown in Fig. S11(a). In the transition region where 0.48 < ϑ/π < 0.51, the entropy for the 24-site cluster increases from 2.4984 to 3.6065, an increment of 1.1081. By contrast, the increment is ∼1.80 for the 32-site case. This indicates that the difference of the entropy in the zigzag phase and the QSL phase should be enlarged as N is increased, and there should be a jump at the transition point as N → ∞. In addition, we also calculate the order parameter M (M) of the zigzag ordering for the two system sizes. As can be seen in Fig. S11(b), the order parameter has a tendency to jump with the increasing of N . In summary, the results of entropy and order parameter both favor the first-order Zigzag-QSL transition. In the main text, we estimate that the transition point ϑ l /π is around but slightly smaller than 0.50. It is the staggered-like Heisenberg (J) interaction that enhances the zigzag phase and pushes the transition point to the Γ limit. To make a reasonable and precise transition point, we turn to the energy derivatives on hexagonal clusters of N = 24 and N = 32. For a finite-size system, there is a peak in the second-order energy derivative. As the system size is increased, the peak will behave as a Gaussian-like wave-packet ( Dirac δ function) for the first-order transition or will diverge for continuous phase transition. The ground-state energy is calculated with an increment of 0.001 and 0.005 (ϑ/π) for the N = 24 and 32 cases, respectively. The energy is extrapolated by a spline-extrapolation method so as to make the energy derivatives be smooth. The first-order (red) and second-order (blue) energy derivatives are shown in Fig. S12. For the N = 24 case (left panel), the peak locates at 0.4985, while it is 0.4980 for the 32-site case (right panel). The transition point is slowly shifted to the left and we expect it will not change too much as N is further increased. Thus, we estimate the transition point as ϑ/π = 0.498 (1), and the pure Γ limit falls in the intermediate phase.
Supplementary Figure S12: (a) The first-order (in red, left) and second-order (in blue, right) energy derivatives of a 24-site hexagonal cluster. The peak position locates at ϑ/π = 0.4985. (b) The same as (a) but for a 32-site cluster. The peak position locates at ϑ/π = 0.4980. Supplementary Note 4: The honeycomb Γ model: energy gap scaling and entanglement entropy scaling

energy gap scaling
In a recent study of the one-dimensional (1D) Γ chain by us 3 , the ground-state energy e g has been shown to have an intimate relation to the boundary condition. Notably, e g exhibits a six-site periodicity with the length of L x . We recall that the gap for vison excitations in the Kitaev honeycomb model shows a three-period structure, as pointed out by Kitaev 4 . We thus speculate that the energy spectrum of the two-dimensional (2D) Γ model should also own a similar but more complicated periodicity.
The unusual size-dependent behavior has an awful impact on the energy and energy gaps shown in Fig. S13, making an accurate extrapolation to the thermodynamic limit intractable. We have calculated the low-lying energy on both XC clusters (blue square) of 8 × 4, 12 × 6, 16 × 8, and 20 × 10, and also YC clusters (red circle) of 8 × 4, 12 × 6, and 16 × 8, see Fig. S13. The XC cluster shown in Fig. 1 of the main text has zigzag open edges, while the YC cluster is a 90 • rotation of the XC cluster and has armchair open edges. In both cases, the energy has an oscillation, and we estimate the ground-state energy is -0.354(3) (more details are presented in Fig. 4 of the main text). As shown in Fig. S13(b), there are similar oscillations in the energy gaps ∆ 1 (open) and ∆ 2 (filled). Because of the overall downward trend in the gaps with the increasing of the system size, we estimate that the lowest gap should be ∆ = 0.00(1) in the thermodynamic limit. This result is in line with the dense energy spectrum on a 24-site hexagonal cluster shown in Fig. 6a of the main text), favoring the gapless ground state in the Γ model. Next, we go beyond the XC/YC cylinders with a fixed L x /L y = 2, and study the energy gap of the Γ model ranging from 1D chain to two-leg honeycomb ladder, and also towards a series of 2×L x ×L y tori of L y = 3 or 4. The geometry of the latter is shown in the inset of Fig. 7 in the main text.
• To begin with, for the 1D isotropic Γ chain, its ground state is found to be a gapless Luttinger liquid 3,5 with emergent SU (2) symmetry.
• Furthermore, we consider a two-leg honeycomb ladder which is a rung-alternating coupling of two isotropic Γ chain. It is a stripe of a honeycomb lattice along its zigzag edges and only contains L x /2 Z bonds. We find that there is a unique ground state with energy E 0 under PBC, followed by a triplet excited state with energy E 1 . There seems to be a continuous spectrum afterwards and the lowest branch has a energy of E 4 . Figure S14(a) shows the energy gap of ∆ 1 = E 1 −E 0 and ∆ 4 = E 4 −E 0 , which go down as L x increased. After an extrapolation of ∆ 4 we find that ∆ 4 < 0.004, which seems to close for long enough ladder.
• Moreover, we study the energy gaps for three-and four-leg tori. Here, PBCs are imposed on both directions so as to remove the possible edge excitations. For L y = 3, we perform the calculation on four different tori with L x = 3, 4, 5, and 6, and find that the gap is around ∼ 0.11. However, for L y = 4, the gap goes down quickly from 0.09 when L x = 3 to 0.015 when L x = 5, see Fig. S14(b). We thus infer that the gap for L y = 4 most probably vanishes as L x → ∞. The results on different clusters are summarized in the table below. We find that the energy gap has a strong cluster dependence, and could vanishes at several cases. We note that this is not the typical character of a gapped system whose gap is usual very stable. In this regard, it is another evidence for the gaplessness of Γ model.
Supplementary Table S1: Energy gap of the pure Γ model under a 1D chain, two-leg honeycomb ladder, and also 2×Lx ×Ly tori of Ly = 3 or 4.

entanglement entropy scaling
As shown in the last subsection, the three-leg cylinder is found to be gapped in the pure Γ model. Such a gapped ground state could also be checked by the entanglement entropy. Figure S15 shows the representative behavior of S(l) on a 2 × 18 × 3 cylinder, which contains six sites along each column. When l is a multiply of 6, 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 red symbols. It is clearly found that the lower branch is very flat in the middle region, in accordance with a gapped system with a central charge of 0.
By contrast, the four-leg cylinder is found to be gapless. For the quasi-one-dimensional conformal invariant critical system under open boundary condition, it is established that the entanglement entropy obeys the following formula 6 where c is the central charge and c is a model-dependent fitting constant. l x is the number of the columns of the subsystem and L x is the length of the cylinder. In Fig. 7 of the main text, we take l x = L x /2 and fit the central charge as S = c 6 ln 2Lx π + c . In that case we find that (c, c ) ≈ (2.92, 1.09), showing that the central charge is approximately 3. Here, instead, we turn to fit the central charge on each individual cylinder according to the formula Eq. (S23). Figure S16 shows the fitting of central charge on cylinders of 2 × 16 × 4, 2 × 24 × 4, and 2 × 32 × 4. We find that the best fitting values are 2.85, 2.88, and 2.90, respectively. Therefore, the central charge converges to 3 for long enough four-leg cylinder.
In the end, we want to address that the significant difference between L y = 3 and 4 suggests an unusual way from multi-leg ladder towards 2D limit. Besides, it may imply that the ground state is likely to own spinon Fermi surface Nevertheless, the central charge 3 on four-leg cylinders may also match with the possibility of a Dirac QSL where three Dirac Fermions are located around M points. We note that central charge on wider cylinders of L y = 5 and 6 should be useful to further clarify this issue. Actually, in the finite-size DMRG calculation, the width of a long cylinder (L x 50) is usually limited to three or four unit-cell, as the precision is not satisfactory for wider cylinders. So it is unclear currently how does the central charge vary with further increasing of width L y . We speculate that this problem might be studied by the variational Monte Carlo calculation on a mean-field Hamiltonian constructed out of Abrikosov fermions or Majorana fermions 9 or the infinite DMRG calculation on the correlation length spectrum 10 . We think that this is an exciting research direction in the future.

Supplementary Note 5: Plaquette order parameter
In the Kitaev honeycomb model, the hexagonal plaquette operatorŴ p = 2 6 S x 1 S y 2 S z 3 S x 4 S y 5 S z 6 commutates with the model andŴ p = ±1 4 . For the Γ model as well as the general bond-modulatedJ-Γ model, [H,Ŵ p ] = 0, soŴ p is no longer a conserved quantity. However, the flux-like density W p = p Ŵ p /N p where N p = N/2 is the number of hexagonal plaquette can be tremendously useful and informative to distinguish different phases.
We calculate the plaquette-plaquette correlation Ŵ pŴq , and define the static plaquette structure factor 11 , where R p is the central position of each plaquette. We find that there is a dominating peak in the Γ point and also a subleading peak at K point of the Brillouin zone. This fact implies that there is no translational symmetry breaking in the honeycomb lattice but with perceptible plaquette correlation. We thus define the plaquette order parameter as P Np = W Np (Q)/N p . The results on the hexagonal clusters of N = 24 and 32 are shown in Fig. S18(a). Due to the dominating contribution from the trivial identity (Ŵ p ) 2 = 1 12 , P Np has a considerable finite-size value, which is approximately 1/ N p , see the horizonal lines in Fig. S18. Therefore, the bared plaquette order parameter is formally defined as The result of the bared plaquette order parameter is shown in Fig. 7 of the main text. It can be found that there is a decreasing plaquette correlation around Γ model. Besides, there is a negative flux-like density W p = −0.25 (2) for Γ model (see Fig. S18(b)). For the quantum bond-modulatedJ-Γ model, there are three distinct phases, including a zigzag order, a stripy order, and also a QSL. The ground state of Γ model belongs to the QSL phase but locates very close to the transition point between the zigzag order and the QSL. The selected contour plots of the SMSF for the three phases are shown in Fig. S19(a)-(c). While the zigzag and stripy phases peak at M and/or M points, the magnetic order at ϑ/π = 0.5 is tiny, and a subleading peak locating at X point in the Brillouin zone appears. This peak could be enhanced by negative third-NN J 3 interaction.   S20 shows the evolution of the magnetic orders versus J 3 in a rather wide region. It could be found that the ferromagnetic (FM) and AFM J 3 model tend to select the FM phase and zigzag phase, respectively, as their ground states on the perturbation of AFM Γ interaction. Between the two, the maximum of M N (Q) appears at X point of the Brillouin zone (see Fig. 3a of the main text). At J 3 = 0, M N (Q) at M becomes comparable to that at X. However, as can be seen from the inset which shows the first derivative of M N (M) versus J 3 , the peak locates at a tiny but nonzero J 3,t ≈ 0.075. This provides further evidence that the ground state of the Γ model, in which J 3 is zero, is not the zigzag ordering.