Nonperturbative Semiclassical Spin Dynamics for Ordered Quantum Magnets

TL;DR

This paper develops non-perturbative methods to analyze spin dynamics in ordered quantum magnets, successfully explaining anomalous continua and bound states observed in experiments, especially in the triangular-lattice antiferromagnet.

Contribution

It introduces the THED approach and combines it with MPS simulations to interpret complex spin excitation spectra in quantum magnets.

Findings

Confirmed two-magnon bound states in the system

Identified the anomalous continuum as a two-magnon resonance

Validated results against inelastic neutron scattering data

Abstract

In ordered quantum magnets where interactions between elementary excitations dominate over their kinetic energy, perturbative approaches often fail, making non-perturbative methods essential to capture spectral features such as bound states and the redistribution of weight within excitation continua. Although an increasing number of experiments report anomalous spin excitation continua in such systems, their microscopic interpretation remains an open challenge. Here, we investigate the spin dynamics of the triangular-lattice antiferromagnet in its 1/3-plateau phase using two complementary non-perturbative approaches: exact diagonalization in a truncated Hilbert space for a gas of elementary excitations (THED) and matrix product state (MPS) simulations. Alongside cross-validation between these methods, we benchmark our results against inelastic neutron scattering (INS) data. The THED analysis confirms the presence of two-magnon bound states and identifies the anomalous scattering continuum observed in both MPS and INS as a two-magnon resonance, arising from hybridization between the bound state and the two-magnon continuum. Furthermore, THED reveals bound states overlapping with the continuum, enriching the interpretation of continuum anomalies. More broadly, THED provides a robust framework for investigating anomalous spin excitation continua and bound-state effects in other materials with gapped spectra. Its combination of accuracy and computational efficiency makes it a powerful tool for extracting reliable microscopic models in semiclassical regimes.

Equations

$$\mathcal{H}=\sum _{⟨i,j⟩}{J}_{ij}(S_i^x{S}_j^x+S_i^y{S}_j^y+{\Delta }_{ij}{S}_i^z{S}_j^z)-{\mu }_B{\mathbf{B}}^T\mathbf{g}\sum _j{\mathbf{S}}_j.$$

$${\mathcal{S}}^{\alpha \alpha }(\mathbf{q},\omega )=\frac{1}{2\pi N_s}{\int }_{-\infty }^{+\infty }dt\sum _{j,l}^{N_s}{e}^{i\left(\omega t-\mathbf{q}\cdot ({\mathbf{r}}_j-{\mathbf{r}}_l)\right)}⟨S_j^{\alpha }(t)S_l^{\alpha }(0)⟩,$$

Full text

Modern condensed matter physics rests on the notion of quasiparticles, which allows for the recasting of complex many-body dynamics in terms of weakly interacting collective excitations. In ordered quantum magnets, these excitations manifest as quantized spin waves, or magnons. When magnon-magnon interactions are weak, diagrammatic perturbation theory (loop expansion) provides a powerful framework for describing their dynamics, accurately capturing features such as spectral renormalization and intrinsic linewidth broadening [1, 2, 3, 4, 5].

As interactions become comparable to or exceed the magnon kinetic energy, the system can no longer be treated as a weakly interacting gas and perturbative approaches break down. In this regime, quasiparticles can form bound states or resonances, undergo strong renormalization, or decay via channels involving bound states [6, 7, 8, 9, 10, 11]. These strong interactions substantially reshape the structure of the multi-magnon continuum, making it essential to accurately capture non-perturbative effects in order to extract reliable microscopic models. Furthermore, the softening of bound states under variation of Hamiltonian parameters often signals the onset of exotic multipolar orderings [12, 11, 13]. Despite a growing number of experiments reporting anomalous spin excitation continua in ordered magnets [4, 14, 15, 16], their microscopic origin remains largely unresolved, partly due to the scarcity of non-perturbative theoretical tools capable of capturing these complex many-body dynamics.

Recent methodological advances in the time evolution of matrix product states (MPS) [17, 18, 19, 20] have enabled the calculation of spectral functions that can be directly compared with experimental measurements, such as inelastic neutron scattering (INS). These techniques have proven highly effective in capturing complex non-perturbative effects in ordered quantum magnets [14, 21, 16]. However, being purely numerical approaches, MPS-based methods often provide limited interpretative insight into the underlying physical mechanisms. Moreover, accurate MPS simulations in two or more spatial dimensions remain computationally challenging, limiting the applicability of these methods for extracting effective microscopic models of real materials.

In this article, we test and extend the truncated Hilbert space exact diagonalization (THED) method, originally developed to capture non-perturbative single-particle effects in the effective spin-one magnet $\mathrm{FeI_{2}}$ [10, 9], to investigate multi-particle continua in more general gapped ordered quantum magnets. By restricting the full Hilbert space to sectors containing a limited number of magnons, THED captures few-body non-perturbative effects while reducing the computational complexity from exponential to polynomial scaling. This approach enables the resolution of non-perturbative spectral features with clear physical interpretation, provided the system possesses a finite single-magnon gap. Our calculations are carried out using a recently developed general-purpose implementation of the method [22].

Here, we test our method on the spin-1/2 triangular-lattice Heisenberg antiferromagnet (TLHAFM) in the presence of a magnetic field, where a robust 1/3-magnetization plateau emerges. This field-induced phase, originally predicted by Chubukov and Golosov [23], hosts a distinctive “up-up-down” (UUD) spin configuration (Fig. 1a) and has since been observed across a wide range of spin-1/2 triangular-lattice compounds [24, 25, 26, 27, 28, 29, 30, 31, 32, 33]. Long recognized as a hallmark of frustrated magnetism, the plateau not only reflects the stabilization of quantum order [34, 35, 36, 37, 38] but also provides a fertile ground for exploring strongly interacting spin dynamics. A persistent puzzle in this context is the presence of an extended excitation continuum above the single-magnon modes that evades explanation by conventional spin-wave theory [4, 15].

While recent MPS simulations have succeeded in reproducing these broad spectral features [14], they offer limited insight into the underlying microscopic processes and the nature of the modes. The origin of the continuum therefore remains an open question, pointing to the need for new theoretical frameworks capable of capturing the inherently non-perturbative nature of spin excitations in this regime.

Beyond its intrinsic interest, the study of excitation dynamics in the 1/3-magnetization plateau phase is also directly relevant to a family of spin-1/2 $J_{1}$-$J_{2}$ triangular-lattice compounds, $A\mathrm{YbSe_{2}}$ (with $A=\mathrm{Na,K,Cs}$) [39, 40, 41, 14], which reside near a quantum critical regime separating magnetically ordered states from a putative spin-liquid ground state [40]. While semiclassical methods are often appropriate in the field-induced UUD phase, capturing the essential features of the excitation spectrum requires incorporating non-perturbative effects. Accurate modeling in this regime is crucial not only for understanding the anomalous spin dynamics but also for reliably extracting microscopic Hamiltonian parameters, which serve as inputs to zero-field studies probing quantum criticality [41].

In this work, we perform parallel simulations using both the THED and MPS methods to enable cross-validation. In the THED approach, we include contributions up to the two-magnon sector to capture the anomalous continua observed in INS experiments, accounting for the finite single-magnon gap generated by quantum fluctuations in the UUD phase [34, 4]. For the MPS simulations, we carry out calculations in a cylindrical geometry with moderate bond dimensions and cluster sizes, taking advantage of the weak entanglement characteristic of the UUD phase. Technical details of both approaches are provided in the “Methods” section.

By combining these two approaches and benchmarking them against INS data, we identify the anomalous continuum in the UUD phase of TLHAFM as a two-magnon resonance arising from the hybridization between a bound state and the multi-magnon continuum. Furthermore, the interpretability afforded by the THED method allows us to uncover bound states embedded within the continuum, which are reminiscent of the bound states in the continuum (BICs) of photonic systems [42, 43].

Beyond the INS examples presented in this work, the predictions of the THED method can be directly compared with a range of other spectroscopic probes of spin dynamics for ordered quantum magnets, including Raman scattering, electron spin resonance, optical conductivity, and resonant inelastic X-ray scattering.

More broadly, the THED framework–performing exact diagonalization while constraining the number of elementary excitations–can be applied to any $N$-level system exhibiting semiclassical behavior. This includes atomic gases in optical lattices [44] serving as quantum simulators, quantum optical systems [45], and orbital degrees of freedom in solids [46, 47]. In these contexts, THED offers a computationally efficient approach to capturing non-perturbative phenomena while keeping the numerical complexity tractable.

Results

The spin Hamiltonian for the triangular lattice antiferromagnet (TLAFM) is given by

[TABLE]

The first term describes the antiferromagnetic (AFM) exchange interaction, with $J_{ij}>0$ acting on the first- and second-nearest neighbor bonds, while $\Delta_{ij}$ is the dimensionless exchange anisotropy. The second term represents the Zeeman interaction, expressed as the contraction of the external magnetic field $\bm{B}$ with the $g$-tensor $\bm{g}$ and the spin operator $\bm{S}_{j}$ of site $j$. Throughout the main text, we either assume that the external field is aligned with the easy-axis (taken to be the $z$-axis), or we employ an isotropic Heisenberg model with $\Delta_{ij}=1$, which preserves the U(1) symmetry of the Hamiltonian for arbitrary field directions. Our main goal is to investigate the low-energy spin excitations of this model in the field-induced three-sublattice UUD ordered phase.

The compound $\mathrm{KYbSe_{2}}$ (KYS) provides a representative realization of the model $\mathcal{H}$ with parameters $J_{1}=0.456\ \text{meV}$, $J_{2}/J_{1}\simeq 0.043$, and $\Delta_{ij}=1$. Figure 1b shows the magnetization curve of KYS under an external magnetic field applied along the $x$-axis, with the 1/3-plateau phase highlighted in orange. Our main results for KYS in this plateau phase at $B=4\penalty 10000\ \mathrm{T}$ are presented in Fig. 1d, where we compare the INS cross-section with simulations performed using THED and MPS methods along the $(H,H,3)$ direction (see Fig. 1c for the momentum path). Notably, both methods reproduce the key spectral features observed in the INS data.

Three bright branches are observed below 1.0 meV, corresponding to collective single-magnon excitations originating from the three sublattices. These modes are well described by including the leading $1/S$ correction to linear spin wave theory (LSWT) [34, 36, 4]. However, as we will show below, the spectral weight between 1.0 and 1.6 meV exhibits a structured two-magnon continuum that cannot be captured by any finite-order $1/S$ expansion [4, 15]. Remarkably, both the THED and MPS methods accurately reproduce these non-perturbative features. A central message of this work is that, beyond capturing such effects, the THED approach provides complementary insights to MPS by elucidating the internal structure of the two-magnon continuum.

To understand the nature of the two-magnon excitations, we begin by considering the classical Ising limit ($\Delta_{ij}\to\infty$). For simplicity, we restrict the discussion to nearest-neighbor interactions with $J_{zz}=J_{1}$, unless otherwise noted. The inclusion of a small $J_{2}/J_{1}$ ratio leads to only quantitative modifications, without altering the qualitative physics. Figure 2 summarizes the relevant excitations involving up to two spin flips in the classical Ising limit.

We first analyze two-spin-flip processes occurring on both the “up” and “down” sublattices. As illustrated in Fig. 2, flipping two neighboring spins incurs a lower energy cost than flipping them farther apart, specifically $\delta E_{c}<\delta E_{a}+\delta E_{b}$, leading to the formation of two-magnon bound states. These excitations carry a total spin change $\Delta S^{z}=0$, making them observable via INS. As we show below, the anomalous continua observed in the INS spectra of materials near the isotropic limit can be traced back to these bound-state excitations in the Ising limit.

Another two-spin-flip process, shown in Fig. 2d, corresponds to a distinct bound state with $\Delta S^{z}=-2$. This quadrupolar excitation is invisible to INS within the current model. However, as demonstrated in Supplementary Sec. IV, such bound states can acquire finite spectral weight if the U(1) symmetry along the easy-axis in Eq. (1) is broken.

As the system evolves from the Ising limit toward the isotropic point ($\Delta=1$), two-magnon bound states become increasingly dispersive and eventually overlap with the continuum formed by pairs of single-magnon excitations. To capture this evolution, we adopt the approach of Ref. [14], tracking the redistribution of spectral weight as a function of the exchange anisotropy $\Delta$. In contrast to previous studies, however, the THED method employed here offers a key advantage: it provides direct access to the many-body eigenstates within the truncated Hilbert space, enabling a more transparent identification of the underlying excitations.

Figure 3 presents a schematic evolution of the two-magnon spectrum extracted from THED simulations as $\Delta$ is varied. Each two-magnon state is classified by its total spin change: $\Delta S^{z}=2$ when both spins flip from “down” to “up,” $\Delta S^{z}=-2$ for two “up” to “down” flips, and $\Delta S^{z}=0$ when one spin flips in each direction (color-coded in Fig. 3).

Finite $xy$ interactions lift the degeneracy of the classical Ising energy $\delta E_{c}$ associated with the $\Delta S^{z}=0$ bound state, splitting it into four distinct branches. The highest-energy branch corresponds to a state with lattice angular momentum $l=0$ and energy $\delta E=\delta E_{c}+J$, while the lowest-energy branch has $l=\pi$ and $\delta E=\delta E_{c}-J$. The remaining two branches form quasi-doublets with angular momenta $l=\pm\pi/3$ and $l=\pm 2\pi/3$, located at $\delta E=\delta E_{c}+J/2$ and $\delta E=\delta E_{c}-J/2$, respectively. In the Ising limit, all four branches appear as sharp bound states, lying between the $\Delta S^{z}=-2$ and $\Delta S^{z}=0$ continua.

As $\Delta$ decreases, the lower two branches are the first to enter the $\Delta S^{z}=-2$ continuum. However, their quantum numbers prohibit hybridization with the continuum, resulting in BICs. By contrast, the upper two branches merge into the $\Delta S^{z}=0$ continuum at smaller $\Delta$, where hybridization is allowed, converting them into two-magnon resonances. At the Heisenberg point, the lower branches persist as BICs or as bound states below the continuum, while the upper branches remain as resonant features.

The schematic picture discussed above is based on identifying the lower and upper bounds of the continuum and isolated states from THED, neglecting their dispersion. The actual dispersions of the two-magnon eigenstates obtained from THED for $\Delta=5$, $2$, and $1$ are shown in the first column of Fig. 4 along the $[H,H,0]$ direction. In the two rightmost columns of Fig. 4, we present the calculated longitudinal dynamical spin structure factor (DSSF), defined as

[TABLE]

with $\alpha=z$ for the model Hamiltonian Eq. (1) with only nearest-neighbor exchange $J=J_{1}$ and $g\mu_{B}B/J_{1}=1.905$, where $\bm{B}$ is applied along the $z$-direction and $N_{s}$ denotes the number of sites. These calculations are carried out using both THED and MPS methods. The longitudinal DSSF predicted here is experimentally accessible: it can be measured directly with polarized INS, or extracted from unpolarized INS provided the $g$-tensor is anisotropic, as illustrated by the examples in Fig. 5.

For $\Delta=5$, we find good agreement between the longitudinal dynamical structure factor $\mathcal{S}^{zz}(\bm{q},\omega)$ computed using THED and MPS. Both methods reveal four nearly flat and well-defined two-magnon branches along the $[H,H,0]$ direction. A real-space analysis of the corresponding Fourier-transformed wavefunctions obtained from THED and presented in the Supplementary Sec. III, shows that given a spin flip on the “up” (“down”) sublattice, the second spin flip on the “down” (“up”) sublattice is strongly localized to the nearest-neighbor sites. This spatial confinement supports the identification of these modes as tightly bound two-magnon states.

For $\Delta=2$, the lower two branches persist as BICs, in line with the discussion above, while the upper two evolve into resonances, appearing as broadened features in both THED and MPS spectra. Upon reaching the isotropic point $\Delta=1$, parts of the lower branches detach from the continuum along the $[H,H,0]$ direction, while others remain embedded. Nevertheless, due to their distinct quantum numbers, they retain their character as sharp BICs. The upper branches continue to manifest as resonances.

The real-space wavefunctions of these resonant states, also analyzed in Supplementary Sec. III, exhibit extended spatial profiles: unlike bound states, whose wavefunction decays rapidly beyond a characteristic binding length, the resonance wavefunctions maintain significant amplitude at large separations, consistent with their delocalized nature.

Figure 4 shows that the agreement between THED and MPS remains robust at the Heisenberg point. However, constant-$\bm{q}$ cuts, provided in Supplementary Figs. S5, S6, and S7, highlight minor discrepancies between the two methods. First, as shown in Supplementary Figs. S5 and S6, THED is not expected to capture the small renormalization of the single-magnon dispersion arising from hybridization with three-magnon states. This explains the slight mismatch in single-magnon peak positions at $\Delta=1$, which could be resolved by extending the THED variational space to include three-magnon states. The smallness of this mismatch suggests that the maximum number of magnons provides a good control parameter for systematically expanding the low-energy DSSF.

For the two-magnon sector revealed by the longitudinal DSSF shown in Fig. 4, the MPS calculation reveals a weak intensity mode around $E/J=1.5$ in the isotropic limit $\Delta=1$, which does not appear in the THED results. This mode exhibits a parabolic dispersion, with its minimum at $H=1/3$ coinciding with that of the two-magnon bound state near $E/J=0.75$. This suggests that the additional mode seen in the MPS data may correspond either to a four-magnon bound state or to the onset of the two two-magnon bound-state continuum. Such an interpretation would also explain its absence in the present THED implementation, which includes only single and two-magnon states. However, as discussed in Supplementary Sec. II, this mode could also be an artifact of the cylindrical geometries used in MPS.

Another discrepancy emerges at $\Delta=2$ in the energy window $4.5\lesssim E/J\lesssim 5.5$. The MPS results display a splitting between two modes for $H<1/3$, which is absent in the THED calculation. Furthermore, while both approaches exhibit a single dominant high-intensity mode for $H>1/3$, the mode identified by MPS appears at higher energy. The detailed analysis in Supplementary Sec. II indicates that these discrepancies may once again stem from the proximity to the four-magnon continuum associated with the two-magnon bound state. At the same time, this analysis shows that finite-size effects due to the cylindrical geometry are pronounced in this region of the spectrum, suggesting that the MPS results may not fully capture the two-dimensional limit of the model.

The above discussion allows us to reinterpret the INS data for KYS. The structured continua observed in Fig. 1d above the single-magnon branches originate from the spectral weight of two-magnon resonances, arising from the hybridization between two-magnon bound states and the two-magnon continuum. Due to their proximity to the single-magnon branches, the spectral weight of the two-magnon bound states overlap with that of the single-magnon modes, making them difficult to distinguish directly in INS (see Supplementary Fig. S5 for detailed line cuts).

A key advantage of the Yb-based compound $\mathrm{CsYbSe_{2}}$ (CYS), described by $\mathcal{H}$ with $J_{1}=0.395\penalty 10000\ \text{meV}$ and $J_{2}/J_{1}\simeq 0.03$, is the pronounced anisotropy of its $\bm{g}$-tensor ($g_{aa}=g_{bb}=3.25$, $g_{cc}=0.2$). This anisotropy enables the direct extraction of the longitudinal DSSF from INS intensities. In Fig. 5a–c, we compare the measured $\mathcal{S}^{xx}(\bm{q},\omega)$ with spins polarized along the $x$-axis at $B=4\penalty 10000\ \mathrm{T}$ (deep inside the UUD phase) to theoretical results obtained from THED and MPS. The agreement corroborates our interpretation: sharp two-magnon bound states dominate the low-energy spectrum, while broader features at higher energies are consistent with two-magnon resonances.

Finally, we turn to the material $\mathrm{K_{2}Co(SeO_{3})_{2}}$ (KCS) in its UUD phase, which realizes the model in Eq. (1) with nearest-neighbor interactions close to the Ising limit. The strong agreement between INS data and theoretical predictions from both THED and MPS simulations provides further validation of the two-magnon bound state picture. To reproduce the experimentally observed positions of the bound-state modes, we refine the exchange anisotropy from the previously reported value $\Delta=J_{zz}/J_{xy}=14.28$ [48, 49] to $\Delta=13.48$, while keeping $J_{xy}=0.217\ \text{meV}$ fixed 111The original estimate of $\Delta$ was obtained in Ref. [49] by fitting magnetization and specific heat using first-order perturbation theory in $J_{xy}$. The small discrepancy likely stems from higher-order corrections not accounted for in that analysis.. This refinement underscores the utility of our approach in accurately extracting microscopic model parameters. Moreover, as shown in the Supplementary Fig. S12, the absence of spectral weight associated with the quadrupolar bound state reinforces the conclusion that KCS is well described by an effective U(1)-symmetric Heisenberg model.

Discussion

We have presented a comprehensive investigation of non-perturbative dynamics in the 1/3-plateau (UUD) phase of triangular lattice antiferromagnets. By combining two complementary numerical methods, truncated Hilbert space exact diagonalization and matrix product states, we demonstrate consistent and robust results, despite the distinct limitations of each approach: truncation errors in THED and finite-size effects in MPS. Benchmarking against inelastic neutron scattering experiments confirms the reliability of both methods for modeling quantum magnets in this regime.

The high expressiveness of MPS makes it a powerful tool for probing non-perturbative phenomena in two-dimensional ordered magnets, while the interpretability of THED allows us to disentangle the structure of anomalous continua and identify underlying bound-state features.

Our results indicate that previously collected INS data on triangular lattice compounds in the UUD phase merit re-examination in light of this refined theoretical framework. Moreover, the role of higher-order $n_{\text{mag}}$-magnon states ($n_{\text{mag}}\geq 3$), omitted in the current THED implementation, underscores the need for methodological improvements, a feasible task given the polynomial scaling of the numerical complexity. At the same time, efforts in advancing the MPS algorithms may enable simulations on wider cylinders, thereby reducing the finite-size effects introduced by the cylindrical geometry. Taken together, the combined use of THED and MPS provides a powerful and broadly applicable framework for studying non-perturbative excitations in ordered quantum magnets with gapped single-magnon spectra. Looking ahead, we anticipate that this dual approach will be especially valuable in contexts where conventional methods face limitations.

In this work, we focus on the canonical example of the triangular lattice in the UUD phase. More generally, the THED method is expected to be applicable to a broader class of gapped quantum magnets that exhibit semiclassical behavior. The two-magnon bound states and resonances discussed here can only be obtained through an infinite-order resummation in the $1/S$ expansion, where $S$ labels the irreducible representations (irreps) of SU(2). For systems with $S\geq 1$ or for coupled entangled units [51], the classical limit of certain Hamiltonians, such as those with large single-ion anisotropy or weak interactions between entangled units, is more appropriately described by direct products of SU($N$) coherent states, where $N$ is the number of levels of each unit. In such cases, the $1/S$ expansion must be replaced by a $1/\lambda_{1}$ expansion, where $\lambda_{1}$ labels the completely symmetric irreps of SU($N$) [52, 53, 54, 55, 56]. In terms for Feynman diagrams, the $1/\lambda_{1}$ expansion is simply the well-known loop expansion. In other words, the order of each diagram in powers of $1/\lambda_{1}$ is linear in the number of loops.

The quadratic Hamiltonian obtained to order $\lambda_{1}$ corresponds to a generalized spin-wave theory that naturally incorporates different types of quasiparticles. For instance, these quasiparticles may correspond to triplon modes in weakly coupled dimer systems [51] or to dipolar and quadrupolar waves in $S=1$ systems [5, 8, 10, 9]. Similar to the $1/S$ expansion, the formation of two-quasiparticle bound states requires an infinite resummation of ladder diagrams with an increasing number of loops. Both perturbative and non-perturbative effects have been demonstrated in effective $S=1$ ordered magnets [5, 8, 10, 9]. In particular, higher-order bound states, such as four- and six-magnon states, have been shown to arise from non-perturbative effects in the excitation spectrum of FeI2 [10, 9].

To demonstrate the broader applicability of the THED approach, we have developed an open-source implementation [22] within the Sunny.jl spin-dynamics framework [57, 58]. This framework provides a user-friendly interface for model definition and many-body calculations, following a workflow analogous to semiclassical linear spin-wave theory. Thanks to its flexibility and computational efficiency, THED offers a powerful tool for investigating gapped quantum magnets with semiclassical behavior in three dimensions, where MPS becomes impractical. Furthermore, as an efficient and accurate forward-scattering approach, THED can be seamlessly integrated into inverse scattering frameworks, particularly in light of recent advances in machine learning algorithms [59].

Acknowledgments

We thank Shang-Shun Zhang, Kipton Barros, David Dahlbom, Martin Mourigal, Chaebin Kim, and Xiaojian Bai for useful discussions. H.Z. and C.D.B. thank Bruce Normand for his insightful comments and critical reading of the manuscript. C.D.B. acknowledges support from the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division, under Award No. DE-SC-0018660. H.Z. acknowledges support from the Lincoln Chair of Excellence in Physics. The work by A.S. is supported by the Quantum Science Center (QSC), a National Quantum Information Science Research Center of the U.S. Department of Energy (DOE). Work at ETHZ was partially supported by a MINT grant of the Swiss National Science Foundation. Data on $\mathrm{K_{2}Co(SeO_{3})_{2}}$ were collected using the AMATERAS spectrometer at J-PARC in Experiment no. 2023B0161. T. X. was supported by the National Key Research and Development Program of China (Grant No. 2024YFA1613100), the National Natural Science Foundation of China (Grant No. 12304187), the Guangzhou Basic and Applied Basic Research Funds (Grant No. 2024A04J4024), and the open research fund of Songshan Lake Materials Laboratory (Grant No. 2023SLABFN30). The neutron scattering experiments of CsYbSe2 were performed on the time-of-flight (ToF) Cold Neutron Chopper Spectrometer (CNCS) at the Spallation Neutron Source at Oak Ridge National Laboratory (ORNL).

Author contributions

C.D.B. and H.Z. conceived the project. H.Z. developed the numerical code for the THED calculations and analyzed the results together with C.D.B.. T.H. and A.M.L. carried out the MPS simulations. A.S. performed the INS experiments and analyzed the data for KYS. T.X. provided the INS data for CYS. M.Z. and A.Z. conducted the INS experiments and analyzed the data for KCS. N.M. and S.O.-K. assisted with neutron scattering experiments at J-PARC for KCS. H.Z., T.H., A.M.L., and C.D.B. wrote the manuscript with input from all authors.