The nature of spin excitations in the one-third magnetization plateau phase of Ba3CoSb2O9

Magnetization plateaus in quantum magnets—where bosonic quasiparticles crystallize into emergent spin superlattices—are spectacular yet simple examples of collective quantum phenomena escaping classical description. While magnetization plateaus have been observed in a number of spin-1/2 antiferromagnets, the description of their magnetic excitations remains an open theoretical and experimental challenge. Here, we investigate the dynamical properties of the triangular-lattice spin-1/2 antiferromagnet Ba3CoSb2O9 in its one-third magnetization plateau phase using a combination of nonlinear spin-wave theory and neutron scattering measurements. The agreement between our theoretical treatment and the experimental data demonstrates that magnons behave semiclassically in the plateau in spite of the purely quantum origin of the underlying magnetic structure. This allows for a quantitative determination of Ba3CoSb2O9 exchange parameters. We discuss the implication of our results to the deviations from semiclassical behavior observed in zero-field spin dynamics of the same material and conclude they must have an intrinsic origin.

Notwithstanding the progress in the search of quantum plateaus, much less is known about their excitation spectra. Given that they are stabilized by quantum fluctuations, it is natural to ask if these fluctuations strongly affect the excitation spectrum. The qualitative difference between the plateau and the classical orderings may appear to invalidate spin-wave theory. For instance, the UUD state in the equilateral TLHAFM is not a classical ground state unless the magnetic field H is fine-tuned 40 . Consequently, a naive spin wave treatment is doomed to instability. On the other hand, spin wave theory builds on the assumption of an ordered moment |〈S r 〉| close to the full moment. Given that a sizable reduction of |〈S r 〉| is unlikely within the plateau because of the gapped nature of the spectrum, a spin wave description could be adequate. Although this may seem in conflict with the order-by-disorder mechanism 1-3 stabilizing the plateau [40][41][42] , this phenomenon is produced by the zero-point energy correction E zp ¼ ð1=2Þ P q ω q þ OðS 0 Þ (ω q is the spin wave dispersion), which does not necessarily produce a large moment size reduction.
Here, one of our goals is to resolve this seemingly contradictory situation. Recently, Alicea et al. developed a method to fix the unphysical spin-wave instability 40 . This proposal awaits experimental verification because the excitation spectrum has not been measured over the entire Brillouin zone for any fluctuationinduced plateau. We demonstrate that the modified nonlinear spin wave (NLSW) approach indeed reproduces the magnetic excitation spectrum of Ba 3 CoSb 2 O 9 within the 1/3 plateau [30][31][32][33][34][35][36][37][38][39] . The excellent agreement between theory and experiment demonstrates the semiclassical nature of magnons within the 1/3 plateau phase, despite the quantum fluctuationinduced nature of the ground state ordering. The resulting model parameters confirm that the anomalous zero-field dynamics reported in two independent experiments 37,39 must be intrinsic and non-classical.

Results
Overview. In this article, we present a comprehensive study of magnon excitations in the 1/3 magnetization plateau phase of a quasi-two-dimensional (quasi-2D) TLHAFM with easy-plane exchange anisotropy. Our study combines NLSW theory with in-field inelastic neutron scattering (INS) measurements of Ba 3 CoSb 2 O 9 . The Hamiltonian is where 〈rr′〉 runs over in-plane nearest-neighbor (NN) sites of the stacked triangular lattice and c 2 corresponds to the interlayer spacing (Fig. 1a). J (J c ) is the antiferromagnetic intralayer (interlayer) NN exchange and 0 ≤ Δ < 1. The magnetic field is in the in-plane (x) direction (we use a spin-space coordinate frame where x and y are in the ab plane and z is parallel to c). h red = g ⊥ μ B H is the reduced field and g ⊥ is the in-plane g-tensor component.
This model describes Ba 3 CoSb 2 O 9 (Fig. 1b), which comprises triangular layers of effective spin 1/2 moments arising from the J ¼ 1=2 Kramers doublet of Co 2+ in a perfect octahedral ligand field. Excited multiplets are separated by a gap of 200-300 K due to spin-orbit coupling, which is much larger than the Néel temperature T N = 3.8 K. Below T = T N , the material develops conventional 120°ordering with wavevector Q = (1/3, 1/3, 1) 30 . Experiments confirmed a 1/3 magnetization plateau for Hjjab (Fig. 1c) 31,33,35,36,38 , which is robust down to the lowest temperatures. We compute the dynamical spin structure factor using NLSW theory in the 1/3 plateau phase. We also provide neutron diffraction evidence of the UUD state within the 1/3 plateau of Ba 3 CoSb 2 O 9 , along with maps of the excitation spectrum obtained from INS.
Quantum-mechanical stabilization of the plateau in quasi-2D TLHAFMs. While experimental observations show that deviations from the ideal 2D TLHAFM are small in Ba 3 Co Sb 2 O 9 33,37,39 , a simple variational analysis shows that any J c > 0 is enough to destabilize the UUD state classically. Thus, a naive spin wave treatment leads to instability for J c > 0. However, the gapped nature of the spectrum 5 implies that this phase must have a finite range of stability in quasi-2D materials. This situation must be quite generic among fluctuation-induced plateaus, as they To put this into a proper semiclassical framework, we apply Alicea et al.'s trick originally applied to a distorted triangular lattice 40 . Basically, we make a detour in the parameter space with the additional 1/S-axis quantifying the quantum effect (Fig. 2). Namely, instead of expanding the Hamiltonian around S → ∞ for the actual model parameters, we start from the special point, J c = 0, h red = 3JS, and a given value of 0 ≤ Δ ≤ 1, where the UUD state is included in the classical ground state manifold. Assuming the spin structure in Fig. 1a, we define for r ∈ A e , B e , A o , and C o and for r ∈ C e , B o . Introducing the Holstein-Primakoff bosons, a y ð Þ μ;r , with 1 ≤ μ ≤ 6 being the sublattice index for A e , B e , C e , A o , B o , and C o in this order, we havẽ andS À r ¼S þ r À Á y for r ∈ μ, truncating higher order terms irrelevant for the quartic interaction. We evaluate magnon selfenergies arising from decoupling of the quartic term. As shown in Fig. 2b, the linear spin wave (LSW) spectrum for J c = 0 and h red = 3JS features two q-linear gapless branches at q = 0, both of which are gapped out by the magnon-magnon interaction (Fig. 2c). Small deviations from J c = 0 and h red = 3JS do not affect the local stability of the UUD state because the gap must close continuously. Thus, we can investigate the excitation spectrum of quasi-2D systems for fields near h red = 3JS. Figure 2d, e show the spectra for h red shifted by ±10% from h red = 3JS, where we still keep J c = 0. For h red < 3JS, a bandtouching and subsequent hybridization appear between the middle and the top bands around q = (1/6,1/6) (Fig. 2d). For h red > 3JS, a level-crossing between the middle and bottom bands appears at around q = 0 (Fig. 2e). A small J c > 0 splits the three branches into six (Fig. 2f, g). Figure 3a, b show the reduction of the sublattice ordered moments for S = 1/2, J c = 0, 0.09J, and selected values of Δ. We find δ S x μ D E =S≲30% throughout the local stability range of the plateau. Thus, our semiclassical approach is fully justified within the plateau phase. Figure 3c, d show the field dependence of the staggered magnetization, which is almost field-independent; a slightly enhanced fieldindependence appears for small Δ. Similarly, while the magnetization is not conserved for Δ ≠ 1, it is nearly pinned at 1/3 for the most part of the plateau (Fig. 3e, f).
As discussed below, most of the broadening stems from the different intensities of the split modes due to finite J c .
Comparing the experiment against the NLSW calculation, we find that the features of the in-plane spectrum in Fig. 5a are roughly captured by the theoretical calculation near the lowfield onset of the plateau in Fig. 2f (h red = 2.7JS ≈ 1.03h red,c1 ). This observation is in accord with the fact that the applied field (μ 0 H = 10.5 T) is close to μ 0 H c1 = 9.8 T 32 . To refine the quantitative comparison, we calculate the scattering intensity where F(q) denotes the magnetic form factor of Co 2+ corrected with the orbital contribution, (γr 0 /2) 2 is a constant,q α ¼ q α = q j j are the diagonal components of the dynamical structure factor evaluated at 10.5 T; off-diagonal components are zero due to symmetry. Defining the UUD order as shown in Fig. 1a, transverse spin fluctuations related to single-magnon excitations appear in S yy and S zz , while longitudinal spin fluctuations corresponding to the two-magnon continuum appear in the inelastic part of S xx , denoted as S jj . Accordingly, I tot (q, ω) can be separated into transverse I ⊥ and longitudinal I jj contributions. To compare with our experiments, the theoretical intensity is convoluted with momentum binning effects (only for I ⊥ ) and empirical instrumental energy resolution. and g ⊥ = 3.95. The agreement between theory and experiment is excellent. When deriving these estimates, J is controlled by the saturation field μ 0 H sat = 32.8 T for Hjjĉ 33 . To obtain the best fit, we also analyzed the field dependence of ω 1 , ω 2 , and ω 3 (Fig. 6).
Remarkably, the calculation in Fig. 5d-f reproduces the dispersions almost quantitatively. It predicts a gapped ω 1 mode, although the gap is below experimental resolution. The smallness of the gap is simply due to proximity to H c1 . For each ω i , the band splitting due to J c yields pairs of poles ω ± i dispersing with a phase difference of π in the out-of-triangular-plane direction (Fig. 5e, f). For each pair, however, one pole has a vanishing intensity for q = (1/2, 1/2, l). Consequently, ω 1 along this direction is free from any extrinsic broadening caused by overlapping branches (Fig. 5e), yielding a relatively sharp spectral line (Fig. 5h). The corresponding bandwidth ≈ 0.2 meV (Fig. 5b) provides a correct estimate for J c . By contrast, for q = (1/3, 1/3, l), all six ω ± i branches have non-zero intensity, which leads to broadened spectra and less obvious dispersion along l (Fig. 5c, f).
The field-dependence of ω 1 -ω 3 at q = (1/3, 1/3, 1) is extracted from constant-q scans at T = 0.1 K for selected fields 10.5-13.5 T within the plateau (Fig. 6a). By fitting the field-dependence of the low-energy branches of ω 1,2 , which become gapless at a plateau edge, we obtain the quoted model parameters. The field dependence is reproduced fairly well (Fig. 6b, c), although the calculation slightly underestimates ω 3 . We find ω 1 and ω 3 (ω 2 ) increase (decreases) almost linearly in H, while the ω 1 and ω 2 branches cross around 12.6 T. The softening of ω 1 (ω 2 ) at the lower (higher) transition field induces the Y-like (V-like) state, respectively 20,35 . The nonlinearity of the first excitation gap near these transitions (visible only in the calculation) is due to the anisotropy; there is no U(1) symmetry along the field direction for Δ ≠ 1.

Discussion
Our work has mapped out the excitation spectrum in the 1/3 plateau-a manifestation of quantum order-by-disorder Exp. Despite the quantum-mechanical origin of the ground state ordering, we have unambiguously demonstrated the semiclassical nature of magnons in this phase. In fact, the calculated reduction of the sublattice magnetization, δS μ = S − |〈S r 〉| with r ∈ μ, is relatively small (Fig. 3b): δS A e ¼ δS A o ¼ 0:083, δS B e ¼ δS C o ¼ 0:073, and δS C e ¼ δS B o ¼ 0:14 at 10.5 T for the quoted model parameters. This is consistent with the very weak intensity of the two-magnon continuum (Fig. 5g, h). This semiclassical behavior is protected by the excitation gap induced by anharmonicity of the spin waves (magnon-magnon interaction). We note that a perfect collinear magnetic order does not break any continuous symmetry even for Δ = 1, i.e., there is no gapless Nambu-Goldstone mode. The collinearity also means that three-magnon processes are not allowed 44 . The gap is robust against perturbations, such as anisotropies, lattice deformations 40 , or biquadratic couplings for S > 1/2 (a ferroquadrupolar coupling can stabilize the plateau even classically 5 ). Thus, we expect the semiclassical nature of the excitation spectrum to be common to other 2D and quasi-2D realizations of fluctuation-induced plateaus, such as the 1/3 plateau in the spin-5/2 material RbFe(MoO 4 ) 2 29 . Meanwhile, it will be interesting to examine the validity of the semiclassical approach in quasi-1D TLHAFMs, such as Cs 2 CuBr 4 23 , where quantum fluctuations are expected to be stronger.
Finally, we discuss the implications of our results for the zerofield dynamical properties of the same material, where recent experiments revealed unexpected phenomena, such as broadening of the magnon peaks indescribable by conventional spin-wave theory, large intensity of the high-energy continuum 37 , and the extension thereof to anomalously high frequencies 37,39 . Specifically, it was reported that magnon spectral-line broadened throughout the entire Brillouin zone, significantly beyond instrumental resolution, and a high frequency (≳ 2 meV) excitation continuum with an almost comparable spectral weight as single-magnon modes 37 . All of these experimental observations indicate strong quantum effects. Given that the spin Hamiltonian has been reliably determined from our study of the plateau phase, it is interesting to reexamine if a semiclassical treatment of this Hamiltonian can account for the zero-field anomalies.
A semiclassical treatment can only explain the line broadening in terms of magnon decay [44][45][46][47][48] . NLSW theory at H = 0 describes the spin fluctuations around the 120°ordered state by incorporating single-to-two magnon decay at the leading order O(S 0 ). At this order, the two-magnon continuum is evaluated by convoluting LSW frequencies. The self-energies include Hartree-Fock decoupling terms, as well as the bubble Feynman diagrams comprising a pair of cubic vertices with the latter computed with the off-shell treatment. The most crucial one corresponds to the single-to-two magnon decay (see the inset of Fig. 7a), where ω H¼0 k denotes the zero-field magnon dispersion. We show the zero-field dynamical structure factor, S tot H¼0 q; ω ð Þ, at the M point for representative parameters in Fig. 7a-d. The NLSW result for the ideal TLHAFM (J c = 0 and Δ = 1) exhibits sizable broadening and a strong two-magnon continuum [45][46][47][48] (see also Fig. 7e). However, a slight deviation from Δ = 1 renders the decay process ineffective because the kinematic condition, ω H¼0 qÀk , can no longer be fulfilled in 2D for any decay vertex over the entire Brillouin zone if Δ≲0:92 45 . This situation can be inferred from the result for J c = 0 and Δ = 0.85, where the two-magnon continuum is pushed to higher frequencies, detached from the single-magnon peaks. In fact, the sharp magnon lines are free from broadening. The suppression of decay results from gapping out one of the two Nambu-Goldstone modes upon lowering the Hamiltonian symmetry from SU(2) to U(1), which greatly reduces the phase space for magnon decay. The interlayer coupling renders the single-magnon peaks even sharper and the continuum even weaker (Fig. 7c, d).
To determine whether the anomalous zero-field spin dynamics can be explained by a conventional 1/S expansion, it is crucial to estimate Δ very accurately. Previous experiments reported Δ = 0.95 (low-field electron spin resonance experiments compared with LSW theory 33 ) and Δ = 0.89 (zero-field INS experiments compared with NLSW theory 37 ). However, the NLSW calculation reported a large renormalization of the magnon bandwidth (≈ 40% reduction relative to the LSW theory) 37 , suggesting that the previous estimates of Δ may be inaccurate. Particularly, given that Δ is extracted by fitting the induced gap / ffiffiffiffiffiffiffiffiffiffiffi ffi 1 À Δ p , the LSW approximation underestimates 1 − Δ (deviation from the isotropic exchange) because it overestimates the proportionality constant 37 . Figure 7d, f show S tot H¼0 q; ω ð Þ for J c /J = 0.09 and Δ = 0.85. We find that S tot H¼0 q; ω ð Þ remains essentially semiclassical, with sharp magnon lines and a weak continuum, which deviates significantly from the recent results of INS experiments 37,39 . We thus conclude that the Hamiltonian that reproduces the plateau dynamics fails to do so at H = 0 within the spin wave theory, even after taking magnon-magnon interactions into account at the 1/S level. We also mention that the breakdown of the kinematic condition for single-to-two magnon decay also implies the breakdown of the c Comparison of the fitted magnon frequencies ω 1 -ω 3 against the NLSW poles. Error bars represent one standard deviation condition for magnon decay into an arbitrary number of magnons 44 . Thus, the semiclassical picture of weakly interacting magnons is likely inadequate to simultaneously explain the lowenergy dispersions and the intrinsic incoherent features (such as the high-intensity continuum and the line-broadening) observed in Ba 3 CoSb 2 O 9 at H = 0. One may wonder if extrinsic effects can explain these experimental observations. It is possible for exchange disorder to produce continuous excitations as in the effective spin-1/2 triangular antiferromagnet YbMgGaO 4 49 . However, our single crystals are the same high-quality samples reported previously 32,37 , essentially free from Co 2+ -Sb 5+ site-disorder. Indeed, our crystals show only one sharp peak at 3.6 K in the zero-field specific heat 32 in contrast to the previous reports of multiple peaks 31 , which may indicate multi-domain structure. Another possible extrinsic effect is the magnon-phonon coupling, that has been invoked to explain the measured spectrum of the spin-3/2 TLHAFM CuCrO 2 50 . However, if that effect were present at zero field, it should also be present in the UUD state. The fact that Eq. (1) reproduces the measured excitation spectrum of the UUD state suggests that the magnon-phonon coupling is negligibly small (a similar line of reasoning can also be applied to the effect of disorder). Indeed, we have also measured the phonon spectrum of Ba 3 CoSb 2 O 9 in zero field by INS and found no strong signal of magnon-phonon coupling. Our results then suggest that the dynamics of the spin-1/2 TLHAFM is dominated by intrinsic quantum mechanical effects that escape a semiclassical spin-wave description. This situation is analogous to the (π, 0) wave-vector anomaly observed in various spin-1/2 square-lattice Heisenberg antiferromagnets [51][52][53][54][55] , but now extending to the entire Brillouin zone in the triangular lattice. Given recent theoretical success on the square-lattice 56 , our results motivate new non-perturbative studies of the spin-1/2 TLHAFM.

Methods
Neutron scattering measurements. The neutron diffraction data under magnetic fields applied in the [1,-1,0] direction were obtained by using CG-4C cold tripleaxis spectrometer with the neutron energy fixed at 5.0 meV at High Flux Isotope Reactor (HFIR), Oak Ridge National Laboratory (ORNL). The nuclear structure of the crystal was determined at the HB-3A four-circle neutron diffractometer at HFIR, ORNL and then was used to fit the nuclear reflections measured at the CG-4C to confirm that the data reduction is valid. Only the scale factor was refined for fitting the nuclear reflections collected at CG-4C and was also used to scale the moment size for the magnetic structure refinement. 14 magnetic Bragg peaks collected at CG-4C at 10 T were used for the magnetic structure refinement. The UUD spin configuration with the spins along the field direction was found to best fit the data. The nuclear and magnetic structure refinements were carried out using FullProf Suite 57 .
Our inelastic neutron scattering experiments were carried out with the multi axis crystal spectrometer (MACS) 58  Constraint on J due to the saturation field. An exact expression for the saturation field for Hjjĉ, H sat , can be obtained from the level crossing condition between the fully polarized state and the ground state in the single-spin-flip sector. From the corresponding expression, we obtain: where g jj ¼ 3:87 and μ 0 H sat = 32.8 T 33 .
Variational analysis on classical instability of the 1/3 plateau in quasi-2D TLHAFMs. We show that the UUD state is not the classical ground state in the presence of the antiferromagnetic interlayer exchange J c > 0. To verify that the classical ground space for J c = 0 acquires accidental degeneracy in the in-plane where the summation of P Δ is taken over the corner-sharing triangles in each layer, with r = (Δ, μ) (μ = A, B, C) denoting the sublattice sites in each triangle. This simply provides an alternative view of each triangular lattice layer (Fig. 8a).x is the unit vector in the x or field direction. The easy-plane anisotropy forces every spin of the classical ground state to lie in the ab plane and the second term in Eq. (8) has no contribution at this level. Hence, for J c = 0, any three-sublattice spin configuration satisfying S z r ¼ 0 and is a classical ground state, where we momentarily regard S Δ,μ as three-component classical spins of length S. Since there are only two conditions corresponding to the x and y components of Eq. (9), whereas three angular variables are needed to specify the three-sublattice state in the ab plane, the classical ground state manifold for J c = 0 retains an accidental degeneracy, similar to the well-known case of the Heisenberg model (Δ = 1) 43 . The UUD state is the classical ground state only for h red = 3JS. The classical instability of the UUD state for J c > 0 can be demonstrated by a variational analysis. The UUD state in the 3D lattice enforces frustration for onethird of the antiferromagnetic interlayer bonds, inducing large variance of the interlayer interaction. As shown in Fig. 1a, only two of the three spin pairs along the c-axis per magnetic unit cell can be antiferromagnetically aligned, as favored by J c , while the last one has to be aligned ferromagnetically. To seek for a better classical solution, we consider a deformation of the spin configuration parameterized by 0 ≤ θ ≤ π at h red = 3JS, such that the spin structure becomes noncollinear within the ab plane (Fig. 8b). Because the magnetization in each layer is fixed at S/3 per spin, the sum of the energies associated with the intralayer interaction and the Zeeman coupling is unchanged under this deformation. In the meantime, the energy per magnetic unit cell of the interlayer coupling is varied as We find that E c (θ) is minimized at θ = π/3 for J c > 0, corresponding to a saddle point. This is a rather good approximation of the actual classical ground state for small J c > 0, as can be demonstrated by direct minimization of the classical energy obtained from Eq. (1). The crucial observation is that the classical ground state differs from the θ = 0 UUD state.
NLSW calculation for the UUD state. We summarize the derivation of the spin wave spectrum in the quasi-2D TLHAFM with easy-plane anisotropy [see Eq. (1)]. As discussed in the main text, we first work on the 2D limit J c = 0 exactly at h red = 3JS, and a given value of 0 ≤ Δ ≤ 1, which are the conditions for the UUD state to be the classical ground state. Defining the UUD state as shown in Fig. 1a, we introduce the Holstein-Primakoff bosons, a ðyÞ μ;r as in Eqs (2)-(4). After performing a Fourier transformation, a μ;k ¼ 1=N mag 1=2P r 2 μ e ÀikÁr a μ;r , where N mag = N/6 is the number of magnetic unit cells (six spins for each) and N is the total number of spins, we obtain the quadratic Hamiltonian as the sum of even layers (sublattices A e -C e ) and odd layers (sublattices A o -C o ) contributions: where the constant term has been dropped. Here, and matrix notation for the quadratic coefficients The excitation spectrum of H 0 LSW retains two k-linear modes at k = 0 (Fig. 2b). Below, we include nonlinear terms to gap out these excitations. At this stage, the nonlinear terms correspond to the mean-field (MF) decoupling of the intra-layer quartic terms. Once we obtain such a MF Hamiltonian with the gapped spectrum, the deviation from the fine-tuned magnetic field h red = 3JS and interlayer coupling (as well as some other perturbation, if any) can be included. Here, the additional term contains both LSW and NLSW terms. To proceed, we first define the following mean-fields (MFs) symmetrized by using translational and rotational invariance:  Fig. 8 Classical instability of the UUD state in the quasi-2D lattice. a Threesublattice structure for a single layer and a decomposition of the intralayer bonds into corner-sharing triangles. b Deformation of the UUD state (see Fig. 1a) parameterized by θ shown in the projection in the ab (or xy) plane; the spins in sublattices B e and C o are unchanged ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04914-1 Here the MF parameters are given as follows. First, those associated with the intralayer coupling are for even layers and for odd layers. Similarly, the new MF parameters associated with the interlayer coupling are Figure 9 shows the Δ-dependence of these MF parameters. Because these MF parameters are real valued, the coefficient matrix in Eq. (20) has the form where α k α y Àk is the six-component vector comprising the annihilation (creation) operators of Bogoliubov bosons. The transformation matrices satisfy U μκ k ¼ ðU μκ Àk Þ Ã and V μκ k ¼ ðV μκ Àk Þ Ã . The poles, ω κ,k , are the square-roots of the eigenvalues of S 2 (P k + Q k )(P k −Q k ).
When calculating the sublattice magnetization, the reduction of the ordered moment relative to the classical value S corresponds to the local magnon density. With the phase factors for each sublattice, Fig. 1a), we have hS x r i ¼ c μ S À ha y μ;r a μ;r i ¼ c for site r in sublattice μ.
The dynamical spin structure factor is defined by where S α q ¼ N À1=2 P r S α r e ÀiqÁr and |n〉 and ω n denote the nth excited state and its excitation energy, respectively. The longitudinal spin component is with Q = (1/3, 1/3,1) and c μ a y μ;kÀq a μ;k ; ð31Þ We truncate the expansions of the transverse spin components at the lowest order: The transverse components of the dynamical structure factor, S ? q; ω ð Þ¼S yy q; ω ð ÞþS zz q; ω ð Þ, reveal the magnon dispersion,