Sharp spectroscopic fingerprints of disorder in an incompressible magnetic state

Domain: physics
Published: 2025-12-06
Editorial outcome: minor_revision

Abstract
Disorder significantly impacts the electronic properties of conducting quantum materials by inducing electron localization and thus altering the local density of states and electric transport. In insulating quantum magnetic materials, the effects of disorder are less understood and can drastically impact fluctuating spin states like quantum spin liquids. In the absence of transport tools, disorder is typically characterized using chemical methods or by semi-classical modeling of spin dynamics. This requires high magnetic fields that may not always be accessible. Here, we show that magnetization plateaus—incompressible states found in many quantum magnets—provide an exquisite platform to uncover small amounts of disorder, regardless of the origin of the plateau. Using optical magneto-spectroscopy on the Ising-Heisenberg triangular-lattice antiferromagnet K2Co(SeO3)2 exhibiting a 1/3 magnetization plateau, we identify sharp spectroscopic lines, the fine structure of which serves as a hallmark signature of disorder. Through analytical and numerical modeling, we show that these fingerprints not only enable us to quantify minute amounts of disorder but also reveal its nature—as dilute vacancies. Remarkably, this model explains all details of the thermomagnetic response of our system, including the existence of multiple plateaus. Our findings provide a new approach to identifying disorder in quantum magnets.

Introduction
Chemical disorder is a fundamental characteristic of solid-state materials and can manifest in various ways, such as impurities and dopants, site mixing and vacancies, bond disorder, and topological defects. In quantum materials, this disorder is imprinted at high temperatures during synthesis and becomes kinetically frozen, often remaining unnoticed until it influences putative emergent phenomena at low temperatures. However, disorder is not always detrimental in quantum systems. A notable example is the quantum Hall effect (QHE), where disorder plays a crucial role in the formation of quantized plateaus in the Hall resistance of devices and materials. This phenomenon has led to the establishment of a new international standard for electrical resistance1,2.

The role of disorder in quantum magnetic materials is less understood, yet it is crucial for a wide range of systems, particularly quantum spin liquids (QSLs). These systems are believed to host topological order and itinerant fractional excitations without exhibiting clear, conventional order-parameter-like signatures. The response, survival, and nature of potential QSLs in the presence of chemical disorder are not fully understood. Increasingly, researchers are recognizing that disorder can mimic the expected behavior of a QSL3–8. But a significant challenge lies in quantifying the level of disorder using chemical techniques and translating this information into a model magnetic Hamiltonian, H, whose emergent quantum many-body phases are not known in advance. A successful approach involves employing electron spin resonance9,10 or fitting inelastic neutron scattering and optical spectroscopy measurements under sufficiently high magnetic fields to polarize the system. In this semi-classical regime, it becomes possible to directly measure H11, including the level of chemical disorder12,13. However, practical limitations arise from the need for large magnetic fields. Furthermore, spin-anisotropic systems of current interest may only reach saturation asymptotically14.

The complex magnetization process of antiferromagnets may offer an alternative strategy. By utilizing quantized magnetization plateaus—when they occur—we show that it is not only possible to extract H at a much lower magnetic field, but these plateaus also exhibit unique, sharp spectroscopic signatures allowing for the detection of disorder with high precision. Similar to the QHE, this sensitivity emerges from the distinctive structure of incompressible magnetic states with quantized magnetization, an effect which has not been widely acknowledged until now.

At first glance, the magnetization process of antiferromagnets is simple. In the paramagnetic regime, Langevin or Brillouin theory describes the response of independent spins S—treated either as classical dipoles or quantum moments—to temperature T and applied magnetic field H15. In the magnetically ordered state, assuming that spins initially lie in the plane perpendicular to the magnetic field direction z, the classical dipoles cannot move toward the field direction. This induces a uniform magnetization M(H) which increases linearly until 〈Sz〉=S for each spin (Fig. 1a). The onset of this fully polarized state defines the saturation field Hsat. A closer inspection reveals that non-linear contributions to the isothermal magnetization encode rich phenomena. For example, quantum corrections result in a convex magnetization in Heisenberg S=1/2 antiferromagnets16,17. In geometrically frustrated systems, fractional magnetization plateaus can occur18–20 with 〈Sz〉=m×NS, where m is a rational number and NS is the total number of spins (Fig. 1b). In systems with biquadratic exchange, quantum and classical effects compete in the stabilization of magnetization plateaus21 while in frustrated S > 1/2 systems with uniaxial anisotropy, a Devil’s staircase with an infinite number of magnetization plateaus is possible at the classical level22. Clearly, magnetization plateaus are a widespread and fascinating occurrence in magnetic systems23. They reflect that the total spin operator is conserved and commutes with the Hamiltonian, i.e., [Stotz, H]=0. As the resulting magnetic excitations are gapped, a reliable determination of the system Hamiltonian H24,25 is possible at much lower magnetic fields than those necessary to reach a fully polarized state. Furthermore, regardless of the underlying classical or quantum mechanisms that stabilize the plateau, the resulting magnetic state is incompressible23; an increase in H does not lead to an increase in M(H). This incompressibility creates a highly sensitive platform for detecting even small amounts of chemical disorder, because spins near a defect can respond strongly to the magnetic field, i.e., these spins are compressible.

To demonstrate this spectacular effect, we focus on the Ising-Heisenberg (or XXZ) triangular-lattice antiferromagnet K2Co(SeO3)2 (KCSO). This compound has long been known26,27 and was recently revisited in detail28. It features a robust magnetization plateau that spans an extraordinarily wide field range, from 2 to 17 T29. The magnetization plateau is flanked by compressible spin supersolid phases at low temperature29–31. These phases are the quantum analogs of non-collinear Y and V phases, where a static magnetic order in the z-direction coexists with phase coherence—but no static order—transverse to z. Although the exact nature of these phases is not yet fully understood, we will refer to them as Y and V in the rest of this article. The crystal structure of KCSO, depicted in Fig. 1c, comprises isolated CoO6 octahedra connected by SeO3 tetrahedra, forming a triangular lattice of Co2+ ions in the ab-plane. These triangular layers are stacked along the c-axis and separated by layers of K ions. Although the effective magnetic moments of KCSO extracted from the Curie-Weiss fit at high temperatures are nearly isotropic, with μeffc=6.14μB ≈ 1.13μeffa (where μB is the Bohrmagneton) corresponding to an effective g-factor of g=7.128, bilinear interactions exhibit a strong Ising character. The out-of-plane nearest-neighbor exchange Jzz is an order of magnitude larger than the in-plane exchange Jxy29,30. Above the critical field μ0Hc1≈2 T (where μ0 is the vacuum permeability) and at or below T=5 K, KCSO stabilizes an archetypal up-up-down (UUD) phase associated with a magnetization plateau at m=1/3. This state persists until the critical field μ0Hc2≈17 T.

Results
Here, we utilize far-infrared magneto-optical spectroscopy (FIRMS) and thermomagnetic measurements to investigate the behavior of KCSO near and within the plateau at m=1/3 (see “Methods”). To our surprise, although our crystals are grown using published protocols and chemical disorder ruled out by advanced X-ray diffraction and electron microscopy methods (see Supplementary Section V and VI), the measurements reveal sharp spectroscopic fingerprints that can be clearly linked to very low levels of chemical disorder. These fingerprints were not seen in previous measurements30. Our optical measurements in the UUD phase reveal two dominant modes accompanied by three satellite modes, as illustrated in Fig. 1d, e. These satellite modes cannot be explained by clean-system models across three simulation techniques: linear spin wave theory (LSWT), Landau-Lifshitz dynamics (LLD), and exact diagonalization (ED). Instead, we develop a disordered-system model that fully captures our experimental observations. In this framework, the energy shift between the dominant and satellite modes reflects the nature and strength of exchange disorder on nearest-neighbor bonds. The large and discrete energy shifts in our data suggest that the disorder in our crystals is akin to missing bonds or spin vacancies, rather than a distribution of weakened bonds. As we illustrate below, a typically small amount of vacancies—as little as 2.3%—is sufficient to produce the observed intensity of the satellite modes. Moreover, magnetization and heat capacity measurements reveal a secondary magnetization plateau state with m<1/3, which is successfully reproduced in classical Monte Carlo (MC) simulations using the vacancy model.

Our FIRMS results, measured in Faraday geometry (H∥c≡z) at T=5 K, are shown in Fig. 1d. At low magnetic fields (0<μ0H<3 T), the system is in the Y phase30, and the spectra exhibit two broad field-dependent modes that cross at μ0H≈3 T. Around 3 T, the system enters the UUD phase and displays much sharper modes, with the two strongest intersecting at μ0H=9.75 T. Three satellite modes, running parallel to the main modes, are also observed. Above 15 T, an additional mode emerges and broadens with increasing field; we associate this mode with the onset of the V phase30. We determine the energy of major modes at each magnetic field, E(H), using Lorentzian fits, and plot the results on Fig. 1e, where the color encodes the parent Y, UUD, or V phase associated with each mode. Now, we address the origin of the two strongest modes in the UUD phase. The magnetic Hamiltonian of KCSO in that field configuration reads

H = Jxy ∑⟨i,j⟩ (Sx i Sx j + Sy i Sy j) + Jzz ∑⟨i,j⟩ Sz i Sz j + μ0 μB gzz H ∑i Sz i (1)

where the first sum runs over nearest neighbor bonds, Jxy > 0 and Jzz > 0 are antiferromagnetic interactions, and Si are spin-1/2 operators. The second term describes the usual Zeeman effect, with magnetic field H in the z-direction, and gzz representing the z-component of the effective g-factor for the Co2+ ions. In the classical Ising limit, the excited states of the UUD phase can be expressed as integer multiples of Jzz:

En = nJzz ± μ0 μB gzz H, (n=0,1,2,3,...) (2)

Due to optical selection rules, only transitions with ΔSz=±1 are visible in FIRMS. Therefore, our experiment primarily detects excited states involving single or triple spin-flip processes (see Supplementary Section II for further discussion). The two strongest modes observed in the UUD phase can be naturally attributed to single spin-flip excitations. In Fig. 2a, solid lines represent the mode energies E(H) expected for single spin-flips allowed by the selection rule. These processes can occur at two distinct spin sites: up spin (major spin) and down spin (minor spin). A minor spin flip costs onsite energy in terms of Jzz only—specifically, 3Jzz for the triangular lattice—since the flipped spin cannot propagate to neighboring spin sites. In contrast, a major spin flip incurs no onsite energy change in the classical limit but is accompanied by an additional quantum correction of 3Jxy/2, reflecting that the flipped spin can propagate to three neighboring up-spin sites (Fig. 2b, c). Linear fits to these two modes (see Supplementary Table 2) yield Jxy=0.20(1) meV, Jzz=2.88(1) meV, and gzz=7.45(5), in good agreement with previous neutron scattering and magnetization studies28–30, but providing a more accurate determination of the g-factor. Remarkably, in the UUD phase, all main and satellite modes appear resolution-limited, as shown in Supplementary Fig. 4. Triple spin-flip excitations, forming both a bound-state and a free particle continuum, are also observed (see Supplementary Section IIA for details).

Our data in the UUD phase exhibit three salient satellite modes that do not correspond to either single or triple spin-flip processes. To account for these excitations, we develop a disorder-bond model in which exchange interactions are modified as:

J'zz = Jzz(1−dzz), J'xy = Jxy(1−dxy) (3)

where dzz and dxy represent the disorder-induced reduction of the respective interaction. Notably, this model yields fine spectroscopic features that are the unmistakable hallmark of the disorder’s nature. Specifically, a major spin-flip excitation can split into two modes, depending on whether the disorder bond is connected to a major spin (mode #1a in Fig. 2) or a minor spin (mode #1b in Fig. 2). The excitation energies are shifted by ±Jzzdzz/2, respectively. Moreover, when the disorder bond is connected to a major spin, there is a quantum correction of −αJxy dxy (where α≃n is the disorder density in the dilute limit), due to the propagation of the flipped spin to two neighboring up-spin sites with intact bonds. In the case of a minor spin-flip process, the disorder only introduces an energy shift of −Jzzdzz/2, since the flipped spin cannot propagate to neighboring spin sites (see detailed discussion in Supplementary Section IID). More generally, we hypothesize that the splitting and multiplicity of spin-flip excitations in incompressible magnetic states serve as a unique indicator of bond disorder.

Using the analytical solution of the disorder-bond model, we perform linear fits to the observed satellite modes (see Supplementary Section IC). We find that the y-intercepts of these three modes are consistent with the model prediction, as illustrated in Fig. 2a. Quantitatively, we obtain dzz=0.93 from the energy shift between modes #1 and #1b, and dzz=0.96 from the energy shift between modes #2 and #2b. The close-to-unity value for dzz points to the complete absence of a bond, rather than a mere modification of the exchange due to disorder. The simplest mechanism to generate missing bonds is through spin-site vacancies, i.e., missing Co2+ ions.

The disorder-bond framework of Eq. (3) interpolates between three distinct models of random disorder: a weakened-bond model, a broken-bond model, and a vacancy model. Figure 3 schematically illustrates the spin and exchange configurations for each case and the associated excitation spectra. In the weakened-bond model, multiple bonds exhibit randomly distributed reduction factors dzz and dxy. The resulting satellite excitations detach from the main modes, with their energy offsets set by the reduction factors. Consequently, for a distribution of weakened bonds, we expect the main spectral lines to be broadened and sharp satellite modes to be absent, depending on the details of the disorder distribution (see Fig. 3b). By contrast, in the broken-bond model, the reduction factors are fixed to unity on each defective bond. This produces sharp satellite modes well separated from the main modes, which also remain sharp (see Fig. 3d). Remarkably, the vacancy model yields an identical spectrum, albeit with distinct intensity, because the six bonds adjacent to a vacancy are effectively removed (see Fig. 3e). Thus, despite their distinct chemical origins, the broken-bond and vacancy models are spectroscopically indistinguishable for an unknown defect density. From a materials standpoint, however, the complete absence of a specific bond is less realistic than the presence of a magnetic-ion vacancy.

Based on the above analysis and exchange parameters, and the fact that all main and satellite modes in the UUD phase are resolution-limited, we simulated the FIRMS spectra with and without vacancy disorder. We employed three different approaches, as summarized in Fig. 4: LSWT, LLD (to account for finite temperature), and ED (to account for quantum corrections). Calculations for a clean system, Fig. 4a–c, show that two sharp single spin-flip modes entirely dominate in the m=1/3 UUD phase, with only weak continua of triple spin-flip excitations in ED at high energies (see Supplementary Section IIA). Simulations including disorder paint a different picture, Fig. 4d–f. To account for disorder in LSWT and LLD, we employ a spin vacancy model with non-adjacent vacancies in supercells up to size 18×18. In ED, due to system size constraints, we harness the spectroscopic equivalence between broken bond and vacancy models and use the former with dzz=1 for two bonds in a 6×3 supercell, which equivalently reproduces the same excitations as the vacancy model (see Supplementary Section IIC). As the vacancy/disorder density is directly related to the intensity ratio between the dominant and satellite excitations, we vary the number of defects (in multiples of 3) in the (18×18) supercell until the simulations are in excellent agreement with the data (Supplementary Section IIB). We find that 2.3(8)% of Co2+ vacancies is the most likely range of concentrations, as it successfully reproduces the intensity of all satellite modes in the UUD phase. A more naive analytical approach to estimate the vacancy concentration (Supplementary Section IIE), which neglects spectroscopic matrix elements, yields 3.2(14)% of vacancies, within error bars of the more elaborate spin-wave-based method.

We now examine the impact of disorder on the thermomagnetic properties of our crystals of KCSO. Typical thermodynamic measurements are sensitive enough to quantify the impurity concentrations in diluted spin systems32,33, and here we extend this concept to a dense spin system. The isothermal magnetization M(H) measured to sub-Kelvin temperatures with the field along the c-axis is shown in Fig. 5a. In the clean case, the magnetization is expected to reach the m=1/3 plateau in a single step (see Fig. 1b for a cartoon), as observed by previous experimental studies29,30 (see also Supplementary Section III). In our case, however, an additional plateau with m<1/3 is observed from μ0H≈2.8 T, and the main plateau never reaches exactly m=1/3. To understand this behavior, we performed classical MC simulations for both our clean and disordered models of KCSO with a 2.5% vacancy concentration close to the representative value obtained from FIRMS measurements. In the clean case, the magnetization saturates at m=1/3 near μ0H≈4 T (Fig. 5b). In contrast, the magnetization in the disordered case shows a pre-plateau around μ0H≈2.8 T, followed by a jump to the primary plateau at 3.4 T (Fig. 5c). The magnetization on both plateaus falls below m=1/3. By tracking the spin configuration at each magnetic field in our MC simulations, we find that the pre-plateau originates from spin flips near vacancy sites (Fig. 5d). With this understanding, our naive analytical approach (Supplementary Section IIE) can be used to estimate the value of the normalized magnetization per spin. In the presence of a vacancy concentration n, it yields (1−n)/3 and (1−n)/3−2n for the magnetization at the main plateau and pre-plateau, respectively. This approach best matches our magnetization measurements for n=2.0% (Fig. 5a), which is remarkably close to the best concentration inferred from FIRMS results, giving further confidence in our vacancy-based disorder model. Thus, the additional plateau we observed serves as a thermodynamic signature of disorder. Taken independently, it lacks the mechanistic insights provided by our spectroscopic approach. But with insight about the disorder mechanism, it yields a surprisingly accurate determination of the disorder concentration.

Heat capacity measurements at low temperatures, with the magnetic field applied along the c-axis, are presented in Fig. 5f, g. The field dependence of specific heat at T=0.2 K reveal that a phase with gapless excitations exists below μ0H≈1 T, which we assign to the Y-phase or its quantum analog. As the field increases, the specific heat is suppressed but reaches a finite value again between μ0H≈2.1 T and 3.5 T (Fig. 5f). This enhancement survives up to T=0.5 K and matches the peaks observed in the field derivative of magnetization (Fig. 5e). This correlation suggests that the low-temperature specific heat can detect the disorder-induced plateau state through the density of states of local spin-flip processes near vacancies, as these become gapless. The value of the entropy release between T≈0.2 and 0.8 K near μ0H≈3 T is consistent with local spin-flip excitations near vacancies (see Supplementary Sec. IV). Above T=1 K, the heat capacity shows only a broad peak that corresponds to the field-driven transition from the Y phase to the UUD phase. These findings emphasize the significance of examining the low-temperature regime when investigating the effects of disorder, even in systems that exhibit Ising-limit behavior.

Methods
Sample preparation
Single crystals of K2Co(SeO3)2 were prepared following the procedure reported in ref. 28. Specifically, high-purity and commercially available K2CO3, Co3O4, and SeO2 were mixed in a molar ratio of 1.2:0.33:2.2 and loaded into an Al2O3 crucible, which was sealed in a vacuumed quartz tube. The tube was heated to T=600°C, maintained at that temperature for 8 h, and then cooled to room temperature at a rate of 6°C/hour. The resulting flat pink crystals were collected in water. The crystallographic c* and c axes are normal to the flat surface of the crystals, as confirmed by backscattering Laue X-ray diffraction. The largest K2Co(SeO3)2 crystal was selected for far-infrared magneto-spectroscopy (FIRMS) measurements and was polished to create parallel planes perpendicular to the c-axis. Two unpolished crystals were used for thermodynamic measurements, which yielded consistent results. Single crystals of Rb2Co(SeO3)2 were prepared similarly.

Far-infrared magneto-spectroscopy
FIRMS measurements were performed using a Bruker VERTEX 80v Fourier-transform infrared spectrometer at the National High Magnetic Field Laboratory (NHMFL, Tallahassee, Florida). The sample was mounted with GE varnish on a nonmagnetic aperture (1.5 mm in diameter) and attached to a top-loading probe. The probe was then inserted into an optical cryostat equipped with a superconducting magnet reaching up to μ0H=17.5 T. Infrared light from a mercury lamp was delivered to the sample via an evacuated light pipe and focused onto the sample surface using a parabolic cone.

To reduce light-induced heating and improve signal-to-noise ratio, a 700 cm−1 low-pass filter was placed in the light path just above the sample. The transmitted infrared light passing through the sample was collected by another pair of parabolic cones and detected by a Si bolometer operating at liquid helium temperature. A schematic of the resulting Faraday transmission configuration, with the light propagating along H∥c, is shown in Supplementary Fig. 1. The lowest effective sample temperature during measurements was Tmin=5 K, as measured by a thermometer mounted on the sample holder.

During the measurements, three different beamsplitters (Multilayer Mylar, Mylar 50 μm, and Mylar 125 μm) were used to cover the relevant energy range. The raw spectra from each (beam splitter) measurement are shown in Supplementary Fig. 2.

Spin dynamics simulation
Simulations of the transmitted frequency-dependent FIRMS spectra I(ω) were performed using the SUNNY.JL package37 and the QUSPIN.PY package38,39 by calculating the diagonal components of the frequency-weighted zero-momentum-transfer dynamical spin structure factor perpendicular to the light transmission direction,

I(ω)=ω[Sxx(q=0,ω)+Syy(q=0,ω)], (4)

where z∥H∥c. The dynamical spin structure factor was simulated using several approaches; in either case, magnetic moments were treated as dipolar spins and the system’s dynamics was governed by the same parent quantum Hamiltonian H, Eq. (1), and its known exchange parameters.

The zero-temperature response (T=0K) was simulated using linear spin wave theory, by which spin-1/2 operators Si are transformed into Holstein-Primakoff bosons bi and the resulting bosonic Hamiltonian, truncated at quadratic order, is diagonalized using known numerical approaches implemented in SUNNY.JL and other programs40,41.

The finite-temperature response was obtained in SUNNY.JL using classical Landau-Lifshitz dynamics (LLD). For these simulations, the spin system is assumed to be a product state |Ω⟩=⊗|Ωi⟩ over SU(2) coherent states |Ωi⟩ representing dipolar spin operators Si, which are replaced with classical vectors Ωi = ⟨Ωi|Si|Ωi⟩. The time-dependent dynamics at finite temperature is calculated using the stochastic Landau-Lifshitz-Gilbert (LLG) equation:

dΩi/dt = γΩi×(Hcl + ξ) + λΩi×(Ωi×Hcl), (5)

where Hcl = ⟨Ω|H|Ω⟩ in the large-S limit, ξ is a temperature-dependent Gaussian white noise, and λ=0.2 sets the coupling strength between the system and the thermal bath. The simulations were performed using a 9×9 super-cell of the conventional chemical unit cell of the crystal with periodic boundary conditions. To match the overall excitation bandwidth of our data, we set the temperature to T=0.5 K in the LLD simulations. The system was first thermalized over 5000 Langevin time steps with a step size of Δt=0.013 meV−1. After thermalization, five spin configurations were sampled using the LLG equation, with each configuration separated by 3000 Langevin time steps to ensure decorrelation.

Exact diagonalization calculations were performed using QUSPIN.PY38,39. The frequency-weighted dynamical spin structure factor, Eq. (4), was calculated using a 6×3 supercell with periodic boundary conditions, employing the Lanczos algorithm42 and the continued fraction method43. The number of Lanczos iterations was set to 200. Due to computational limitations, quenched chemical disorder was modeled using weak bonds rather than directly introducing vacancies. A more detailed discussion can be found in the Supplementary Section IIC.

Classical Monte Carlo simulation
Field-dependent magnetization was calculated from the standard Monte Carlo Metropolis algorithm using SUNNY.JL package37. The supercell size for the Monte Carlo simulation was set to 18×18×1. For each magnetic field, we perform 20000 spin-flip steps for the annealing at a given temperature. After the annealing process, 20 samples of spin configurations are averaged after 100 spin flips for each sample.

Thermomagnetic measurements
Isothermal magnetization measurements were performed using a Quantum Design MPMS-3 SQUID magnetometer operated in DC mode with the 3He insert option (temperatures from T=0.4 K to 4 K). The magnetic field was applied along the c-axis. Heat capacity measurements were conducted using a Quantum Design Dynacool PPMS equipped the Dilution Refrigerator unit. Magnetic fields up to μ0H=5 T were also applied along the c-axis to map the temperature- and field-dependent specific heat C(T,H) using the thermal relaxation method over a temperature range from T=200 mK to 4 K. Measurements above T=800 mK were performed at fixed magnetic fields with a high density of data points along the temperature axis, while measurements below T=800 mK were conducted at fixed temperatures with a higher density of data points along the magnetic field axis.

Single-crystal X-ray diffraction
Single-crystal X-ray diffraction (SCXRD) was performed on a specimen with dimensions 0.15×0.1×0.08 mm3. The crystal was affixed to a nylon loop using Paratone oil and examined with a Rigaku XtaLAB Synergy, Dualflex, Hypix SCXRD diffractometer, operated at room temperature.

Crystallographic data were acquired using the ω-scan method with MoKα radiation (λ=0.71073 Å) from a micro-focus sealed X-ray tube operating at ΔV=50 kV and I=1 mA. Experimental parameters, including the total number of runs and images, were determined algorithmically based on strategy computations performed by CrysAlisPro software (version 1.171.42.101a, Rigaku OD, 2023). A detailed analysis of the SCXRD data is provided in Supplementary Section V.

Scanning transmission electron microscopy
A focused ion beam from a Thermo Fisher Helios 5 CX was used to lift out a specimen from a K2Co(SeO3)2 single crystal. The specimen was then polished to a thickness of <50 nm to achieve atomic resolution using a Hitachi HD-2700 scanning transmission electron microscope (STEM). The microscope was operated at ΔV=200 keV with a convergence angle of 27 mrad and a spatial resolution of Δr≈1.3 Å. Annular dark-field (ADF) images were collected and compared with the K2Co(SeO3)2 lattice structure to determine the zone axis. A selected ADF image is shown in Supplementary Fig. 10.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This is an interesting paper on an important topic. The authors propose a novel method for probing disorder in a quantum magnetic system by analyzing its behavior within a magnetization plateau phase. The combination of high-resolution spectroscopy, thermodynamic measurements, and numerical modeling makes for a compelling and cohesive study. I am not sure the manuscript is suitable for publication in Nature Communications in its current state. However, I have several specific questions and comments which I believe, if addressed, would clarify key aspects of the analysis and significantly strengthen the conclusions.\n1- The central interpretation relies on a model involving a low concentration (2.5%) of magnetic site vacancies. However, this concentration is below the detection threshold of the structural measurements presented. It would be helpful to know whether the authors can estimate an uncertainty on this value, which appears to be inferred primarily from visual agreement between simulated and measured spectral intensities in the spin-wave calculations.\n2- In the methods section, it is mentioned that the crystals were separated from the flux using water. Could this step introduce surface modifications or an elevated density of defects near the surface, possibly influencing the spectroscopic measurements?\n3- The vacancy model provides a strong match to the measured spectra. Has the possibility of other forms of disorder - such as a distribution in exchange strengths rather than missing sites - been considered? If so, would such a scenario lead to distinct spectral signatures, or could it mimic the observed satellites?\n4- In the thermodynamic data, the low-temperature specific heat feature is attributed to localized excitations near vacancies. Could the authors provide a rough estimate of the entropy associated with this feature, based on their assumed vacancy density, and compare it with the experimentally integrated specific heat data?\n5- In the spectroscopic analysis shown in Figure 3, the simulated response contains only perfectly sharp modes. In contrast, the measured spectra exhibit notable linewidths. A brief discussion on the physical origins of this broadening—such as lifetime effects, unresolved fine structure, or instrumental contributions—would be valuable. Moreover, would this broadening affect the comparison of relative intensities between modes?\n6- In Figure 1e and the supplementary, there is a small but noticeable difference between the fitted slope values for the positive- and negative-field branches of the same excitation. What is the justification for using an averaged slope in the main figure? Could this asymmetry reflect an intrinsic anisotropy or a secondary interaction?\n7- In the supplementary material (specifically Equation S7), a quantum correction term is introduced that modifies the energy of one of the satellite modes. In the experimental data, the satellite and parent modes appear nearly parallel over the full field range. Does this constrain the correction parameter or place any limits on the acceptable disorder fraction?\nOverall, in my opinion, while the authors make a compelling case, the entire argument of \"metrological precision\" rests on a single proposed disorder type that remains indirectly inferred. In my opinion, the authors need to temper their claims to align with the indirect nature of their evidence, explicitly discuss the limitations and assumptions of their model, and strengthen the argument for the uniqueness of their interpretation. This is not minor clarifications and will require a rewriting of the manuscript's core narrative.", "reviewer_2": "The manuscript \"Sharp spectroscopic fingerprints of disorder in an incompressible magnetic state\" by C. Kim et al. presents a combined microscopic and thermodynamic study of the magnetization plateau in the up-up-down (UUD) state of the Ising triangular lattice antiferromagnet K2Co(SeO3)2. The authors observe sharp resonance modes within the UUD state that are not reproducible in a clean system model. They develop a new model incorporating spin disorder, which successfully captures these spectroscopic features and concurrently explains minor anomalies in the specific heat and magnetization data.\nIn my opinion, this is a very good work. The data are of high quality and have been properly analyzed. The manuscript is well-written and easy to follow, with a few exceptions mentioned below. The Supplementary Information provides all relevant additional data required to support the main findings. The authors performed very careful spectroscopic measurements and applied various state-of-the-art numerical methods to simulate the magnetic dynamics and developed a clear physical interpretation. They convincingly demonstrate the leading role of chemical disorder in the formation of the additional sharp resonance features within the plateau state. The results presented will be of great interest to the broad readership of Nature Communications and particularly to scientists working in the field of quantum magnetism. I recommend acceptance after minor revision, which primarily concerns the style of presentation and includes some optional suggestions for additional analysis.\nRemarks:\n1. In Fig. 1c, the authors represent the anisotropy of the exchange tensor using an elongated ellipsoid. While this visualization may be intuitive for experienced users of software like SpinW, it is likely more confusing than informative for the general condensed matter audience.\n2. \"In this semi-classical regime, becomes possible to directly measure H [9], including the level of chemical disorder [10]. \" Reference [10] does not seem right here because author do provide an accurate estimate of the disorder strength in YbMgGaO using the field-polarized INS data.\n3. The two high-energy modes present in Fig. 1d but absent in Fig. 1e may be confusing at first glance. While they are addressed in the main text and SI, the authors should clarify directly in the figure caption that these features are also understood within their modeling framework (SI, Sec. IIa,b).\n4. The authors discuss two distinct models for disorder—\"broken-bonds\" and \"site-vacancies\"—and transition between them in the main text. I found this confusing until reaching SI Sec. IIc, where they demonstrate the equivalence of these models in generating the dynamical response. Given that the two scenarios are physically different and that the analysis identifies site-vacancies as the more probable mechanism, I strongly recommend the authors restructure the narrative to prioritize this interpretation throughout the main text. The site-vacancy model should be presented as the primary physical model. To avoid conceptual ambiguity for the reader the broken-bonds mechanism should be mentioned much less. Probably only in application to the ED simulations, where it was used as a technical trick to help with computations.\n5. Magnetization plateaus. In Figs 4a and S9 authors show magnetization data which do not reach 1/3. The magnetization is normalized to Co2+. Was the 2.5 % site-vacancies incuded in the normalization?\n6. The observation of an additional magnetization plateau at ~2-4 T and its subsequent physical interpretation is a particularly exciting result. Could there be additional spectral features associated with (presumably local) resonances of the spins close to the vacancy?\n7. The samples of K2Co(SeO3)2 studied in Refs [26,27] do not show anomalies near the plateau [Fig. S9] and presumably are cleaner. It would be very interesting (probably in the future) to show that the additional spectral lines do not appear in those samples and thus provide a direct experimental evidence for the disorder-based theory.\n8. While the XRD and STEM data did not directly observe atomic-scale vacancies, the spectroscopic and thermodynamic evidence for disorder is compelling. This suggests the disorder may not be in the form of vacancies but could instead arise from isovalent substitution of Co²⁺ by a non-magnetic ion of a similar size. If this is the case, techniques sensitive to elemental identity, such as neutron diffraction or resonance x-ray spectroscopy, could be used to identify the specific impurity in future studies.\n9. As a small remark to the general conclusion - there is an alternative thermodynamics-based approach to the disorder problem in quantum magnets when starting on the opposite end - measurements of specific heat of magnetic dimers in spin diluted systems. It was successfully applied to several systems (Yb:YAlO3 [PhysRevB.101.24515] and Yb:NaLuSe2 [PhysRevB.106.024404]) and allowed to accurately detect very low concentration of magnetic impurities with high fidelity.", "reviewer_3": "Kim et al. present a comprehensive study on the magnetic properties of a frustrated quantum spin system. They probe the textbook triangular lattice antiferromagnet K2Co(SeO3)3 by a coherent combination of optical spectroscopy, magnetization and specific heat measurements and consistently reproduce the results by theoretical modeling. In the optical data they find two main modes in the far-infrared that are accompanied by three satellite peaks, which they assign to magnetic defects. Thus they conclude a vacancy spin density of 2.5%, which are assigned to Co sites as pairwise alignment prohibits propagation of flipped spins to neighboring sites due to Pauli principle. From the theoretical model both a clean system and one with small disorder levels are compared, yielding absence of satellites in the former and good agreement of the latter with experiment – again very reasonable. The same comparison between experiment and theory is also made for magnetization and specific heat, revealing the particular influence of disorder on these bulk properties probed at zero frequency. For magnetization the impurity spins show up as a finite slope within a magnetization plateau, i.e. which is nicely explained in terms of compressible (impurity) and gapped (bulk) response. Finally, spectroscopic fingerprints of disorder are seen also in a sister compound (Rb2Co(SeO3)3), providing the results more generality.\nOverall, I share the excitement of the authors about the in-depth determination of magnetic disorder type and density by spectroscopic probes. While the optical measurements surely provide information way beyond standard probes (magnetization, specific heat), the detection of the impurity density from a slope in an otherwise flat magnetization plateau is an innovative tool. Below my comments:\na) On page 3 the authors claim 2.5% of vacancy spins as a ‘typically undetectable amount’. While I would grossly agree with this statement for standard probes (magnetization, specific heat), there are certainly other methods that can probe even smaller amounts of magnetic disorder. In particular electron spin resonance has proven highly sensitive, e.g. in organic spin liquid systems [J. Phys. Soc. Jpn. 65, 1340–1354 (1996), Science 372, 276–279 (2021)]. It is very nice that the authors present a new way of spectroscopic disorder probing, but also other established spectroscopic probes should be acknowledged.\nb) Following up on the point above, 2.5% of disorder can have profound effects on magnetic properties, even concealing the underlying magnetic ground states (see e.g. [PRB 101, 140401(R) (2020), PRB 108, 104410 (2023)]). Thus, I’d not agree that 2.5% is a small amount, which should be toned down, albeit it is less than in herbertsmithite with 5-10% of disorder [J. Am. Chem. Soc. 132, 16185–16190 (2010)].\nc) In Fig. 4 the temperature is indicated with color scale bars. While this is a generally nice approach, especially useful for large data sets, here only few curves are shown, which would allow to write down a legend with temperatures within the panels. This way the reader directly has the full information and does not have to guess the precise temperature value from a color scale.\nd) Lastly, I appreciate a lot that the authors probed a sister compound to verify the findings. From the manuscript and supplement I understood that several samples have been measured. Are there experimental data on several samples, such that one can quantitatively compare the disorder density from the intensity of optical peaks with magnetization (and also specific heat)? While the manuscript in its present form does fit well to Nature Communications, a quantitative comparison of different experimental probes on several samples and with theoretical modeling, may even raise the bar to Nature Physics."}, "author_response": "Date: November 1, 2025\nManuscript: NCOMMS-25-63285-T Kim\nPoint-by-point response:\nReviewer 1\nGeneral Comment 1. This is an interesting paper on an important topic. The authors propose a novel\nmethod for probing disorder in a quantum magnetic system by analyzing its behavior within a magneti-\nzation plateau phase. The combination of high-resolution spectroscopy, thermodynamic measurements,\nand numerical modeling makes for a compelling and cohesive study. I am not sure the manuscript is\nsuitable for publication in Nature Communications in its current state. However, I have several specific\nquestions and comments which I believe, if addressed, would clarify key aspects of the analysis and\nsignificantly strengthen the conclusions.\nOur Answer: Thank you.\nComment 1.1. The central interpretation relies on a model involving a low concentration (2.5%) of\nmagnetic site vacancies. However, this concentration is below the detection threshold of the structural\nmeasurements presented. It would be helpful to know whether the authors can estimate an uncertainty\non this value, which appears to be inferred primarily from visual agreement between simulated and\nmeasured spectral intensities in the spin-wave calculations.\nOur Answer: Thank you for your question regarding the uncertainty in the vacancy concentration. Be-\ncausethesupercellsizesetstheavailablevacancyfractioninoursimulations,finelysamplingconcentration\nspace is computationally expensive. Since the experimental intensity ratios between the main peaks and\nsatellitepeaksarebetweenthevaluespredictedbysimulationswith3.1%and1.5%vacancyconcentration,\nwe estimate the vacancy concentration in our sample to be 2.3(8)%, noting that 2.5% is the most repre-\nsentative. In the revised manuscript, we explain how this value was obtained and include an additional\nback-of-the-envelope calculation that yields a similar result of 3.2(14)%. Furthermore, we have extended\nour back-of-the-envelope analysis of the isothermal magnetization to estimate the vacancy concentration\nfromthereductioninUUDmagnetizationcomparedtoacleansystem. Thisanalysisalsoyieldsconsistent\nresults.\nComment 1.2. In the methods section, it is mentioned that the crystals were separated from the flux\nusing water. Could this step introduce surface modifications or an elevated density of defects near the\nsurface, possibly influencing the spectroscopic measurements?\nOur Answer: Thank you for your question regarding the possible origin of the defects. After separation,\nthe samples were polished (using the standard polishing procedure) to enable FIRMS measurements.\nGiven this preparation, it is unlikely that the observed spectroscopic signatures arise from surface de-\nfects introduced by the flux separation with water. Moreover, our thermodynamic measurements were\nperformed on unpolished specimens. Because the optical and thermodynamic results are mutually con-\nsistent with vacancy-type disorder, we conclude that the vacancy signal does not originate from extrinsic\n1 / 9\nsample-preparationprocedures. WehaveupdatedourMethodssectiontoexplainthepolishingprocedure.\nComment 1.3. The vacancy model provides a strong match to the measured spectra. Has the possibility\nof other forms of disorder - such as a distribution in exchange strengths rather than missing sites -\nbeen considered? If so, would such a scenario lead to distinct spectral signatures, or could it mimic the\nobserved satellites?\nOur Answer: Thank you for your question regarding the disorder model. This is an important point, as\nreal materials often exhibit some degree of bond disorder, which can give rise to a distribution of exchange\ncouplings. In our case, however, such a scenario is inconsistent with the experimental data (or, at most,\nnot dominant) due to the sharpness and the discrete energy shift of the observed satellite peaks. First,\nif the exchange couplings were broadly distributed, then the satellite peaks would emerge as a collection\nof sidebands whose individual energy separation directly follows the distribution of exchange couplings.\nThis distribution would produce a continuum-like spectral weight filling the energy range between the\nmain and satellite peaks, rather than the distinct satellites that are observed. In the weak bond-disorder\nlimit, this would effectively broaden the main peak, making it appear smeared rather than sharp. As we\nexplain in Comment 1.5, our observed main and satellite peaks are both at or close to the experimental\nresolution limit. Second, in order for a distribution of exchange strengths to reproduce the energy shift\nof the observed satellite peaks, the distribution would need to only include entirely suppressed exchange\ninteractions (i.e., 100% reduction). For these reasons, the sharp satellites observed in the data are incom-\npatible with a model based on randomly distributed exchange strengths. We have added this discussion\nto the revised manuscript.\nComment 1.4. In the thermodynamic data, the low-temperature specific heat feature is attributed to\nlocalizedexcitationsnearvacancies. Couldtheauthorsprovidearoughestimateoftheentropyassociated\nwith this feature, based on their assumed vacancy density, and compare it with the experimentally\nintegrated specific heat data?\nOur Answer: Thank you for the comment on the entropy analysis. From our vacancy model, the expected\nentropy release due to the field-alignment of spins surrounding a vacancy is as follows:\n∆S = (Vacancy %)×(# of U spin in UUD phase)×(# of spin flip near the vacancy)×Rln2\n= 0.023×2/3×3×Rln2 ≃ 0.046×Rln2,\nwhere R denotes the universal gas constant.\nIn the revised manuscript, we compare this prediction with the experimental entropy release between\nT = 0.2 K and T ≈ 0.8 K for different applied magnetic fields in the UUD phase. The experimental\nentropy change, at its maximum near µ H ≈ 3.2 T, reaches ≈ 0.05 × Rln2 (see Fig. R1), consistent\n0\nwith the above expectation. We thank the reviewer for the illuminating comment, which improved the\nmanuscript. We have added this discussion and Fig. R1 to the revised manuscript.\nComment 1.5. In the spectroscopic analysis shown in Figure 3, the simulated response contains only\nperfectly sharp modes. In contrast, the measured spectra exhibit notable linewidths. A brief discussion\non the physical origins of this broadening - such as lifetime effects, unresolved fine structure, or instru-\nmental contributions - would be valuable. Moreover, would this broadening affect the comparison of\nrelative intensities between modes?\nOur Answer:\nThank you for the comment about the linewidth analysis. In Fig. R2 below, we compare the fitted\n2 / 9\nFig.R1: EntropyofKCSOfromourspecificheatdata. Theminimumtemperaturewas0.2K.Thedashed\nline indicates our expected entropy from the localized spin flip near the vacancy.\nFig. R2: Fitted full-width half-maximum (FWHM) of the main peaks and satellite peaks. The numbering\nconvention of the peaks (i.e., Peak 1 and Peak 2) is the same as the main text. The blue circles indicate\nthe FWHM of the main peaks, while the cyan triangles show the FWHM of the satellite peaks.\nfull-width half-maximum (FWHM) of the main peaks and satellite peaks in the UUD phase. We find that\nthe FWHM of main and satellite peaks are similar and close to the best case scenario for the instrumental\nresolution estimated to be around 0.3–0.5 cm−1 (or 0.04–0.06 meV). This observation implies that the\nleading origin of the satellite peaks cannot be broadly distributed bond disorder. We have added this\ndiscussion and Fig. R2 to the revised manuscript.\n3 / 9\nComment 1.6. In Figure 1e and the supplementary, there is a small but noticeable difference between\nthe fitted slope values for the positive- and negative-field branches of the same excitation. What is the\njustification for using an averaged slope in the main figure? Could this asymmetry reflect an intrinsic\nanisotropy or a secondary interaction?\nOur Answer: Thank you for the comment about the difference in the fitted slope values. We believe the\ndifference originates from systematic errors in the fitting, in particular the fitting range. In the original\nversion of the manuscript, we only considered statistical error bars reported by the fitting routine. In the\nreevaluatedversion,wereportsystematicerrorsonthefittedslopes. Nowthefittedslopesinbothpositive-\nand negative-field branches are consistent within error bars. Of course, we cannot exclude that additional\ntermsinthematerialsHamiltonian,suchassubleadinganisotropiesorinteractions,mayimpacttheslopes.\nComment 1.7. In the supplementary material (specifically Equation S7), a quantum correction term is\nintroduced that modifies the energy of one of the satellite modes. In the experimental data, the satellite\nand parent modes appear nearly parallel over the full field range. Does this constrain the correction\nparameter or place any limits on the acceptable disorder fraction?\nOur Answer: Thank you for the comment regarding the quantum correction. The reviewer is correct\nthat a quantum correction only impacts the 1a mode. However, the quantum correction only corrects the\nenergy (i.e., the zero-field intercept) of this mode, not the slope (as shown in the Supplementary Table 2).\nConsequently, the satellite 1a mode is expected to run parallel to the parent mode throughout the entire\nfield range.\nConcluding Comment 1. Overall, in my opinion, while the authors make a compelling case, the entire\nargument of \"metrological precision\" rests on a single proposed disorder type that remains indirectly\ninferred. In my opinion, the authors need to temper their claims to align with the indirect nature of\ntheir evidence, explicitly discuss the limitations and assumptions of their model, and strengthen the\nargument for the uniqueness of their interpretation. This is not minor clarifications and will require a\nrewriting of the manuscript’s core narrative.\nOur Answer: Thank you for the insightful comments and constructive criticism. We have taken into\naccount the reviewer’s valuable comments and greatly extended the discussion about the choice of the\ndisordermodel. Themanuscriptis writtenin amorelogicalway now. Inparticular, we added an estimate\nof the uncertainty of the vacancy density and an explicit comparison of the FWHM of observed excita-\ntions in the UUD phase. As we stated above in Comment 1.3, these refinements show that our vacancy\ndisorder model is not one of many possibilities to explain the data among those disorder types, but the\nonlysimplemodelthatcancomprehensivelyexplainourspectroscopicandthermodynamicmeasurements.\n4 / 9\nReviewer 2\nGeneral Comment 2. The manuscript \"Sharp spectroscopic fingerprints of disorder in an incompressible\nmagnetic state\" by C. Kim et al. presents a combined microscopic and thermodynamic study of the\nmagnetization plateau in the up-up-down (UUD) state of the Ising triangular lattice antiferromagnet\nK2Co(SeO3)2. The authors observe sharp resonance modes within the UUD state that are not\nreproducible in a clean system model. They develop a new model incorporating spin disorder, which\nsuccessfully captures these spectroscopic features and concurrently explains minor anomalies in the\nspecific heat and magnetization data.\nIn my opinion, this is a very good work. The data are of high quality and have been properly ana-\nlyzed. The manuscript is well-written and easy to follow, with a few exceptions mentioned below. The\nSupplementary Information provides all relevant additional data required to support the main findings.\nThe authors performed very careful spectroscopic measurements and applied various state-of-the-art\nnumerical methods to simulate the magnetic dynamics and developed a clear physical interpretation.\nThey convincingly demonstrate the leading role of chemical disorder in the formation of the additional\nsharp resonance features within the plateau state. The results presented will be of great interest to\nthe broad readership of Nature Communications and particularly to scientists working in the field of\nquantum magnetism. I recommend acceptance after minor revision, which primarily concerns the style\nof presentation and includes some optional suggestions for additional analysis.\nOur Answer: Thank you.\nComment 2.1. InFig. 1c, theauthorsrepresenttheanisotropyoftheexchangetensorusinganelongated\nellipsoid. While this visualization may be intuitive for experienced users of software like SpinW, it is\nlikely more confusing than informative for the general condensed matter audience.\nOur Answer: Thank you for the comment about the visualization of the anisotropy. We understand that\nthe ellipsoid can confuse the general audience with different physical quantities, such as a nematic order\nparameter and so on. We removed the ellipsoid in our revised manuscript.\nComment 2.2. \"In this semi-classical regime, becomes possible to directly measure H [9], including the\nlevel of chemical disorder [10].\" Reference [10] does not seem right here because author do provide an\naccurate estimate of the disorder strength in YbMgGaO using the field-polarized INS data.\nOur Answer: Thank you. In the revised manuscript, we cite two YbMgGaO papers investigating the\n4\neffects of disorder [Li et al. Phys. Rev. Lett. 118, 107202 (2017) and Steinhardt et al. Phys. Rev. Res.\n3, 033050 (2021)].\nComment 2.3. The two high-energy modes present in Fig. 1d but absent in Fig. 1e may be confusing at\nfirst glance. While they are addressed in the main text and SI, the authors should clarify directly in the\nfigure caption that these features are also understood within their modeling framework (SI, Sec. IIa,b).\nOur Answer: Thank you for the comment about the triple spin-flip excitations. We now show the fitted\ntriple spin-flip excitations with purple color in the revised Fig. 1e.\n5 / 9\nComment2.4. Theauthorsdiscusstwodistinctmodelsfordisorder-\"broken-bonds\"and\"site-vacancies\"-\nand transition between them in the main text. I found this confusing until reaching SI Sec. IIc, where\nthey demonstrate the equivalence of these models in generating the dynamical response. Given that the\ntwo scenarios are physically different and that the analysis identifies site-vacancies as the more probable\nmechanism, I strongly recommend the authors restructure the narrative to prioritize this interpretation\nthroughout the main text. The site-vacancy model should be presented as the primary physical model.\nTo avoid conceptual ambiguity for the reader, the broken-bonds mechanism should be mentioned much\nless. Probably only in application to the ED simulations, where it was used as a technical trick to help\nwith computations.\nOur Answer: Thank you for the suggestion to prioritize the vacancy model in the main text. In response\nto this comment, as well as Reviewer 1’s feedback, we have modified the manuscript in depth. The key\nidea we now emphasize is that the “broken-bond” and “site-vacancy” models are identical/equivalent in\nthe spectroscopic sense. To explain this clearly, we have added a new figure in the main text to compare\nthree distinct disorder models: random bond disorder, broken-bond, and vacancy models.\nComment 2.5. Magnetization plateaus. In Figs 4a and S9 authors show magnetization data which do\nnot reach 1/3. The magnetization is normalized to Co2+. Was the 2.5 % site-vacancies incuded in the\nnormalization?\nOur Answer: Thank you for the question about the magnetization normalization. The reviewer is correct\nthat the magnetization was normalized to the Co2+ ion. So, the site vacancy effect was not considered in\nthe normalization. We have corrected this, however, due to the tiny vacancy concentration, the change in\nnormalization is not clearly visible on the figures.\nComment 2.6. The observation of an additional magnetization plateau at 2-4 T and its subsequent\nphysical interpretation is a particularly exciting result. Could there be additional spectral features\nassociated with (presumably local) resonances of the spins close to the vacancy?\nOur Answer: Thank you for the question about the additional spectral features of the vacancy. A natural\nquestion is the finite-momentum dispersion of the vacancy-induced excitations. In the UUD phase, since\nthe spin-flip process remains localized near the vacancy, the finite-momentum dispersion of these excita-\ntions is expected to be almost flat. We have not yet explored how near-vacancy excitations behave in\nother phases of the material, such as the proposed supersolid phases, which could be an interesting topic\nfor future research.\nComment 2.7. The samples of K2Co(SeO3)2 studied in Refs [26,27] do not show anomalies near the\nplateau [Fig. S9] and presumably are cleaner. It would be very interesting (probably in the future)\nto show that the additional spectral lines do not appear in those samples and thus provide a direct\nexperimental evidence for the disorder-based theory.\nOur Answer: Thankyouforthecommentonthesamplequality. Beforethesubmissionofthemanuscript,\nwe discussed with one of the authors of the recent study of K2Co(SeO3)2, featuring a clean sample,\narXiv:2402.15869 (2024). The discussion indicated that the Johns Hopkins team also observed satellite\npeaks in some samples, but measured several samples and successfully picked a clean sample without\nany satellite peaks for their study. The time-domain terahertz spectroscopy (TDTS) results for the clean\nsample are now included in the SI of arXiv:2402.15869 (2024). These results are compared to our FIRMS\nresults in Fig. R3, where we note that satellite peaks are absent in TDTS for the clean sample. The\nabsence of visible satellite peaks in these measurements strengthens our claim that the satellite modes we\n6 / 9\nobserve are the direct and unequivocal experimental evidence of the disorder.\nFig. R3: Comparison of magneto-optical spectroscopy between our disordered sample (a) and clean-limit\nsample (b). Panel (b) is adopted from the Johns Hopkins study [arXiv:2402.15869 (2024)].\nComment 2.8. While the XRD and STEM data did not directly observe atomic-scale vacancies, the\nspectroscopic and thermodynamic evidence for disorder is compelling. This suggests the disorder may\nnot be in the form of vacancies but could instead arise from isovalent substitution of Co2+ by a non-\nmagnetic ion of a similar size. If this is the case, techniques sensitive to elemental identity, such as\nneutron diffraction or resonance x-ray spectroscopy, could be used to identify the specific impurity in\nfuture studies.\nOur Answer:\nWe thank the reviewer for this interesting suggestion. While we agree that substitutional disorder is an\nintriguing possibility, we consider it unlikely in our case, given the details of the sample synthesis pro-\ncess. Nevertheless, we agree that targeted substitution studies would be highly valuable. In particular,\nintroducing non-magnetic ions of similar ionic size to Co2+, such as Zn2+, could provide an informative\navenue for future research.\nComment 2.9. As a small remark to the general conclusion - there is an alternative thermodynamics-\nbased approach to the disorder problem in quantum magnets when starting on the opposite end -\nmeasurements of specific heat of magnetic dimers in spin diluted systems. It was successfully applied\nto several systems (Yb:YAlO3 [PhysRevB.101.24515] and Yb:NaLuSe2 [PhysRevB.106.024404]) and al-\nlowed to accurately detect very low concentration of magnetic impurities with high fidelity.\nOur Answer: Thank you for the suggestion about the references. This is a nice suggestion to broaden the\nscope of the study to the characterization and detection of magnetic “color centers” in quantum magnets.\nWe have cited these papers in the revised manuscript.\n7 / 9\nReviewer 3\nGeneral Comment 3. Kim et al. present a comprehensive study on the magnetic properties of a frus-\ntrated quantum spin system. They probe the textbook triangular lattice antiferromagnet K2Co(SeO3)3\nby a coherent combination of optical spectroscopy, magnetization and specific heat measurements and\nconsistently reproduce the results by theoretical modeling. In the optical data they find two main modes\nin the far-infrared that are accompanied by three satellite peaks, which they assign to magnetic defects.\nThus they conclude a vacancy spin density of 2.5%, which are assigned to Co sites as pairwise alignment\nprohibits propagation of flipped spins to neighboring sites due to Pauli principle. From the theoretical\nmodelboth aclean system andone with smalldisorder levels are compared, yielding absence ofsatellites\nin the former and good agreement of the latter with experiment - again very reasonable. The same com-\nparison between experiment and theory is also made for magnetization and specific heat, revealing the\nparticular influence of disorder on these bulk properties probed at zero frequency. For magnetization the\nimpurity spins show up as a finite slope within a magnetization plateau, i.e. which is nicely explained\nin terms of compressible (impurity) and gapped (bulk) response. Finally, spectroscopic fingerprints of\ndisorder are seen also in a sister compound (Rb2Co(SeO3)3), providing the results more generality.\nOverall, I share the excitement of the authors about the in-depth determination of magnetic disorder\ntype and density by spectroscopic probes. While the optical measurements surely provide information\nway beyond standard probes (magnetization, specific heat), the detection of the impurity density from\na slope in an otherwise flat magnetization plateau is an innovative tool. Below my comments:\nOur Answer: Thank you.\nComment 3.1. On page 3 the authors claim 2.5% of vacancy spins as a ’typically undetectable amount’.\nWhileIwouldgrosslyagreewiththisstatementforstandardprobes(magnetization, specificheat), there\nare certainly other methods that can probe even smaller amounts of magnetic disorder. In particular\nelectron spin resonance has proven highly sensitive, e.g. in organic spin liquid systems [J. Phys. Soc.\nJpn. 65, 1340-1354 (1996), Science 372, 276-279 (2021)]. It is very nice that the authors present a\nnew way of spectroscopic disorder probing, but also other established spectroscopic probes should be\nacknowledged.\nOur Answer: Wethankthereviewerforpointingthisout. Indeed,whileavacancyconcentrationof∼2.5%\nwould be undetectable in standard bulk probes such as magnetization or specific heat, we fully agree that\nmore sensitive spectroscopic techniques—including electron spin resonance (ESR, a sister technique of\nour FIRMS measurement) and nuclear magnetic resonance (NMR)—can detect even smaller amounts\nof magnetic disorder. In particular, the reviewer’s suggested references nicely illustrate the remarkable\nsensitivity of ESR in probing dilute spin defects in quantum magnets. We have now acknowledged this\npoint in the revised manuscript and cited the suggested references accordingly.\nComment 3.2. Following up on the point above, 2.5% of disorder can have profound effects on magnetic\nproperties, evenconcealingtheunderlyingmagneticgroundstates(seee.g. [PRB101, 140401(R)(2020),\nPRB108,104410(2023)]). Thus,I’dnotagreethat2.5%isasmallamount,whichshouldbetoneddown,\nalbeit it is less than in herbertsmithite with 5-10% of disorder [J. Am. Chem. Soc. 132, 16185-16190\n(2010)].\nOur Answer: We thank the reviewer for this important clarification. We agree that even a 2.5% level of\ndisorder can substantially affect magnetic properties and, in some cases, even obscure the intrinsic ground\nstate, as demonstrated in suggested references. We have cited these useful references and revised the\n8 / 9\nmanuscript accordingly.\nComment 3.3. In Fig. 4 the temperature is indicated with color scale bars. While this is a generally\nnice approach, especially useful for large data sets, here only few curves are shown, which would allow\nto write down a legend with temperatures within the panels. This way the reader directly has the full\ninformation and does not have to guess the precise temperature value from a color scale.\nOur Answer: We thank the reviewer for this helpful suggestion. In the revised manuscript, we have added\nexplicit temperature legends to each panel of Fig. 4, so that the precise temperature values can be read\ndirectly rather than estimated from the color scale bar.\nComment 3.4. Lastly, Iappreciatealotthattheauthorsprobedasistercompoundtoverifythefindings.\nFrom the manuscript and supplement I understood that several samples have been measured. Are there\nexperimental data on several samples, such that one can quantitatively compare the disorder density\nfrom the intensity of optical peaks with magnetization (and also specific heat)? While the manuscript\nin its present form does fit well to Nature Communications, a quantitative comparison of different\nexperimental probes on several samples and with theoretical modeling, may even raise the bar to Nature\nPhysics.\nOur Answer: Thank you. Unfortunately, we measured only one batch of samples: the largest crystal\nwas used for FIRMS measurements, a smaller piece for magnetization, and another for specific heat. Our\nvacancy model provides a consistent explanation for all experimental results. We have clarified this point\nin the revised Methods section.\n9 / 9"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The revised manuscript represents a careful and thoughtful response to previous reviewer concerns. While some assumptions remain only qualitatively justified, the authors have appropriately tempered their claims, clarified methodological limitations, and strengthened the manuscript with additional analysis. I now believe the paper is suitable for publication.", "reviewer_2": "The authors have successfully addressed all of my questions and have significantly improved the quality of the manuscript. I believe the current version is now suitable for publication in Nature Communications.", "reviewer_3": "The Authors have adequately addressed my requests and I recommend publication in Nature Communications.\nIn view of my comment 3.4, I encourage the authors to perform sample-dependent experiments (different samples from the same batch and also comparison of different batches) in future studies, which will provide in important insights into the disorder density."}, "author_response": "Date: November 20, 2025\nManuscript: NCOMMS-25-63285-A Kim\nPoint-by-point response:\nReviewer 1\nGeneral Comment 1. The revised manuscript represents a careful and thoughtful response to previous\nreviewer concerns. While some assumptions remain only qualitatively justified, the authors have appro-\npriately tempered their claims, clarified methodological limitations, and strengthened the manuscript\nwith additional analysis. I now believe the paper is suitable for publication.\nOur Answer: Thank you.\nReviewer 2\nGeneral Comment 2. The authors have successfully addressed all of my questions and have significantly\nimproved the quality of the manuscript. I believe the current version is now suitable for publication in\nNature Communications.\nOur Answer: Thank you.\nReviewer 3\nGeneral Comment 3. The Authors have adequately addressed my requests and I recommend publication\nin Nature Communications. In view of my comment 3.4, I encourage the authors to perform sample-\ndependent experiments (different samples from the same batch and also comparison of different batches)\nin future studies, which will provide in important insights into the disorder density.\nOur Answer: Thank you for the comments and the nudge to continue our studies on different samples.\nWe are planning to study the sample dependence of vacancy excitations in the near future.\n1 / 1"}]

Decision letter