Introduction

Magnets with significant spin–orbital couplings (SOCs) have become a new playground for quantum magnetism in recent years thanks to their potential of hosting novel quantum phases of matter1. A prominent example is the Kitaev model, an exactly solvable spin model that features a topological quantum spin liquid (QSL) ground state2,3,4,5. Microscopically, such a model could emerge from magnetic insulators with competing, spin anisotropic exchange interactions3. The essential prerequisites are that (a) the magnetic ions have spin–orbital entangled, Kramers degenerate ground states, and (b) they are arranged on a suitable lattice, with the two-dimensional (2D) honeycomb lattice being the simplest example2,3,4,5. So far, the search for the material incarnations of the Kitaev model has been focused on 4d/5d transition mental-based systems due to their relatively strong SOCs4,5. The examples include H3LiIr2O6, α-Li2IrO3, α-Na2IrO3, and α-RuCl34,5,6,7,8,9,10,11,12,13,14,15. H3LiIr2O6 was initially thought to be a QSL, but later studies argue for a random singlet state resulted from the quenched disorder induced by the mobile protons6,7. As for α-Li2IrO3, α-Na2IrO3, and α-RuCl3, all of them magnetically order due to the non-Kitaev interactions8,9,10,11,12,13,14,15. To describe these materials, the Kitaev model has been extended to a generalized Heisenberg–Kitaev (H–K) model with five symmetry-allowed terms, Kitaev term K, off-diagonal symmetric exchange term Γ and Γ′, nearest-neighbor (NN) Heisenberg coupling J, and the third NN Heisenberg coupling J38,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24. α-RuCl3 is unique among these materials in that its effective Hamiltonian features a dominant Kitaev term12,15,19,20,21. Remarkably, applying an in-plane magnetic field around 7 T greatly suppresses its magnetic order and induces a potential QSL state23,25,26,27,28,29,30,31,32,33,34.

More recently, the theoretical studies suggest that Kitaev physics could also be found in honeycomb magnets made of cobalt, a 3d transition metal35,36,37. In 3d7 Co (t2g5eg2) compounds, the spin-active eg electrons not only generate new spin–orbit exchange channels of t2g–eg and eg–eg in addition to the t2g–t2g channel operating in d5 systems with t2g-only electrons but also produce the Kitaev interaction almost entirely from the t2g–eg process37. Moreover, the t2g–eg and eg–eg contributions to Γ and Γ′ are of opposite signs and largely cancel each other37. This makes the cobaltates good candidates for realizing the Kitaev model35,36,37. A pertinent material is the honeycomb magnet Na2Co2TeO6 (NCTO)35,36,37. The other example we are aware of is BaCo2(AsO4)2, which also shows a similar field-induced disordered state at a low magnetic field38.

As we shall present below, our studies show that NCTO fulfills all the prerequisites for Kitaev physics despite its seeming differences from the preceding 4d/5d compounds. Furthermore, we observe a QSL-like spin disordered state in an applied magnetic field B ⊥ b-axis (parallel to Co-Co bond) and 7.5 T < B < 10.5 T. This result suggests NCTO could be a novel honeycomb magnet that exhibits a field-induced disordered state and may broaden our horizon in the quest for Kitaev materials.

Results

Specific heat and magnetic susceptibility

The magnetic specific heat at 0 T measured on the polycrystalline sample shows a sharp peak at TN = 25 K, and its derivative shows another two anomalies at TF = 15 K and T* =7 K, Fig. 1a. Their values are consistent with the previous reports39,40,41,42,43,44. With an increasing magnetic field, the sharp peak at TN shifts to lower temperatures and becomes broader. Eventually, it is indiscernible above 7.5 T44, which indicates that the magnetic ordering at TN has been suppressed, and the system enters a magnetically disordered state. Supplementary Fig. 2 presents the heat capacity as a function of temperature at various fields, and the coincident heat capacity suggests no missing entropy above 50 K. The magnetic entropy Smag is obtained after subtracting the lattice contribution calculated by the Debye–Einstein (DE) model using Supplementary S(1) (Supplementary Fig. 2). Integrating Cmag/T over temperature from 100 Kelvin down to 2 Kelvin yields the residual magnetic entropy Smag(0 T) = 3.05 J/(mol K), which is about 27% of the theoretical total magnetic entropy, Fig. 1b. Remarkably, the residual entropy increases with increasing fields in NCTO. This unusual behavior may be contrasted with other magnets that are known to possess finite residual magnetic entropy, where the magnetic field tends to reduce the residual entropy rather than enhance it45,46,47,48.

Fig. 1: Heat capacity of Na2Co2TeO6.
figure 1

a Temperature dependence of the magnetic specific heat of NCTO for different magnetic fields up to 14 T. The inset shows the differential magnetic specific heat dCmag/dT, and three transitions are observed as TN, TF, and T*, respectively. b The magnetic entropy Smag(T) at 0, 9, and 12 T. The olive dashed line refers to Smag(T → ∞) calculated with effective spin Jeff = 1/2 for Co2+. In order to better observe the residual magnetic entropy, the Smag(T) curves are shifted vertically so that their maxima coincide with the Smag(T → ∞). The black, red, and blue dashed line refer to the possible residual magnetic entropy below 2 K under different magnetic field. Compared with the total magnetic entropy Smag(T → ∞), the percentage of residual magnetic entropy in different magnetic fields is shown in the figure. c A zoom-in of magnetic entropy below TN at various magnetic fields. Note the data are not shifted vertically as opposed to (b). The inset shows the field dependence of the Smag at 17 K.

Full size image

Another noteworthy feature is that the field dependence of the Smag at 17 K (the insert of Fig. 1c), reaches a maximum near 9 T. This non-monotonic behavior of Smag with increasing fields mirrors that of the magnetic specific heat Cmag/T (Fig. 1a), which suggests the energy gap of the magnetic excitations closes at first and then reopens with increasing fields. It is worth mentioning that this field-dependent non-monotonic behavior was also observed in the α-RuCl3 as considered a signal of entering the filed-induced non-Abelian QSL state25,27,29.

Figure 2 shows the DC magnetic susceptibility χ(T) measured by the zero-field cooling process on a single crystalline sample with B ⊥ b-axis in the ab plane. In low fields, three anomalies at TN, TF, and T* are observed as shown in Fig. 2a, which mirror the specific heat anomalies. With increasing fields, the peak with the antiferromagnetic (AFM) characteristics at TN shifts to lower temperatures and becomes flattened, and the other two anomalies become weaker, Fig. 2b. With B > 7.5 T, all anomalies are indiscernible. This is consistent with the disappearance of the specific heat peak with B > 7.5 T and suggests that the system enters a disordered phase.

Fig. 2: Magnetic susceptibility of Na2Co2TeO6 for B ⊥ b-axis.
figure 2

a, b Temperature dependence of susceptibility χ(T) in NCTO measured at various magnetic fields. The inset of Fig. 2b shows the honeycomb lattices of Co viewed along the c-axis. The cartoon shows the moments to be in the basal plane and parallel to the b-axis40, 41, and the TeO6 octahedra sit at the center of each honeycomb unit. The first NN ferromagnetic (FM) interaction J originated from a collaboration between Co–O–Co (red bond) superexchange interactions and Co–Co (blue bond) direct exchange interactions. The second NN J2 has multiple superexchange interaction pathways, which may lead to a wide range of variations, mainly determined by J and J3. The third NN AFM interaction J3 arises from the existence of the Te atom in the center of the honeycomb lattice, which leads to the unique superexchange interaction pathway Co–O–Te–O–Co (golden bond). The olive arrow refers to the direction of applied magnetic field B in the magnetization. c The isothermal magnetization M(B) with the applied magnetic field B ⊥ b, B // b, and B // c-axis at 0.5 K, respectively. The dashed line indicates the Van-Vleck paramagnetic background. d The differential isothermal magnetization as functions of fields dM/dB vs. B at different temperatures with B ⊥ b.

Full size image

The magnetization (M) measured at 0.5 K with B ⊥ b-axis shows an anomaly around 6 T (Fig. 2c), which leads to a peak on the dM/dB curve (Fig. 2d). The peak shifts to a lower field and becomes weaker with increasing temperature. Above 7.5 T, the dM/dB data show two opposing temperature-dependent behaviors: the intensity of the dM/dB decreases with increasing temperature below 10.5 T, whereas it increases with a temperature above 10.5 T. The dM/dB curves at different temperatures all cross approximately at 10.5 T. We also note that for temperatures below 6 K and 7.5 T < B < 10.5 T, when the system is in the disordered state as suggested by the specific heat and susceptibility data, the dM/dB curves qualitatively coincide with each other.

We interpret our magnetization data heuristically as follows18. Between 7.5 and 10.5 T, the system is in a spin-disordered state with short-range spin correlations. On the one hand, the short-range spin correlations increase the dM/dB intensity with increasing temperatures since these thermal fluctuations thermally activate the spins that can be flipped by the magnetic field. On the other hand, the increasing temperature thermalizes the already polarized spins to decrease the dM/dB intensity. Therefore, the temperature-independent dM/dB intensity below 6 K and 7.5 T < B < 10.5 T indicates there exist strong short-range spin correlations for such low temperatures in the spin disordered state. Above 6 K and 7.5 T < B < 10.5 T, the thermal fluctuations quickly thermalize the already polarized spins so that the dM/dB intensity decreases with increasing temperature. By contrast, the system crosses over to the polarized state above 10.5 T. Thus, the temperature-independent dM/dB intensity at 10.5 T reflects a characteristic field at which these two competing effects compensate each other. We, therefore, take 10.5 T as the crossover field from the correlated spin-disordered state to the spin-polarized state. The magnetization curve further suggests that the saturation field for the B ⊥ b-axis case is around BS = 12.5 T with saturation magnetization MS = 2.01μB/Co2+ obtained by subtracting off the Van-Vleck paramagnetic contribution49 (Fig. 2c).

High-field electron spin resonance (ESR)

As shown in Fig. 3, the high-field ESR data measured at 2 K exhibits a rich excitation spectrum, in which the modes A–C were observed in the ordered phase and the mode D was only detected with B > 6 T. The modes A–C can be described by conventional AFM resonance modes and single-magnon excitations at the Γ point33,50, the related intensities become weaker and disappear above TN with the increasing temperature (Supplementary Fig. 3). The inset of Fig. 3k shows the electron paramagnetic resonance (EPR) spectra measured at 50 K with 214 GHz. From the resonance fields obtained by Supplementary S(2), the calculated g-factors are gab = 4.13 and gc = 2.3, respectively. These values are also consistent with the magnetization data. The saturation magnetization MS = 2.01μB/Co2+ for B ⊥ b-axis, Fig. 2c, leads to a gJ = 4.02 for the pseudospin-1/2 case. Both the ESR and the magnetization data corroborate the effective spin-1/2 model for Co2 + ions. The higher energy crystal field levels of Co2+ are thermally inactive in the temperature range considered in this work. The strongly anisotropic g factors for Co2+ ions in the octahedral environment confirm the strong SOC and magnetocrystalline anisotropy, which can usually drive a magnetic insulator to open a spin-wave energy gap. The extracted resonance data from Fig. 3a–j are summarized in the frequency-field diagram shown in Fig. 3k, and the extrapolation of the frequency-field dependences of modes A–C to zero field reveal a gap ΔE ≈ 100 GHz ≈0.4 meV.

Fig. 3: High-field ESR of Na2Co2TeO6.
figure 3

a–j Frequency-dependence of ESR at 2 K; k ESR frequency-field diagram of NCTO at 2 K. The unit conversion with meV (1 meV ≈241.8 GHz) is shown on the right axis. The inset shows the EPR spectra measured at 50 K with 214 GHz. The red lines are the fitting line by Supplementary S(2). The olive squares and magenta stars are obtained from the EPR data at 50 K.

Full size image

Around 7–8 T, the AFM modes of A and B approach the EPR line with gab = 4.13, which again suggests a field-driven magnetically disordered state with paramagnetic-like behavior in the ab-plane. The field-frequency curves of the A and B modes intersects with the EPR line near 6 T, which may be heuristically interpreted as the field at which the spin gap closes. This interpretation is also consistent with the specific heat data, which suggests the gap closes near 6 T. The field-frequency curve of the C mode gradually approaches the EPR line with gc = 2.3, which indicates that its polarization is along the c-axis.

The ESR measurement reveals another mode that is not directly connected to the aforementioned AFM resonance modes, which we dub as D mode. It only appears when B > 6 T. Its excitation energy shows a linear-field dependence with a slope of ~0.15 meV/T, from which we deduce an effective g ≈ 2.6. This effective g-factor is between gab and gc. The D mode must be associated with a magnetic excitation that only exists or becomes visible in the high-field spin disordered phase. Comparing with the other three modes, the D mode is much weaker, and its linewidth is broader. Its microscopic origin is unclear at the moment; however, we note its close resemblance with the ESR data of α-RuCl3, where new modes with linear field-energy relationship emerge in the spin disordered state26,28,33.

Phase diagram

A phase diagram was constructed in Fig. 4 by combining the critical temperatures and fields obtained above. There are four regimes: (i) with B < 6 T, three magnetic transitions occur with decreasing temperatures. The transition at TN should be the one from paramagnetic to zigzag AFM ordering. The ones at TF and T* could be the adjustments of the canting moments of the zigzag order; (ii) For 6 T < B < 7.5 T, there is a coexistence of low-field magnetically ordered state and high-field spin disordered state; (iii) most interestingly, within 7.5 T < B < 10.5 T, the system enters a spin disordered state; (iv) with B > 10.5 T, the system begins to enter the polarized state and becomes fully saturated above 12.5 T.

Fig. 4: Phase diagram of Na2Co2TeO6.
figure 4

Temperature–magnetic field phase diagram for NCTO. The stars symbol the first-ordered phase boundary. The color background is used only as a simple guide.

Full size image

Spin-wave excitations

Figure 5a presents the inelastic neutron scattering (INS) measurement with incident energy Ei = 11.4 meV for NCTO at 3 K and two bands (from 0.4 to 2.9 meV and from 5.9 to 7.1 meV) are observed. The magnetic mode shows an apparent minimum near Q = 0.72 Å −1, which is close to the magnitude of the M point of the honeycomb reciprocal lattice as expected for a 2D magnetic system. Moreover, the concave shapes are observed at the scattering edges, similar to the magnon excitations observed in other honeycomb lattice magnets with zigzag AFM order of α-Na2IrO351 and α-RuCl312. As shown in Fig. 5c, d, the INS spectra show a gap of 0.4 meV at the M point. This gap is comparable in energy scale with the ESR gap, which corresponds to the excitation energy at the Γ point.

Fig. 5: Powder INS of Na2Co2TeO6.
figure 5

a Powder INS using 11.4 meV incident neutron for NCTO at 3 K on spectrometer HRC. The high-temperature data (95 K) was subtracted as background and the thermal effect was simulated by the Bose factor. b The calculated powder-averaged scattering including the magnetic form factor with K = Γ =0.125 meV, J = −2.325 meV, and J3 = 2.5 meV. The black and white dashed lines label M and K points at the first Brillouin zone with Q = 0.72 and 0.83 Å −1, respectively. c Constant-E cuts over the energy ranges centered at the locations labeled dashed lines in Fig. 5d. The solid blue line is constant-E cuts over the energy ranges centered at 0.66 meV from our LSWT. The absence of structured scattering below 0.4 meV confirms the gap in the magnetic excitation spectrum. d Powder inelastic neutron scattering using 5 Å (~3.27 meV) incident neutrons for NCTO at T = 3 K on spectrometer NEAT II, HZB, Deutschland.

Full size image

The diffraction refinement suggested that the magnetic moments were along the crystallographic b-axis with a possible small canting toward the c-axis40,41. Since the inter-layer interaction is relatively weak, NCTO could be treated as a single-layer 2D compound at a first approximation. The super–super-exchange Co–O–Te–O–Co pathway produces a significant third-neighbor exchange interaction J341. Meanwhile, the NN Co ions can interact with each other through two 90° Co–O–Co super-exchange paths or direct AFM exchange interactions35,41. The cancellation between the different hopping contributions from d–d and d–p–d orbits can weaken the ferromagnetic NN J52,53. Together, J and J3 can stabilize a zigzag magnetic order. The second neighbor J2 is likely to be relatively weak in that it tends to destabilize the zigzag order. However, as J1,2,3 is isotropic in the spin space, they cannot select the direction of the magnetic moments; this is achieved by the spin anisotropic interactions such as the K and Γ terms mentioned in the Introduction.

The linear spin-wave theory (LSWT) is performed to analyze the INS spectra with the following exchange Hamiltonians17,18,19,20,22,23

$${{\rm H}}=\mathop{\sum}\limits_{\langle i,j\rangle }[J{{{{{{\bf{S}}}}}}}_{i}\cdot {{{{{{\bf{S}}}}}}}_{j}+K{S}_{i}^{\gamma }{S}_{j}^{\gamma }+{\Gamma }({S}_{i}^{\alpha }{S}_{j}^{\beta }+{S}_{i}^{\beta }{S}_{j}^{\alpha })]+{J}_{3}\mathop{\sum}\limits_{\langle \langle \langle i,j\rangle \rangle \rangle }{{{{{{\bf{S}}}}}}}_{i}\cdot {{{{{{\bf{S}}}}}}}_{j}$$
(1)

where, <i, j > denotes NN sites, Si and Sj are effective spin-1/2 operators at sites i and j, respectively, α and β are perpendicular to the Kitaev spin axis γ. When the J2 and I (Ising exchange interaction) are close to zero (Supplementary 4.1 Symmetries and model), the zigzag AFM order will be more stable. Finally, the powder-averaged scattering numerical results are presented in Fig. 5b. With K = Γ =0.125 meV, J = −2.175 meV, and J3 = 2.5 meV, the calculated dispersion can reproduce qualitatively the experimental data. While the INS cannot access the excitation energy gap at the Γ point, the LSWT calculation suggests an energy gap on the same order of magnitude as the M point gap, in qualitative agreement with the ESR results.

Discussion

The characteristic behaviors for the field-induced spin disordered state in NCTO, such as the disappearance of the peaks observed on specific heat and χ(T), the field dependence of magnetic entropy with a maximum near 9 T below 17 K, and the existence of low energy excitations at low temperatures indicated by the dM/dB curves, are all similar to those observed for the field-induced disordered state above about 7 T in α-RuCl327,29,30,54,55,56. These behaviors have been believed to be evidence that α-RuCl3 enters a Kitaev QSL state under fields27,29,30. The ESR measurements of α-RuCl3 also show extra modes in the field-induced disordered state with a linearly increasing energy gap with increasing magnetic fields26,33, similar to that of the D mode observed for NCTO. For α-RuCl3, it has been suggested that when the magnetic field and gap become large enough, it can overcome the energy scale related to the residual Heisenberg interactions so that a QSL emerges23,26,28,32,33. While the exact nature of this field-induced disordered state in NCTO needs further studies to determine, the similarities between the disordered states of α-RuCl3 and NCTO suggest the possibility of the latter being a QSL state.

We note several recent reports on the magnetic excitations of NCTO with different interpretations57,58,59. Ref. 57 performs the INS of NCTO on the same facility as this work but at the higher incident neutron beam energy of 16.54 meV. The resulted INS spectra are qualitatively similar to this work. Moreover, we obtained additional low energy spectra of NCTO with the incident energy of 3.27 meV and achieved better instrument resolution by using another neutron spectrometer NEAT II, HZB, Deutschland. The result (Fig. 5d) shows an energy gap of about 0.4 meV at the M-point. It is noticed that Ref. 57 applies a different model on the INS spectra. Ref. 57 suggests a large J (−1.5 meV) and a large AFM K (3.3 meV), whereas our simulation leads to a large J (−2.325 meV) and a small AFM K (0.125 meV). Surprisingly, both could qualitatively reproduce the INS spectra despite different choices for the relative magnitude of the K term.

Meanwhile, ref. 58 suggested a small J (−0.1 meV) and a large FM K (−9 meV); ref. 59 compared two sets of exchange interactions, which are J (−0.2 meV), K (−7 meV) and J (−3.2 meV), K (2.7 meV), respectively. In ref. 59, it was pointed out that the calculated spectra with the two models, one with FM K and the other with AFM K, are almost indistinguishable. Therefore, at this stage, it is difficult to judge which set of exchange interactions is more accurate.

The main reason for the tension among these different fitting parameters, such as a large K term vs. a small K term, or FM vs. AFM K term, could be due to the fact that all these INS measurements were performed on polycrystalline samples. The powder average effect on the weak signals makes it difficult to determine the Hamiltonian parameters for the complex systems with competing interactions and frustration. Single-crystal neutron scattering measurement in the future is critical to clarify the magnetic structure and dynamics of NCTO while it could be challenging due to the small size of single crystals grown by the flux method44.

In summary, the most significant finding from our studies on NCTO presented here is a QSL-like disordered state induced under fields with 7.5 T < B < 10.5 T. Therefore, NCTO is a novel example of an effective spin-1/2 honeycomb magnet that hosts a field-induced spin disordered state. Its origin, in addition to the related spin structure and spin dynamics, calls for future experimental work on single crystals and theoretical studies.

Methods

Sample preparation and characterization

NCTO polycrystalline was prepared by a solid-state reaction method. At first, Na2CO3 (Alfa, 99.997%), Co3O4 (Alfa, 99.7%), and TeO2 (Alfa, 99.99%) were mixed in a stoichiometric molar ratio with 5% excess Na2CO3, and fully ground; then, the mixture was loaded in an alumina crucible and sintered at 850 °C in the air for 40 h. The high-quality single crystal was grown by the flux method. The polycrystalline sample of NCTO was mixed with the flux of Na2O and TeO2 in a molar ratio of 1:0.5:2 and gradually heated to 900 °C at 3 °C/min in the 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 to cool to room temperature. To confirm the structure and purity of the sample, powder X-ray diffraction (XRD) measurement was performed with a HUBER imaging-plate Guinier camera 670, using Cu Kα1 radiation (λ = 1.54051 Å). The XRD patterns were refined with the Rietveld method using the conventional refinement program FullProf (Supplementary Fig. 1). The magnetic properties were checked through measurements as a function of temperature (T) and magnetic field (B) using a vibrating sample magnetometer in the physical properties measurement system (PPMS, Quantum Design). Besides, magnetization measurements were also carried out using a high-sensitive Hall sensor magnetometer60,61,62 for the temperature range from 0.4 to 30 K. The heat capacity measurements were carried out using the relaxation time method in the PPMS.

High-field ESR

The high-field ESR measurements were performed by the high-field high-frequency electron magnetic resonance spectrometer with 25 T water-cooled resistive magnet (field sweep range: 0–25 T, frequency range: 50–690 GHz, and temperature range: 2–300 K).

Inelastic neutron scattering

The INS of powder NCTO was performed using the High-Resolution Chopper Spectrometer (HRC) at the J-PARC63. The HRC delivers high-resolution and relatively high-energy neutrons for a wide range of studies on materials dynamics. Spectra were collected at various temperatures by operating in high-flux mode (energy resolution of ≥2%Ei) with Ei = 11.4 meV. The INS of powder NCTO was also carried out on the recently upgraded cold-neutron direct-geometry time-of-flight spectrometer NEAT II, HZB, Deutschland64,65. Each sample was packed in aluminum cans filled with He exchange gas. Each scan was counted for around 6 h with the incident neutron energy Ei = 3 Å (~9 meV) and 5 Å (~3.27 meV).