Evidence for field induced quantum spin liquid behavior in a spin-1/2 honeycomb magnet

One of the most important issues in modern condensed matter physics is the realization of fractionalized excitations, such as the Majorana excitations in the Kitaev quantum spin liquid. To this aim, the 3d-based Kitaev material Na2Co2TeO6 is a promising candidate whose magnetic phase diagram of B // a* contains a field-induced intermediate magnetically disordered phase within 7.5 T<|B|<10 T. The experimental observations, including the restoration of the crystalline point group symmetry in the angle-dependent torque and the coexisting magnon excitations and spinon-continuum in the inelastic neutron scattering spectrum, provide strong evidence that this disordered phase is a field induced quantum spin liquid with partially polarized spins. Our variational Monte Carlo simulation with the effective K-J1-{\Gamma}-{\Gamma}'-J3 model reproduces the experimental data and further supports this conclusion.


INTRODUCTION
Quantum spin liquids (QSLs) are exotic phases of matter resulting from competing interactions or geometric frustration.Due to the long-range quantum entanglement in the QSL ground states, interesting phenomena can arise, such as the collective excitations with fractional quantum numbers and the emergence of gauge fluctuations. 1,2−8 Meanwhile, anisotropic interactions resultant from spin-orbital coupling (SOC) have attracted more and more interests. 1,9−32 A paradigmatic example is the exactly solvable Kitaev model on the honeycomb lattice which hosts QSL ground state and Majorana-fermion-like elementary excitations. 8−31,34−38 However, owing to the existence of non-Kitaev interactions, all of these materials exhibit zigzag antiferromagnetic (AFM) order at low temperatures.6][17][18][19]23,24,36 Unlike the Ru-and Ir-materials, 19,36,38 the 3d orbitals in the Co-based materials are more compact and the contributions from the SOC channels t 2g -e g and e g -e g can weaken the Γ and Γ ′ exchanges. 24,38 On of the most representative Kitaev materials among the 3d-cobalt magnets is the Na 2 Co 2 TeO 6 (NCTO), 1,24,29,34,35,39−44 in which the honeycomb layers are formed by the magnetic Co 2+ ions surrounded by the O 2− octahedrons, Figures 1A-B.The principal reciprocal vectors a* (crystallographic vector a) direction is parallel (perpendicular) to the Co-Co bond, which corresponds to the [ 110] ([ 12 1]) direction in the spin coordinate, Supplementary Figure S1. NCTO preents a zigzag AFM order below T N ≈ 26 K and another two anomalies at T F ≈ 15 K and T* ≈ 7 K, Figures S2A-B.1,35,41−44 The zigzag order in NCTO can be easily suppressed by a magnetic field parallel to the a*-axis, leading to an intermediate field-induced magnetically disordered phase above ˜7.5 T before entering a trivial polarized phase near 10 T. 1,44 The exact nature of this intermediate phase, most intriguingly, whether it belongs to a QSL, is still illusive and deserves further investigation. 1 The study of zero-field spin-wave excitations of NCTO indicates that while the AFM thirdneighbor Heisenberg exchange interaction J 3 is fairly large, 1,29,33 the Kitaev term cannot be ignored.24,29 Several theoretical works have estimated the value of K, however, it varies from large to small, even its sign from ferromagnetic to AFM. 1,24,29,34,35,38,41,45,46 It is reasonable to expect that the SOC caused bond-dependent interactions, including the Kitaev term, play an important role in understanding the rich phase diagram of NCTO in magnetic fields, Figure1C.
In the present work, we studied the nature of the field-induced intermediate magnetically disordered phase of a single-crystal NCTO via magnetic torque and inelastic neutron scattering (INS) spectroscopy.Under low temperatures and low fields, the torque is very weak and exhibits a 2-fold (C 2 ) symmetric angular dependence, which confirms the AFM long-range order.The AFM order vanishes above 7.5 T as the lattice 6-fold (C 6 ) symmetry is restored in the angular dependence of the torque, which is verified by the disappearance of Bragg peaks at the M-point at B = 8 T, Figures 1D-E.The material enters the polarized phase at 10 T where a phase transition is observed in the differential magnetic susceptibility as well as the differential torque, Figure 1G.The region between 7.5 T and 10 T is shown to be a field-induced QSL phase with partial spin polarization and strong quantum fluctuations.With an 8 T magnetic field along the a* direction, the intensity of INS spectra at the M-point is suppressed, while gapped spin-wave bands show up at 1.5 meV˜2.5 meV and 3 meV˜4 meV in the vicinity of the Γ-point (resulting from the partial polarization of the spins) and an intense 'Λ' shape spinon continuum appears at 4 meV˜8 meV.These features are consistent with a theoretically computed dynamical structure factor of a field-induced partially polarized QSL phase.

MATERIALS AND METHODS
Sample preparation and characterization.The high-quality single crystals were grown by the flux method.The polycrystalline sample of NCTO was mixed with the flux of Na 2 O and TeO 2 in molar ratio of 1:0.5:2 and gradually heated to 900 • C at 3 • C/min in air after grinding.The sample was retained at 900 • C for 30 h, and was cooled to a temperature of 500 • C at the rate of 3 • C/h.The furnace was then shut down 1 .
Magnetization and heat capacity.The magnetization measurements were performed by using a vibrating sample magnetometer (VSM) in the physical properties measurement system (PPMS Dynacool-9 system, Quantum Design) with field up to 9 T. The heat capacity measurements were carried out using the relaxation method in another PPMS with field up to 13 T.The magnetization and heat capacity could be found in the supplementary Figure S2.
Magnetic torque.The magnetic torque measurements were carried out using piezo-resistive sensor made by Quantum Design, external bridge excitation and Lock-in amplifier readout were utilized.An oriented NCTO single crystal was mounted onto the sensor.The magnetic field was applied in the ab plane, as illustrated in Supplementary Figure S3.Both angular and magnetic field dependent torque measurements were carried out.The low temperature and magnetic field environment were provided by either a Quantum Design PPMS-9 or a top-loading 18T-320 mK system.
Inelastic neutron scattering.INS experiments were performed using the SEQUOIA time-offlight spectrometer at the Spallation Neutron Source, Oak Ridge National Laboratory, USA. 47,48out 0.559 g samples were fixed on an aluminum sheet with 3 × 6 × 0.05 cm 3 in size, and co-aligned in the (HHL) scattering plane with B // a*, Supplementary Figure S8.The sample was inserted in a liquid-helium cryostat, reaching a base temperature of T = 2 K. Measurements at 2 K with applied magnetic fields B = 0 T and 8 T were performed by rotating the sample in steps of 1 • with E i = 18 meV and choppers in high-resolution mode, yielding a full-width at halfmaximum (FWHM) elastic energy resolution of about 0.41 meV.When the magnet was removed, we collected again the INS data at 0 T and 4.9 K, which also were performed by rotating the sample in steps of 1 • with E i = 18 meV and choppers in high-resolution mode.In order to subtract the background, the INS data were collected at 90 K with or without magnet.
Variational Monte Carlo simulation.The VMC method is a variational approach using Gutzwiller projected mean field states as trial wave functions of spin models.The mean field state is obtained in the slave particle representation, where the spin operators are represented in bilinear form of fermions under a particle number constraint.The mean field parameters are not obtained self-consistently, but are treated as variational parameters whose optimal values are determined by minimizing the trial energy.The trial energy and physical quantities (including the correlation functions) of the Gutzwiller projected state are obtained using Monte Carlo simulations.

Magnetic torque.
The magnetic torque of a sample τ = µ 0 VM×B is highly sensitive to the external magnetic field B when the induced magnetization M is not aligned with B, where µ 0 denotes the permeability of the vacuum and V the volume of the sample.Therefore, the torque in a uniform B is a direct detection of the magnetic anisotropy.
where n is the unit vector along the field direction, the quantity 1 B dτ dB contains the information of the off-diagonal differential magnetic susceptibility and is thus helpful for locating the phase boundaries, Figure 1G.At low temperatures, three phase transitions can be identified by the anomalies in the field derivative ( 1 B dτ dB ) where B C1 is the transition field from the zigzag-phase to an intermediate region labeled as 'X', B C2 is the critical field from the X-region to a magnetically disordered phase, and the B C3 is the threshold of the trivial polarized phase. 1 The values of the critical field strength B C1 , B C2 and B C3 slightly vary with the angle θ, but the features of the three transitions are qualitatively unchanged.The anomalies of the 1 B dτ dB curves become weak when θ approaches 0 • (but the three critical fields are still consistent with the differential susceptibility dM/dB curve at θ = 0 • in Figure 1G).Therefore, for clarity we choose the critical fields at θ = 10.6 • to construct the temperaturefield phase diagram, Figure 1C.The critical fields obtained from specific heat measurements with B // a*(Figure S2E), are comparable with the ones obtained by the torque measurements at θ = 10.6 • .Notice that the field-dependent magnetic torque and magnetization both show obvious hysteresis near B C1 , Figure S5, indicating the transition is first-order.This is further verified by the clear hysteresis loop in the angular dependence of the torque around 6 T, Figure S4.
−51 Since a* is the easy axis, the induced magnetization M is parallel to a* if B // a* (i.e. for θ = n × 60 • , n is an integer), see Figure 2B-D.As the space group of NCTO is P6 3 22 (No.182) whose point group is D 6 , τ (θ) should exhibit a C 6 symmetry (namely 2π/6 periodicity) if there is no symmetry breaking.As shown in Figure 2A, τ (θ) only shows a C 2 symmetry for B = 3 T.This indicates a rotation symmetry breaking in NCTO (from C 6 to C 2 ) which confirms the AFM long-range order in weak magnetic fields below T N .Since the thermal fluctuations tend to melt the symmetry breaking orders, the symmetry of τ (θ) is expected to increase with increasing temperature and eventually reaches the C 6 in the paramagnetic state above T N .However, as shown in Figures 2E & H, above T N , the symmetry is still C 2 with a different orientation.A possible reason for these inconsistences is that the magnetic field is not perfectly lying in the ab-plane (the c-direction is not strictly parallel to the rotation axis, Figure S3), thus the absolute value of the angle between the field and the c-axis oscillates with a 2-fold periodicity. 52Since the effective in-plane and out-of-plane g-factors are different, g ab = 4.13 and g c = 2.3, 1 the oscillation of the field component along the c-direction results in the two-fold periodic pattern in τ (θ).
Meanwhile, strong magnetic field and quantum spin fluctuations can also suppress the zigzag order and restore the symmetry.As shown in Figures 2B-D & G, in the angle-dependent torque data the 6-fold symmetry indeed shows up above B C1 = 6 T at low temperatures with coexisting C 2 symmetry.The C 6 symmetry becomes almost perfect when the AFM order is completely suppressed at B C2 = 7.5 T. Above B C3 = 10 T, the magnitude of the torque decreases with field strength for a polarized phase with diminished quantum fluctuations.The most interesting physics falls in the region between B C2 and B C3 , a field-induced disordered state with fairly strong quantum fluctuations which is likely to be a QSL phase.Later we will provide theoretical and further experimental evidences to verify the QSL phase.The region X between B C1 and B C2 is considered as a phase with coexisting AFM and topological order, Figure S6.It should be mentioned that the τ (θ) pattern of the QSL region (above B C2 ) still does not show a strict C 6 symmetry with some mild amplitude modulation of 2π/2 period, Figures 2C & F. Those 2-period Fourier components are the same as that of the high-temperature paramagnetic state, Figure S6, hence, this should also be the issue of field alignment mentioned above.This field-induced intermediate QSL phase is supported by the Variational Monte Carlo (VMC) simulation with a K-J 1 -Γ-Γ ′ -J 3 model (J 1 is the first-neighbor Heisenberg exchange, the values of the parameters will be discussed later), where four phases are obtained with B // a*-axis including the zigzag phase, an intermediated phase with coexisting magnetic order and topological order, the filed-induced QSL phase and the polarized trivial phase, Figure S10.Especially, fixing the field's strength and varying the field's direction in the QSL phase, the induced magnetization M is parallel to B as B is along the a*or a-direction.When the field is deviated from a or a*, M contains nonzero component in a direction perpendicular to the field, which gives rise to nonzero magnetic torque.As shown in Supplementary Figure S12, the simulated magnetization indeed exhibits a 6-fold periodicity in the QSL region, which is consistent with the experimental data shown in Figures 2C & G.
Neutron scattering.To further verify the field induced QSL behavior, we performed scattering measurements at 2 K in the (HHL) plane with B // a*-direction ([K, -K, 0]) at 0 T and 8 T. As shown in Figure 1D, at zero field, the magnetic Bragg reflections can be observed at the M-point (such as [-1/2, 0, 0] and [0, -1/2, 0]), which presents the zigzag AFM order.At 8 T, these magnetic Bragg reflections at the M-points completely disappear but a new Bragg reflection appears at the Γ-point for the partial polarization, Figure 1E.More interestingly, the applied field also dramatically changes the spin excitation spectrum.
Figure 3A presents the momentum dependence of INS intensity integrated from 1.5 to 2.5 meV at 4.9 K and 0 T. The ring-shaped spectra are clearly seen around the M-points, which can be identified as magnon excitations in the zigzag ordered ground state.At the Γ-point, some excitations also show up with smaller weight compared to the M-point.Figure 3B  show up.These features are consistent with experiment in Figure 3C. 34th increasing field, the largest intensity of magnetic excitations shifts from the M-point to the Γ-point.At 8 T, around the Γ-point, a band of concave shape shows up at 1.5 meV ˜2.5 meV and another band of convex shape appears at 3 ˜4 meV.These two bands, which look like the upper lip and the lower lip, are constituted by single-particle-like magnon excitations.Away from the Γ-point, the weights of the two magnon bands decay rapidly.Instead, a large piece of continuum is observed at higher energy, indicating the existence of fractional excitations beyond the linear-spin-wave theory which predicts only two magnon bands.The continuum extends to the whole Brillouin zone, and its lower edge is overlapping with the upper-lip shaped magnon band.
Along the Γ-X-K-M-Γ path, the bright weights of the continuum in the energy range 4 ˜8 meV form a 'Λ' shape.The pattern of momentum-energy distribution and the fairly strong intensity rule out the possibility of two-magnon continuum.Thus, the continuum is most likely formed by twospinon fractionalized excitations.The coexistence of (incomplete) single-particle-like magnon bands and fractionalized continuum is the most exciting observation of the present work.From the strong continuum excitations and the C 6 symmetry of the torque, we infer that NCTO enters a field-induced QSL phase with partial spin polarization and strong quantum fluctuations at 8 T and low temperature.
The dynamical structure factor of the field-induced QSL phase obtained from VMC simulation, Figure3F, captures most of the important features of the neutron experiment.(1) Both single-particle-like magnon bands and the spinon-continuum are obtained.The magnon modes are dispersive in-gap two-spinon bound states, which form two bands.Above the magnon bands a continuum is formed by fractionalized spinons.The energy ranges of the magnon band and the continuum agree with experiment.(2) In the vicinity of the Γ-point, the two magnon bands form the shape of a lower lip and an upper lip.Similar lip-structure also exists in the linear spin-wave dispersion and is resulting from the significant J 3 interactions.Nonzero magnon weights appear at the M-and M 1 -points with energies close to zero.These features agree with the experiment.
(3) The lower edge of the continuum is overlapping with the upper magnon band.From 4 meV to 8 meV, the bright weights of the continuum form a shape of 'Λ', which qualitatively agrees with the experiment.( 4) The phase has 4-fold topological degenerate ground states on a torus as the emergence of deconfined Z 2 gauge fluctuations and the Z 2 QSL nature of the low energy physics.The deconfined Z 2 gauge charges, namely, the Majorana-fermion like spinons, give rise to the continuum spectrum in the dynamical structure and interpret the experimental weights at 4 ˜8 meV.Especially, in the intermediate field region the linear spin wave spectrum based on a fully polarized state contains imaginary part around the M-points, Figure S14C, hence, this phase is distinct from the fully polarized phase and is beyond the description of linear-spin-wave theory.
To verify the validity of the VMC computations, we performed analytic calculations using random phase approximations (RPA), Figure S15-S16, and the RPA results qualitatively agree with those of the VMC.

DISCUSSION AND CONCLUSION
In our theoretical simulation, we adopted the parameters J 1 = -1.54meV, J 3 = 1.32 meV, K = 1.408 meV, Γ = -1.32meV, and Γ ′ = 0.88 meV, which are equivalent to J 1 = 0.066 meV, J 3 = 1.32 meV, K = -3.399meV, Γ = 0.286 meV, and Γ ′ = 0.077 meV via the dual transformation.This set of parameters is adopted from the tx+ model, 46 but with enlarged J 3 and globally multiplied by a constant.8][39][40]46,53 In the 3d Cobased honeycomb geometry, the hopping integral associated with the 90 • e g -ligand hybridization plays a significant role through the large σ-type hopping process t σ P d , which is particularly relevant for the third neighbor J 3 super-exchanges in honeycomb materials. 24,38Moreover, the ratio ∆/λ between the trigonal crystal field (∆) and the SOC (λ) can regulate the spin-orbit entanglement.24 With the increasing of ∆/λ, the orbital degeneracy is lifted and the spin-orbit entanglement is suppressed.Powder INS analysis of λ = 21 meV and ∆ = 13meV with a small ratio ∆/λ ∼ 0.62 indicates that the spin and orbit are highly entangled.34 In summary, based on the magnetic torque and neutron scattering experiments, we studied the

S2. Magnetization and specific heat
The single-crystal DC magnetic susceptibility χ(T), Fig. S2A, measured by the zero-field cooling at 0.01 T and the single-crystal specific heat C P (T)/T at 0 T, Fig. S2B, show two sharp anomalies at T N ≈ 26 K with zigzag antiferromagnetic order and T F ≈ 15 K with a spin reorientation transition in the ab plane.A weak anomaly of the χ(T) and C P (T)/T curves can be also seen around T * ≈ 7 K, which may be related to the modulation of magnetic domain.[3][4][5] Due to the magnetic-exchange anisotropy in NCTO, the magnetic order state exhibits a distinct magnetic anisotropic behavior along the different crystal-axis directions,   the C mag /T curves.

S3. Magnetic torque
Figure S4 shows the angle-dependent magnetic torque at 2 K.The backlash of the mechanical rotation is about 5 degree, marked by the blue horizontal error bars in Fig. S4.Therefore, we observe that the angle-dependent magnetic torque curves will have hysteresis obtained by counterclockwise and clockwise rotations.Surprisingly, a larger direction deviation of more than 5 degrees is seen between about 5 T and 7 T, which cannot be attributed to the error caused by the instrument.We note that an obvious hysteresis can be seen in the field dependence of magnetic torque and magnetization, Fig. S5, indicating the first-order phase transition around B C1 that is just between 5 T and 7 T. Hence, such large angular deviation of Fig. S4 between about 5 T and 7 T should be related to the first-order phase transition around B C1 .
Actually, the 'X phase' is not "another disordered phase".Instead, according to our experimental data and theoretical calculation, we identify it with a distinct phase with coexisting zigzag order phase and Z 2 topological order.The Z 2 topological order is reflected by the fact that there is a four-fold ground state degeneracy on a torus, indicating the deconfined Z 2 gauge fluctuations.
Firstly, the transition from the ordered phase to the intermediate disordered phase is of first order (which is consistent with our theory).The evidence for this first order phase transition is provided in the supplemental material, where obvious hysteresis structure is observed in the torque data.In Fig. S4, the angular dependence of the torque curves shows a big hysteresis with field around 6 T when the sample (or equivalently the field direction) is rotating clockwise or counterclockwise.Furthermore, in Fig. S5, the field-dependent torque curves (with fixed field intensity) also have a hysteresis structure around 6T (with fixed angle) when the field strength is increasing or decreasing.
Secondly, the angular dependence of the torque shows the same (the phases of the 2-period Fourier components are the same, see Fig. S6) 2-fold symmetric structure for field intensities B=3T, 6T, 7T, indicating the existence of zigzag order up to 7T.On the other hand, when the field strength exceeds 6T, the torque data eventually establish a 6-fold symmetry, indicating the appearance of disordered state.Especially, between B C1 =6T and B C2 =7.5T, the angular dependence of the torque contains both 2-fold and 6-fold symmetric components.
Thirdly, the first-order nature of the transition between the magnetically ordered phase and the disordered phase had also been observed in previous studies.
[ 5, 6]   From the above three facts, the most reasonable possibility is that between B C1 and B C2 is

E x p . f o r i n c r e a s i n g f i e l d E x p . f o r d e c r e a s i n g f i e l d F i t t i n g F i t t i n g
M ( m B / f .u .)   a phase with coexisting zigzag order and strong QSL fluctuations.Indeed, this phase is very interesting and to further confirm the coexistence of the two types of orders in this field region, more experimental explorations are needed.But this is not the focus of the present study.Our main interest is focused on the field-induced QSL phase between B C2 and B C3 .

E x p . f o r i n c r e a s i n g f i e l d E x p . f o r d e c r e a s i n g f i e l d F i t t i n g F i t t i n
In order to draw the magnetic phase diagram in the Fig. 1(c

S4. Inelastic neutron scattering
In the main text, the data of Fig. 1D-E

A. For the Ground States
The VMC approach is based on the spinon representation by introducing two species of fermions The spin operators are written in quadratic forms of the fermionic spinons , where m ≡ x, y, z, and σ m is Pauli matrix.It is convenient to introduce the matrix operator T such that the spin operators can also be written as Since the spin operators are invariant under a local SU(2) transformation ψ i → ψ i W i , the fermionic spinon representation has an SU(2) gauge 'symmetry'.The spin interactions are rewritten in terms of interacting fermionic operators and are further decoupled into a non-interacting mean-field Hamiltonian.In a quantum spin liquid (QSL) phase, all the physical symmetries are preserved in the ground state, and the SU(2) gauge symmetry may reduce to U(1) or Z 2 .For instance, the gauge symmetry of the Kitaev spin liquid is Z 2 .The mean-field Hamiltonian of a QSL phase does not necessarily preserve the physical symmetries, but should be invariant under an extended group called the 'projective symmetry group' (PSG).[9]   We adopt the same PSG of the Kitaev spin liquid, [10, 11] and the corresponding QSL mean field Hamiltonian of the K-J 1 -Γ-Γ ′ -J 3 model reads where ρ a , ρ c , ρ d , and ρ f come from the nearest neighbor interactions, t 3 and ∆ 3 come from the J 3 Heisenberg exchange interactions, σ m stands for the spin operation, τ m denotes the gauge operation, and λ λ λ is imposed to ensure the SU(2) gauge invariance (the λ z component is the Lagrangian multiplier for the particle number constraint).
Furthermore, to describe the long-range magnetic order, we introduce a background field M M M i to enforce the symmetry-breaking.The ordering pattern is contained by the single-Q Q Q ansatz [12] with where Q Q Q is the ordering momentum, ê e e x,y,z are the local spin axes, and ϕ is the canting angle.The ordering momentum Q Q Q is adopted either from the classical ground state.For a given Q Q Q, the local axes ê e e x,y,z are fixed as they are in the classical state, while M and ϕ are treated as variational parameters.
Hence, the complete trial mean-field Hamiltonian for the K-J 1 -Γ-Γ ′ -J 3 model in an external magnetic field reads ⟩ is computed using Monte Carlo sampling, and the optimal parameters R R R are determined by minimizing the energy E(R R R).
In the variational process, we first choose a classical metastable configuration with Q Q Q, then optimize the energy of the projected state.By comparing the optimal energies of different trial Q Q Qs, we obtain the approximate ground state.When the variational parameter M M M i is zero, the magnetic order disappear, and system goes into QSL phase.

Magnitude dependence of the magnetic field
The phase diagram of the K-J 1 -Γ-Γ ′ -J 3 model with K = 1.408meV,J 1 /|K| = −1.09,Γ/|K| = −0.94,Γ ′ /|K| = 0.63, J 3 /|K| = 0.94 under applied in-plane magnetic field applied along the a * = 1 √ 2 (x − ŷ)-direction is shown in Fig. S10(A), where four phases appear.The first phase is a classical phase with long-ranged zigzag order, the third phase and the fourth phase are the QSL phase and the trivial polarized phase, respectively.The second phase is an interesting phase with coexisting zigzag type magnetic order and Z 2 topological order (the QSL phase is also topologically ordered).The Z 2 topological order is reflected by the four-fold ground state degeneracy on a torus.
Remarks: The adopted interaction parameters in this work are proportional to the tx+ set of parameter in Ref. [13] , namely, J 1 =-3.5meV,K =3.2meV, Γ =-3.0meV,Γ ′ =2meV, excepted that the J 3 is enlarged to 3meV to ensure that the spin wave dispersion in the disordered phase is of concave shape at the Γ point (see section S6).Furthermore, to compare the excitation spectrum (see section S5 B) with experiment, we then divide all the parameters by a factor 2.5.Comparing the VMC result with the phase diagram of Fig. 1(a) in the main text, the theoretically obtained critical magnetic fields in Fig. S10, namelyB C2 = 0.47|K|/µ B ≈ 11.4T and B C3 = 0.77|K|/µ B ≈ 18.7T are slightly larger than the experimental values.This is possibly caused by the simplified K-J 1 -Γ-Γ ′ -J 3 effective model where only two-body interactions are considered.The ignored multi-spin interactions may suppress the zigzag order and reduce the critical magnetic field.
Noticing that the partially polarized QSL phase and the polarized trivial phase have the same symmetry.The difference is that the QSL phase is Z 2 deconfined such that the spions and gauge fluxes are elementary excitations and the polarized trivial phase is Z 2 confined such that the elementary excitations are magnons.The deconfinement of the Z 2 gauge field can also be reflected from the degeneracy of the ground states on a torus.If the Z 2 gauge field is deconfined, then Z 2 fluxes are allowed low energy excitations.Consequently, there will be four-fold degenerate ground states on a torus (given that the Chern number of the mean field ground state is zero).On the other hand, if the Z 2 gauge field is confined, then there will be a single ground state on a torus.
The four-fold degeneracy in the Z 2 deconfined QSL phase can be understood as the following.
Notice that the interactions are short-ranged, the energy density is a local quantity, but the global fluxes in the two holes are non-local.Hence in the thermodynamic limit the inserting of π-flux in one of the holes does not change the energy density of the system.The existence or absence of Z 2 fluxes through the two holes results in four-fold degenerate ground states on the torus.To judge the GSD, we calculate the fidelity matrix of the above four states with ρ αβ = ⟨ψ α |ψ β ⟩.
The four eigenvalues of the matrix ρ reflects the GSD.If all of the eigenvalues are of order 1, as shown in Fig. S11A for the partially polarized QSL phase (the fact that last eigenvalue is quit small may owe to the small Z 2 -flux excitation gap and the finite size effect), then the four states are linearly independent and consequently the GSD is four.The four-fold degeneracy of the ground states on a torus indicates the deconfinement of the Z 2 gauge field.On the other hand, if one of the eigenvalue is close to 4 and the other three are close to 0, as shown in Fig. S11B for the polarized trivial phase, then the GSD is one which indicates the Z 2 confinement.
In the partially polarized QSL phase, there is no symmetry breaking.Therefore, when changing the angle between B B B and a * , it is expected that the angle dependence of the off-diagonal magnetization should respect all the symmetries of the Honeycomb lattice.The numerical result is shown in Fig. S12, which indicates a 6-fold periodicity and is consistent with the expectation.x, y, z is defined by where n B (ω) is the Bose-Einstein distribution and χ ′′ mn (k k k, ω) = lim δ→0 Im χ mn (k k k, ω + iδ) is the imaginary part of χ mn (k k k, ω) which is the analytic continuitioin of the finite temperature susceptibility χ mn (k k k, iω).As an example, we illustrate the calculation of χ +− (k k k, iω) in the following ⟨T τ c α † ↑,q q q (τ )c α ↓,q q q+k k k (τ )c β † ↓,p p p+k k k (0)c β ↑,p p p (0)⟩e iωτ dτ where S ± = S x ±iS y and the superscript α, β sums over the A,B sub-lattices, T τ means time order and ⟨⟩ stands for thermal average.If we abbreviate the spin, sub-lattice and particle-hole indices of C-fermions (here ↑,q q q (τ )c α ↓,q q q+k k k (τ )c β † ↓,p p p+k k k β ↑,p p p (0)⟩ = ⟨T τ c a † q q q (τ )c b q q q+k k k (τ )c c † p p p+k k k (0)c d p p p (0)⟩, using Wick theorem, we have, ⟨T τ c a † q q q (τ )c b q q q+k k k (τ )c c † p p p+k k k (0)c d p p p (0)⟩e iωτ dτ = 1 V β 0 q q qp p p,{abcd} −⟨T τ c d p p p (0)c a † q q q (τ )⟩⟨T τ c b q q q+k k k (τ )c c † p p p+k k k (0)⟩ + ⟨T τ c c+4 p p p+k k k (0)c a † q q q (τ )⟩⟨T τ c b q q q+k k k (τ )c d+4, † p p p (0)⟩ e iωτ dτ = 1 V β 0 q q q,{abcd} −⟨T τ c d q q q (0)c a † q q q (τ )⟩⟨T τ c b q q q+k k k (τ )c c † q q q+k k k (0)⟩ + ⟨T τ c c+4 −q q q (0)c a † q q q (τ )⟩⟨T τ c b q q q+k k k (τ )c d+4, † −(q q q+k k k) (0)⟩ e iωτ dτ At mean field lever, the C fermions are diagonalized into the Bogoliubov particles, and the dynamic structure factors can be evaluated numerically.The complete dynamic structure factor for a system with 30 × 30 × 2 = 1800 sites is shown in Fig. S13A, where the spinon continuum can be clearly seen in the range ω < 3meV.

VMC results of the DSF
Now we try to construct the low-energy excitations via Gutzwiller projected spinon excited states |k k k; (k k k − q q q) m , (q q q) n ⟩ = P G f † k k k−q q q,m f † q q q,n |Ψ mf ⟩, where f k k k−q q q,m and f q q q,n are eigen particles of the mean-field Hamiltonian.Here m and n are the band indices, k k k and q q q are the lattice momentum according to the translation operators.
The above two-spinon excitations form a continuum and do not correctly describe the lowenergy excitations of the system.The failure of the state (S3) in describing the low-energy excitations is owing to the strong gauge interactions between the spinons which are not correctly addressed.To partially solve this problem, we diagonalize the original Hamiltonian in the subspace spanned by the two-spinon excitation continuum.Specifically, we calculate the matrix elements of the matrix H (k k k) with where k k k is the total momentum of the state with a pair of spinon excitations.
Since the states in the two-spinon continuum are not orthogonal, we need to calculate the metric matrix g formed by the overlap of the states, Hence, the eigenproblem of H (k k k) should be calculated by where the eigenvalues ϵ 1 , ϵ 2 , ... are the 'renormalized' energy of the excitations.
The renormalized eigenfunction |k k k⟩ rn is given by |k k k⟩ rn = q q q∈BZ F (qmn)|k k k; (k k k − q q q) m , (q q q) n ⟩, where F (qmn) is the eigenvector of the matrix H (k k k).
The more reliable DSF can be calculated using the renormalized energy and eigenstates.For instance, the zz-component of the DSF is expressed as S zz (q q q, ω) = n Ψ q q q n |S z q q q |Ψ 0 2 δ (ω − E q q q n + E 0 ) where |Ψ 0 ⟩ is the variational ground state with energy E 0 and |Ψ q q q n ⟩ is the n-th 'renormalized' twospinon excited state with momentum q q q and energy E q q q n .We note that S z q q q = 1 √ L r r r exp[iq q q • r r r]S z r r r is the Fourier-transformed spin operator for the component z, with L being the number of sites in the system.
However, due to the 'renormalizing' process, the computational complexity increases dramatically with system size.The data shown in Fig. S13B is the result for a system with size 6 × 6 × 2 = 72 sites.In Fig. S13B, the magnon modes as bound states of spinons can be clearly seen.However, due to finite size effect, the excitation gap of the magnons and the energy of the spinon continuum are exaggerated.
Noticing that the spinons are deconfined in the low-energy limit, the mean field results are qualitatively correct (but the bound states cannot be obtained at mean field level).Therefore, in the main text, we have combined the mean field result (with the spinon continuum) and the VMC result (with the magnon modes) to estimate the dynamic structure factor at large size limit.

S6. Linear spin wave
A. Ordered phase In the zigzag phase, the magnetic unit cell contains 4 sites, as shown in Fig. S14A.To obtain the linear spin wave (LSW) dispersion, we firstly adopt a new spin frame such that the new z-axis are parallel to the zigzag order, namely where R α is a SO(3) matrix with α =R,B stands for the red/blue sublattice index.Notice that the new z-axes in the red sublatice and the blue sublattice are opposite to each other.Thus the original K-J 1 -Γ-Γ ′ -J 3 model is transformed into a new one under the above rotations, where (with H ij the original interaction matrix on the ij-bond), B B B ′ i = R i B B B with R i = R R if i belong to the red sublattice and R i = R B for the blue sublattice, and g = 4.13. [5] In the new spin frame, we adopt the Holstein-Primakoff transformation S ′x i ∼ 1 2 (a † i +a i ), S ′y i ∼ i 2 (a † i − a i ), S ′z i = 1 2 − a † i a i on the A-sublattice and S where being the magnon creation operator on A (B) and red/blue sublattices.Diagonalizing the above Hamiltonian using the Bosonic Bogoliubov transformation, we obtain the magnon excitations spectrum at zero field, as shown in Fig. S14B.

B. Disordered phases
In the disordered phases, the spins are (partially) polarized due to the magnetic field.One can also calculate the spin wave spectrum assuming that the fully polarized state is the classical ground state.To this end, we perform an uniform orthogonal transformation to rotate the new z-axis to

S7. Random phase approximations
Using the determined model parameters, we extend our analysis to include spin fluctuations beyond the mean-field level by employing the random phase approximation (RPA).[14, 15] The RPA provides a means to incorporate the dynamical behavior of the interacting system by applying corrections to the bare susceptibility χ(0).The expression for the RPA-corrected susceptibility is given by χ(q, iω n ) = [I + V (q)χ (0) (q, iω n )] −1 χ (0) (q, iω n ) (S8) where I denotes the identity matrix and the kernel V(q) is defined as: In this equation, U is the strength of a phenomenological Hubbard term introduced for penalty of double occupancy of spinons, and J(q) denotes the Fourier transformation of the exchange interactions.At zero temperature, the bare spin susceptibility χ(0) is given by: χ (0)i,j µ,ν (q, ω) = meV and E ≈ 3 meV.These results agree well with both experimental data and VMC results.Our results also reveal a dispersive resonance mode centered at Γ point between 4 and 6 meV, a feature also captured by experimental data, but with strong damping.

FIG. 1 :
FIG. 1: Structure, magnetic properties and temperature-field phase diagram of NCTO.(A) Threedimensional stacking of the Co honeycomb layers.Honeycomb network shows the Co-Co bonds (red/blue/green) and edge-shared CoO 6 octahedra (the black, grey, cyan and golden spheres represent Co1, Co2, Te and O atoms, respectively).The Co1 and Co2 honeycomb layers present the ABAB-type layer stacking along the c axis.The dark red refers to the direction of moments that are in the ab plane and parallel to b-axis, indicating a zigzag AFM ground state. 42,43(B) A honeycomb network with three selected adjacent edge-shared CoO 6 octahedra.In the P6 3 22 structure, Co ions form a perfect honeycomb lattice with an equal 92.17 • Co-O-Co bond angle and the nearest-neighbor Co-Co bond length d Co−Co = 3.05 Å. (C) Temperature-field phase diagram with B parallel to Co-Co bonds.The phase boundaries are deduced from the temperature-dependent magnetic specific heat C mag /T with B // a* (Supplementary Figure S2E) and field-dependent differential magnetic torque 1 B

FIG. 2 :
FIG. 2: Symmetry evolution of magnetic torque with temperature and magnetic field.(A)-(D) Polar plots of magnetic torque τ (θ) at different temperatures with the fixed magnetic fields.(E)-(G) Polar plots of magnetic torque τ (θ) at different magnetic fields with the fixed temperatures.(H) Temperature-dependent amplitude of C6 and C2 symmetry obtained by Fourier transform of angle-dependent magnetic torque at different fields.

FIG. 3 :
FIG. 3: Spin-excitation spectra using fixed incident energy E i = 18 meV with B ∥ a*.Constant-energy scattering at 0 T (A) and 8 T (B), respectively, integrated over L = [-2.5,2.5] and E = [1.5, 2.5] meV, projected on the reciprocal honeycomb plane defined by the perpendicular directions [H, H,0] and [K, -K, 0].The white dashed lines represent the Brillouin zone boundaries.(C) and (D) Spin-excitation spectra along high symmetry momentum directions Γ-X-K-M-Γ-M 1 at 5 K for zero field and 2 K for 8 T, respectively.The color bar indicates scattering intensity with arbitrary unit in linear scale.The dark red dashed line indicates that an intense 'Λ' shape spinon continuum appears at 4 meV ˜8 meV in Figure (D).(E) and (F) Calculated dynamic structure factor for the zigzag AFM order and field-induced QSL with partially polarized spins, respectively.The black regions lack detector coverage.The dark red dashed line shows the same 'Λ' shape spinon continuum, which is compared to the experiment at 8T.The data were collected using the SEQUOIA chopper spectrometer at the Spallation Neutron Source (SNS)..
FIG. S1: of the crystal structure of Na 2 Co 2 TeO 6 .(A) Definitions of the crystallographic axes (a, b, c) and the Ising spin axes (x, y, z).(B) Honeycomb network showing the Co-Co bonds (red/blue/green) and CoO 6 octahedras (the black, cyan and golden spheres represent Co, Te and O atoms, respectively).The black arrows are the two principal reciprocal vectors a * and b * in first Brillouin zone.Magnetic field B is applied within the honeycomb plane.θ is the angle between B and a * and with respect to the real-space orientation of the Co-Co bonds.The black and red fonts represent the coordinate vector in reciprocal space and spin-axes coordinate, respectively.Directions that are equivalent to [1,0,0] and [1,1,0] correspond to angles θ = n × 60 • and θ = n × 60 • + 30 • , respectively (where n is an integer) in reciprocal space coordinate.

Fig. S2A .
Fig.S2A.However, the χ(T) curves also present a strongly anisotropic behavior above T N in the ab plane and c-axis, which is mainly attributed to the strongly anisotropic g factor originating from the low-spin state of Co 2+ .Hence, the spin-1/2 model is an effective, low-energy description for the Co 2+ ions in the temeprature range considered in this work, which is consistent with our previous report.

[ 5 ]
FIG. S2: Macroscopic physical properties of Na 2 Co 2 TeO 6 .(A) Temperature dependence of susceptibility χ(T) in NCTO measured at 0.01 T along the different crystal-axis directions.(B) Temperature dependence of specific heat under 0 T for NCTO and Na 2 Zn 2 TeO 6 .(C) The isothermal magnetization M(B) with the applied magnetic field B // a * -axis at selected temperatures.(DD) The differential isothermal magnetization as functions of fields dM/dB vs. B, originating from the data of (C).The solid lines represent the fitted curves by polynomial in (C) and (D).The critical fields B C1 and B C2 are determined by the dM/dB curves, which are comparable with the critical fields in the main text.(E) Temperature dependence of magnetic specific heat in NCTO under selected fields.

Figure
Figure S2(E) shows the magnetic specific heat C mag /T curves with B // a * -axis, subtracting phonon contributions measured on an isostructural Na 2 Zn 2 TeO 6 reference crystal that is shown in the Fig. S2B.The critical temperatures in Fig. 1(A) in the main text is marked by the anomalies in

FIG. S4 :
FIG. S4: Angle-dependent magnetic torque without subtracting background at 2 K with different fields.The blue horizontal error bars indicate the instrumental backlash (∼ 5 degrees) when reversing rotation.The black solid lines represent the increasing angles from 0 to 360 degrees.The red solid lines represent the decreasing angles from 360 to 0 degrees.
FIG. S5: Field dependence of magnetic torque and magnetization.Field dependence of magnetic torque (A) and differential magnetic torque (B) in NCTO measured at 3 K with magnetic field 10.6 degrees away from a * -aixs.Field dependence of magnetization (C) and differential magnetization (D) in NCTO measured at 2 K along the field B // a * -axis.

FIG. S7 :
FIG. S7: Isothermal magnetic torque with field close to a * -axis.Field dependence of magnetic torque (A) and differential magnetic torque (B) in NCTO measured at different temperatures along the angle θ = 10.6°.
) of the main text, we measure the field dependence of magnetic torque along the angle θ = 10.6 • at different temperatures, Fig. S7.The superimposed color map of Fig. 1(c) is from the differential 1 B dτ dB of Fig. S5.
Fig. S8B.The software suite DAVE was used to visualize the time-of-flight neutron scattering data and the space group P6 3 22 with a = b = 5.227 Å, c = 11.2231Åwas also adopted to analyze these INS results.Figures S13A-B show the rod-like dependence on the out-of-plane momentum component L, demonstrating the inter-layer interaction is relatively weak, NCTO could be regarded as a 2D Co 2+ -based honeycomb magnet.

[ 8 ]
FIG. S9: Rod-like dependence on the out-of-plane momentum component L. (A) and (B) Constantenergy scattering at 0 T and 8 T, respectively, integrated over H = [-0.1,0.1] and E = [1.5, 2.5] meV, projected on the reciprocal honeycomb plane defined by the perpendicular directions [0, 0, L] and [K, -K, 0].The color bar indicates scattering intensity in arbitrary unit in linear scale.The black regions lack detector coverage.
where B B B is the effective Zeeman field due to the external magnetic field B B B, but B B B is not necessarily equal to B B B. Then we perform Gutzwiller projection to the mean-field ground state |Ψ mf (R R R)⟩ to enforce the particle number constraint.The projected states |Ψ(R R R)⟩ = P G |Ψ mf (R R R)⟩ provide a series of trial wave functions depending on the choice of the mean-field Hamiltonian H mf (R R R), where P G denotes a Gutzwiller projection and R R R are treated as variational parameters.The energy of the trial state

1 √ 2
FIG. S10: (A) Calculated phase diagram with varying field strength and (B) the spinon dispersion in the QSL phase.The field is applied along the a * = 1 √ 2 (x − ŷ)-direction, and the in-plane Landé g factor is adopted as g = 4.13.
Since the inserting of a π flux in one of the holes is equivalent to change the boundary condition of the mean field Hamiltonian in the corresponding direction from the periodic one to anti-periodic one (vice versa).If we note periodic and anti-periodic boundary conditions as + and − respectively, then there are four different combinations of boundary conditions for the xand ydirections, namely (+, +), (+, −), (−, +), (−, −).The four boundary conditions results in four projected states, namely |ψ α ⟩ = P G |α⟩ mf where |α⟩ mf denotes the mean field ground state with boundary condition α.If the ground state degeneracy(GSD) of the spin model on the torus is four, then |ψ (+,+) ⟩, |ψ (+,−) ⟩, |ψ (−,+) ⟩, |ψ (−,−) ⟩ should be linearly independent.Otherwise, if the GSD is 1, then the above four states are essentially the same.

2 .
FIG. S11:Eigenvalues of the fidelity matrix for the ground states on a torus.(A) In the QSL phase, all the eigenvalues are approaching 1 with the increasing size, indicating that in thermodynamic limit the GSD is 4 and the Z 2 gauge field is deconfined.(B) In the polarized phase, the GSD is one, indicating Z 2 confinement.
FIG.S12: Calculated angle-dependent off-diagonal magnetization with fixed field strength.The angle-dependence exhibits a 6-fold periodicity.

2
− b † i b i on the B-sublattice, where the bosons satisfy the usual commutation relations [a i , a † j ] = δ ij , [b i , b † j ] = δ ij and [a i , b † j ] = 0. Substituting the above formulas into the rotated Hamiltonian (S5) and keeping the quadratic terms, we obtain the following Hamiltonian on the Fourier bases, FIG. S14: Linear spin wave dispersion.(A) The Zigzag order on honeycomb lattice.(B) Linear spin wave of the Zigzag phase at zero field.(C) Linear spin wave of the disordered phase at µ B B/|K| = 0.56.(D) Linear spin wave of the polarized phase at µ B B/|K| = 0.95.
FIG. S15: Calculated dynamic structure factors.The dynamic structure factors obtained by RPA with 30 × 30 × 2 sites at U = 6.