The nature of spin excitations in the one-third magnetization plateau phase of Ba$_3$CoSb$_2$O$_9$
Y. Kamiya L. Ge Tao Hong Y. Qiu D. L. Quintero-Castro Z. Lu H. B. Cao M. Matsuda E. S. Choi C. D. Batista M. Mourigal H. D. Zhou J. Ma

Reinforced stochastic gradient descent for deep neural network learning
    Taro Toyoizumi
    December 30, 2023
=======================================================================

23pt



  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 Ba_3CoSb_2O_9 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 Ba_3CoSb_2O_9 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.



18pt

Quantum fluctuations favor collinear spin order in frustrated magnets <cit.>, which can be qualitatively different from the classical limit (S →∞) <cit.>. In particular, quantum effects can produce magnetization plateaus <cit.>, where the magnetization is pinned at a fraction of its saturation value. Magnetization plateaus can be interpreted as crystalline states of bosonic particles, and are naturally stabilized by easy-axis exchange anisotropy, which acts as strong off-site repulsion <cit.>. However, the situation is less evident and more intriguing for isotropic Heisenberg magnets, which typically have no plateaus in the classical limit. In a seminal work, Chubukov and Golosov predicted the 1/3 magnetization plateau in the quantum triangular lattice Heisenberg antiferromagnet (TLHAFM), corresponding to an up-up-down (UUD) state <cit.>. Their predictions were confirmed by numerical studies <cit.> and extended to plateaus in other models <cit.>. Experimentally, the 1/3 plateau has been observed in the spin-1/2 isosceles triangular lattice material Cs_2CuBr_4, <cit.> as well as in the equilateral triangular lattice materials RbFe(MoO_4)_2 (S=5/2) <cit.> and Ba_3CoSb_2O_9 (effective S= 1/2) <cit.>.


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 <cit.>. 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 |⟨𝐒_𝐫⟩| close to the full moment. Given that a sizable reduction of |⟨𝐒_𝐫⟩| 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 <cit.> stabilizing the plateau <cit.>, this phenomenon is produced by the zero-point energy correction E_ zp = (1/2)∑_𝐪ω_𝐪 + O(S^0) (ω_𝐪 is the spin wave dispersion), which does not necessarily produce a large moment size reduction.

  
One of our goals is to resolve this seemingly contradictory situation. Recently, Alicea developed a method to fix the unphysical spin-wave instability <cit.>. This proposal awaits experimental verification because the excitation spectrum has not been measured over the entire Brillouin zone for any fluctuation-induced plateau. We will demonstrate that the modified nonlinear spin wave (NLSW) approach indeed reproduces the magnetic excitation spectrum of Ba_3CoSb_2O_9 within the 1/3 plateau <cit.>. The resulting model parameters confirm that the anomalous zero-field dynamics reported in two independent experiments <cit.> must be intrinsic and non-classical.


In this , 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_3CoSb_2O_9. The Hamiltonian is

    = 
      J ∑_⟨𝐫𝐫'⟩( S^x_𝐫 S^x_𝐫' + S^y_𝐫 S^y_𝐫' + Δ S^z_𝐫 S^z_𝐫')
      
       
      + J_c ∑_𝐫( S^x_𝐫 S^x_𝐫+𝐜/2 + S^y_𝐫 S^y_𝐫+𝐜/2 + Δ S^z_𝐫 S^z_𝐫+𝐜/2)
      - ∑_𝐫 S_𝐫^x,

where ⟨𝐫𝐫'⟩ runs over in-plane nearest-neighbor (NN) sites of the stacked triangular lattice and 𝐜/2 corresponds to the interlayer spacing (Fig. <ref>a). 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). = g_⊥ H is the reduced field and g_⊥ is the in-plane g-tensor component.





This model describes Ba_3CoSb_2O_9 (Fig. <ref>b), which comprises triangular layers of effective spin 1/2 moments arising from the 𝒥 = 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 𝐐 = (1/3,1/3,1) <cit.>. Experiments confirmed a 1/3 magnetization plateau for 𝐇∥ ab (Fig. <ref>c) <cit.>, 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_3CoSb_2O_9, along with maps of the excitation spectrum obtained from INS. The excellent agreement between theory and experiment demonstrates the semiclassical nature of magnons within the 1/3 plateau phase, despite the quantum fluctuation-induced nature of the ground state ordering.




§ RESULTS



 §.§ 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_3CoSb_2O_9 <cit.>, 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 <cit.> 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 normally require special conditions to be a classical ground state <cit.>.





To put this into a proper semiclassical framework, we apply Alicea 's trick originally applied to a distorted triangular lattice <cit.>. Basically, we make a “detour” in the parameter space with the additional 1/S-axis quantifying the quantum effect (Fig. <ref>). Instead of expanding the Hamiltonian around S →∞ for the actual model parameters, we start from the special point, J_c = 0, = 3 J S, and a given value of 0 ≤Δ≤ 1, that includes the UUD state in its classical ground state manifold. Assuming the spin structure in Fig. <ref>a, we define

    S^x_𝐫 = S̃^z_𝐫,  
      S^y_𝐫 = S̃^y_𝐫,  
      S^z_𝐫 = -S̃^x_𝐫,

for 𝐫∈ A_e, B_e, A_o, and C_o and

    S^x_𝐫 = -S̃^z_𝐫,  
      S^y_𝐫 = S̃^y_𝐫,  
      S^z_𝐫 = S̃^x_𝐫,

for 𝐫∈ C_e, B_o. Introducing the Holstein-Primakoff bosons, a^(†)_μ,𝐫, 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 have

    S̃^z_𝐫   = S - a^†_μ,𝐫 a^ _μ,𝐫,
      
    S̃^+_𝐫   = S̃^x_𝐫 + i S̃^y_𝐫≈√(2S)( 1 - a^†_μ,𝐫 a^ _μ,𝐫/4S) a^ _μ,𝐫,

and S̃^-_𝐫 = ( S̃^+_𝐫)^† for 𝐫∈μ, truncating higher order terms irrelevant for the quartic interaction. We evaluate magnon self-energies arising from decoupling of the quartic term.





As shown in Fig. <ref>b, the linear spin wave (LSW) spectrum for J_c = 0 and = 3 J S features two q-linear gapless branches at 𝐪 = 0, both of which are gapped out by the magnon-magnon interaction (Fig. <ref>c). Small deviations from J_c = 0 and = 3 J S 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 = 3 J S. Figures <ref>d and <ref>e show the spectra for  shifted by ± 10% from = 3 J S, where we still keep J_c = 0. For < 3 J S, a band-touching and subsequent hybridization appear between the middle and the top bands around 𝐪=(1/6,1/6) (Fig. <ref>d). For > 3 J S, a level-crossing between the middle and bottom bands appears at around 𝐪=0 (Fig. <ref>e). A small J_c > 0 splits the three branches into six (Figs. <ref>f and <ref>g). Figures <ref>a and <ref>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_μ⟩|/S ≲ 30% throughout the local stability range of the plateau. Thus, our semiclassical approach is fully justified within the plateau phase. Figures <ref>c and <ref>d show the field dependence of the staggered magnetization,

    M_UUD = 1/6(⟨S^x_A_e⟩ + ⟨S^x_B_e⟩ - ⟨S^x_C_e⟩ + ⟨S^x_A_o⟩ - ⟨S^x_B_o⟩ + ⟨S^x_C_o⟩),

which is almost field-independent; a slightly enhanced field-independence 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 (Figs. <ref>e and <ref>f).







 §.§ UUD state in Ba_3CoSb_2O_9.

Next we show experimental evidence for the UUD state in Ba_3CoSb_2O_9 by neutron diffraction measurements within the plateau phase for field applied along the [1,-1,0] direction. We used the same single crystals reported in Refs. <cit.>, grown by the traveling-solvent floating-zone technique and characterized by neutron diffraction, magnetic susceptibility, and heat capacity measurements. The space group is P6_3/mmc, with the lattice constants a= b = 5.8562 Åand c = 14.4561 Å. The site-disorder between Co^2+ and Sb^5+ is negligible with the standard deviation of 1%, as reported elsewhere <cit.>. The magnetic and structural properties are consistent with previous reports and confirm high quality of the crystals <cit.>. These crystals were oriented in the (h,h,l) scattering plane. The magnetic Bragg peaks at (1/3, 1/3, 0) and (1/3, 1/3, 1) were measured at T = 1.5 K (Figs. <ref>a and <ref>b). The large intensity at both (1/3, 1/3, 0) and (1/3, 1/3, 1) confirms the UUD state at μ_0 H ≥ 9.8 T <cit.> (Fig. <ref>c). The estimated ordered moment is 1.65(3)  at 10 T and 1.80(9)  at 10.9 T. They correspond to 85(2)% and 93(5)% of the full moment <cit.>, roughly coinciding with the predicted range (Fig. <ref>d). This diffraction pattern can be contrasted with that of the 120^∘ state, characterized by a combination of the large intenisty at (1/3, 1/3, 1) and lack of one at (1/3, 1/3, 0). Our diffraction result is fully consistent with previous nuclear magnetic resonance (NMR) <cit.> and magnetization measurements <cit.>.










 §.§ Excitation spectrum.

We now turn to the dynamical properties in the UUD phase. Figures <ref>a–<ref>c show the INS intensity I(𝐪,ω) ≡ k_ i/k_ f (d^2σ/dΩdE_ f) along high-symmetry directions. The applied magnetic field μ_0 H = 10.5 T is relatively close to the transition field μ_0 H_c1 = 9.8 T <cit.> bordering on the low-field coplanar ordered phase <cit.>, while the temperature T = 0.5 K is low enough compared to T_N≈ 5 K <cit.> for the UUD phase at this magnetic field. The in-plane dispersion shown in Fig. <ref>a comprises a seemingly gapless branch at 𝐪=(1/3, 1/3, -1) (Fig. <ref>c), and two gapped modes centered around 1.6 meV and 2.7 meV. Due to the interlayer coupling, each mode corresponds to two non-degenerate branches. As their splitting is below the instrumental resolution, we simply refer to them as ω_1, ω_2 and ω_3, unless otherwise mentioned (Fig. <ref>). The dispersions along the c-direction are nearly flat, as shown in Figs. <ref>b and <ref>c for q=(1/2, 1/2, l) and q=(1/3, 1/3, l), respectively, reflecting the quasi-2D lattice <cit.>. Among the spin wave modes along 𝐪=(1/2, 1/2, l) and 𝐪=(1/3, 1/3, l) in Figs. <ref>b and <ref>c, ω_1
for 𝐪=(1/2, 1/2, l)
displays a relatively sharp spectral line. 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. <ref>a are roughly captured by the theoretical calculation near the low-field onset of the plateau in Fig. <ref>f (= 2.7 J S ≈ 1.03). 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 <cit.>. To refine the quantitative comparison, we calculate the scattering intensity I_tot(𝐪,ω) ≡ (γ r_0/2)^2 | F(𝐪) |^2 ∑_α (1 - q̂^αq̂^α) g^2_α𝒮^αα(𝐪,ω) where F(𝐪) denotes the magnetic form factor of Co^2+ corrected with the orbital contribution, (γ r_0/2)^2 is a constant, q̂^α = q^α/|𝐪|, and 𝒮^αα(𝐪,ω) 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. <ref>a, transverse spin fluctuations related to single-magnon excitations appear in 𝒮^yy and 𝒮^zz, while longitudinal spin fluctuations corresponding to the two-magnon continuum appear in the inelastic part of 𝒮^xx, denoted as 𝒮_∥. Accordingly, I_tot(𝐪,ω) can be separated into transverse I_⊥ and longitudinal I_∥ contributions. To compare with our experiments, the theoretical intensity is convoluted with momentum binning effects (only for I_⊥) and empirical instrumental energy resolution. Figures <ref>d–<ref>f show the calculated I_⊥(𝐪,ω), along the same high-symmetry paths as the experimental results in Figs. <ref>a–<ref>c, for , , , and . 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 𝐇∥ĉ <cit.>. To obtain the best fit, we also analyzed the field dependence of ω_1, ω_2 and ω_3 (Fig. <ref>).


Remarkably, the calculation in Figs. <ref>d–<ref>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 (Figs. <ref>e and <ref>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. <ref>e), yielding a relatively sharp spectral line (Fig. <ref>h). The corresponding bandwidth ≈ 0.2 meV (Fig. <ref>b) 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 (Figs. <ref>c and <ref>f).


The field-dependence of ω_1–ω_3 at 𝐪= (1/3, 1/3, 1) is extracted from constant-q scans at T = 0.1 K for selected fields 10.5 T–13.5 T within the plateau (Fig. <ref>a). 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 (Figs. <ref>b and <ref>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 <cit.>. 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 effect—in Ba_3CoSb_2O_9. 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 - |⟨𝐒_𝐫⟩| with 𝐫∈μ, is relatively small (Fig. <ref>b): δ 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 (Figs. <ref>g and <ref>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 <cit.>. The gap is robust against perturbations, such as anisotropies, lattice deformations <cit.>, or biquadratic couplings for S > 1/2 (a ferroquadrupolar coupling can stabilize the plateau even classically <cit.>). 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 <cit.>. Meanwhile, it will be interesting to examine the validity of the semiclassical approach in quasi-1D TLHAFMs, such as Cs_2CuBr_4 <cit.>, where quantum fluctuations are expected to be stronger.


Finally, we discuss the implications of our results for the zero-field 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 <cit.>, and the extension thereof to anomalously high frequencies <cit.>. Specifically, Ref. <cit.> reported 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. 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 <cit.>. 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 Γ_3 ∼ O(S^1/2) <cit.>, 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. <ref>a),

    Σ(𝐪,ω)
      = 1/2N∑_𝐤|Γ_3(𝐤, 𝐪 - 𝐤; 𝐪)|^2/ω - ω^H=0_𝐤 - ω^H=0_𝐪 - 𝐤 + i0,

where ω^H=0_𝐤 denotes the zero-field magnon dispersion. We show the zero-field dynamical structure factor, 𝒮_H=0^tot(𝐪,ω), at the M point for representative parameters in Figs. <ref>a–<ref>d. The NLSW result for the ideal TLHAFM (J_c = 0 and Δ = 1) exhibits sizable broadening and a strong two-magnon continuum <cit.> (see also Fig. <ref>e). However, a slight deviation from Δ = 1 renders the decay process ineffective because the kinematic condition, ω^H=0_𝐪 = ω^H=0_𝐤 + ω^H=0_𝐪 - 𝐤, can no longer be fulfilled in 2D for any decay vertex over the entire Brillouin zone if Δ≲ 0.92 <cit.>. 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 (Figs <ref>c and <ref>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 <cit.>) and Δ = 0.89 (zero-field INS experiments compared with NLSW theory <cit.>). However, the NLSW calculation reported a large renormalization of the magnon bandwidth (≈ 40% reduction relative to the LSW theory) <cit.>, suggesting that the previous estimates of Δ may be inaccurate. Particularly, given that Δ is extracted by fitting the induced gap ∝√(1 - Δ), the LSW approximation underestimates 1 - Δ (deviation from the isotropic exchange) because it overestimates the proportionality constant <cit.>.


Figures <ref>d and <ref>f show 𝒮_H=0^tot(𝐪,ω) for and . We find that 𝒮_H=0^tot(𝐪,ω) remains essentially semiclassical, with sharp magnon lines and a weak continuum, which deviates significantly from the recent results of INS experiments <cit.>. 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 condition for magnon decay into an arbitrary number of magnons <cit.>. Thus, the semiclassical picture of weakly interacting magnons is likely inadequate to simultaneously explain the low-energy dispersions and the intrinsic incoherent features (such as the high-intensity continuum and the line-broadening) observed in Ba_3CoSb_2O_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 <cit.>. However, our single crystals are the high-quality samples used in Refs. <cit.>, 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 <cit.> in contrast to the previous reports of multiple peaks <cit.>, 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 <cit.>. However, if that effect were present at zero field, it should also be present in the  UUD state. The fact that Eq. (<ref>) 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_3CoSb_2O_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, <cit.> but now extending to the entire Brillouin zone in the triangular lattice. Given recent theoretical success on the square-lattice <cit.>, 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 triple-axis 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 <cit.>.


Our inelastic neutron scattering experiments were carried out with the Multi Axis Crystal Spectrometer (MACS) <cit.> at NIST Center for Neutron Research (NCNR), NIST, and the cold neutron triple-axis spectrometer (V2-FLEXX) <cit.> at Helmholtz-Zentrum Berlin (HZB). The final energies were fixed at 3 meV and 5 meV on the MACS and 3.0 meV on V2-FLEXX.




 §.§ Constraint on J due to the saturation field.

An exact expression for the saturation field for 𝐇∥ĉ, 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:

    J = 
        g_∥ H_sat S^-1/
        3 + 6 Δ + 2 (1 + Δ) (J_c/J)
      ,

where g_∥ = 3.87 and μ_0 H_sat = 32.8 T <cit.>.




 §.§ 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 magnetic field, we rewrite Eq. (<ref>) as

    =
      J/2∑_( 𝐒_,A + 𝐒_,B + 𝐒_,C - /3J𝐱̂)^2
      
       
      - (1 - Δ) J ∑_⟨𝐫𝐫'⟩ S^z_𝐫 S^z_𝐫'
      + J_c ∑_𝐫( S^x_𝐫 S^x_𝐫+𝐜̂/2 + S^y_𝐫 S^y_𝐫+𝐜̂/2 + Δ S^z_𝐫 S^z_𝐫+𝐜̂/2)  + const.,

where the summation of ∑_ is taken over the corner-sharing triangles in each layer, with 𝐫 = (,μ) (μ=A,B,C) denoting the sublattice sites in each triangle. This simply provides an alternative view of each triangular lattice layer (Fig. <ref>a). 𝐱̂ 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. (<ref>) has no contribution. Hence, for J_c = 0, any three-sublattice spin configuration satisfying S^z_𝐫 = 0 and

    𝐒_,A + 𝐒_,B + 𝐒_,C = /3J𝐱̂,  ∀,

is a classical ground state, where we momentarily regard 𝐒_,μ as three-component classical spins of length S. Since there are only two conditions corresponding to the x and y components of Eq. (<ref>), 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) <cit.>. The UUD state is the classical ground state only for = 3 J S.





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 one third of the antiferromagnetic interlayer bonds, inducing large variance of the interlayer interaction. As shown in Fig. <ref>a, only two of the three spin pairs along the 𝐜 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 = 3 J S, such that the spin structure becomes noncollinear within the ab plane (Fig. <ref>b). 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

    E_c(θ) = 2 J_c S^2 ( cos2θ - 2 cosθ).

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. (<ref>). 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. (<ref>)]. As discussed in the main text, we first work on the 2D limit J_c = 0 exactly at = 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. <ref>a, we introduce the Holstein-Primakoff bosons, a^(†)_μ,𝐫 as in Eqs. (<ref>)–(<ref>). After performing a Fourier transformation,
a^ _μ,𝐤 = (1/)^1/2∑_𝐫∈μ e^-i𝐤·𝐫 a_μ,𝐫,
where = 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 ) and odd layers (sublattices A_o–C_o) contributions:

    =  +  ,

where the constant term has been dropped. Here, 

    = S/2∑_𝐤∈RBZ[  (𝐚^†_𝐤)^T (𝐚^ _-𝐤)^T ][ H^0_11,𝐤 H^0_12,𝐤; H^0_21,𝐤 H^0_22,𝐤 ][  𝐚^ _𝐤; 𝐚^†_-𝐤 ],

with H^0_11,𝐤 = H^0_22,𝐤, H^0_12,𝐤 = H^0_21,𝐤, where the summation over 𝐤 is taken in the reduced Brillouin zone (RBZ) corresponding to the magnetic unit cell of the UUD state. From now on, we will denote this summation as ∑_𝐤. We have introduced vector notation for the operators,

    𝐚^_𝐤 
      =  
      [ a^ _A_e,𝐤; a^ _B_e,𝐤; a^ _C_e,𝐤 ]≡[ a^ _1,𝐤; a^ _2,𝐤; a^ _3,𝐤 ],  𝐚^†_-𝐤 
      = 
      [ a^†_A_e,-𝐤; a^†_B_e,-𝐤; a^†_C_e,-𝐤 ]≡[ a^†_1,-𝐤; a^†_2,-𝐤; a^†_3,-𝐤 ],

and matrix notation for the quadratic coefficients,

    H^0_11,𝐤   =
      [               S^-1  3/2J(1 + Δ) γ^_ k 3/2J(1 - Δ) γ^_- k; 3/2J(1 + Δ) γ^_- k               S^-1  3/2J(1 - Δ) γ^_ k;  3/2J(1 - Δ) γ^_ k 3/2J(1 - Δ) γ^_- k          6J - S^-1 ],
      [1]
      
    
      H^0_12,𝐤   =
      [                   0  -3/2J(1 - Δ) γ^_ k -3/2J(1 + Δ) γ^_- k; -3/2J(1 - Δ) γ^_- k                   0  -3/2J(1 + Δ) γ^_ k;  -3/2J(1 + Δ) γ^_ k -3/2J(1 + Δ) γ^_- k                   0 ],

with γ_𝐤 = 1/3 (e^i𝐤·𝐚 + e^i𝐤·𝐛 + e^-i𝐤·(𝐚 + 𝐛)).
Similarly, we have

    = S/2∑_𝐤[  (𝐚̅^†_𝐤)^T (𝐚̅^ _-𝐤)^T ][ H̅^0_11,𝐤 H̅^0_12,𝐤; H̅^0_21,𝐤 H̅^0_22,𝐤 ][  𝐚̅^ _𝐤; 𝐚̅^†_-𝐤 ],

with

    𝐚̅^_𝐤 
      =  
      [ a^ _A_o,𝐤; a^ _B_o,𝐤; a^ _C_o,𝐤 ]≡[ a^ _4,𝐤; a^ _5,𝐤; a^ _6,𝐤 ],  𝐚̅^†_-𝐤 
      = 
      [ a^†_A_o,-𝐤; a^†_B_o,-𝐤; a^†_C_o,-𝐤 ]≡[ a^†_4,-𝐤; a^†_5,-𝐤; a^†_6,-𝐤 ],

and

    H̅^0_11,𝐤   =
      H̅^0_22,𝐤
      = 
      [ 0 1 0; 0 0 1; 1 0 0 ]
      H^0_11,𝐤[ 0 0 1; 1 0 0; 0 1 0 ],
      
    H̅^0_12,𝐤   =
      H̅^0_21,𝐤
      =
      [ 0 1 0; 0 0 1; 1 0 0 ]
      H^0_12,𝐤[ 0 0 1; 1 0 0; 0 1 0 ].



The excitation spectrum of  retains two relativistic modes at 𝐤 = 0 (Fig. <ref>b). 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 = 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:

    ρ^_μ   = 1/∑_𝐫∈μa^†_μ,𝐫 a^ _μ,𝐫,
      
    δ^_μ   = 1/∑_𝐫∈μ(a^ _μ,𝐫)^2,
      [1]
      
    ξ^_μν   = 1/3∑_𝐫∈μ∑_η̂_μν^a^†_μ,𝐫 a^ _ν,𝐫+η̂_μν^,
      [1]
      
    ζ^_μν   = 1/3∑_𝐫∈μ∑_η̂_μν^a^ _μ,𝐫 a^ _ν,𝐫+η̂_μν^,

where η̂_μν^ represents the in-plane displacement vector connecting sites 𝐫∈μ to a nearest neighbor site in sublattice ν. The mean values ... are evaluated with the ground state of . The MFs for odd (even) layers are obtained  from those for even (odd) layers as

    ρ^_A_o = ρ^_B_e,  ρ^_B_o = ρ^_C_e,  ρ^_C_o = ρ^_A_e,  [1]
      
       δ^_A_o = δ^_B_e,  δ^_B_o = δ^_C_e,  δ^_C_o = δ^_A_e,  [1]
      
       ξ^_A_oB_o = ξ^_B_eC_e,  ξ^_B_oC_o = ξ^_C_eA_e,  ξ^_C_oA_o = ξ^_A_eB_e,  [1]
      
       ζ^_A_oB_o = ζ^_B_eC_e,  ζ^_B_oC_o = ζ^_C_eA_e,  ζ^_C_oA_o = ζ^_A_eB_e.

By collecting all the contributions mentioned above, we obtain the NLSW Hamiltonian,
  
    = S/2∑_𝐤[   (𝐚^†_𝐤)^T  (𝐚̅^†_𝐤)^T  (𝐚^ _-𝐤)^T (𝐚̅^ _-𝐤)^T ][       𝐤^       𝐤^       𝐤'       𝐤';  (𝐤^)_^†       𝐤^ (-𝐤')_^T       𝐤'; (-𝐤')_^∗ (-𝐤')_^∗ (-𝐤^)_^∗ (-𝐤^)_^∗;  (𝐤')_^† (-𝐤')_^∗ (-𝐤^)_^T (-𝐤^)_^∗ ][   𝐚^ _𝐤;  𝐚̅^ _𝐤;  𝐚^†_-𝐤; 𝐚̅^†_-𝐤 ],

  where
  
    𝐤^   
        = H^0_11,𝐤
        +
        [ -2 J_c      0      0;      0  2 J_c      0;      0      0  2 J_c ]
        + S^-1[ μ^A_e_MF + 2 J_c ρ^_A_o     (t^A_eB_e_MF)^∗γ^_𝐤        t^C_eA_e_MFγ^_-𝐤;        t^A_eB_e_MFγ^_-𝐤  μ^B_e_MF -2 J_c ρ^_B_o     (t^B_eC_e_MF)^∗γ^_𝐤;     (t^C_eA_e_MF)^∗γ^_𝐤        t^B_eC_e_MFγ^_-𝐤 μ^C_e_MF - 2 J_c ρ^_C_o ],
        [1]
        
    𝐤^   
        = H̅^0_11,𝐤
        +
        [ -2 J_c      0      0;      0  2 J_c      0;      0      0  2 J_c ]
        + S^-1[ μ^A_o_MF + 2 J_c ρ^_A_e     (t^A_oB_o_MF)^∗γ^_𝐤        t^C_oA_o_MFγ^_-𝐤;        t^A_oB_o_MFγ^_-𝐤 μ^B_o_MF - 2 J_c ρ^_B_e     (t^B_oC_o_MF)^∗γ^_𝐤;     (t^C_oA_o_MF)^∗γ^_𝐤        t^B_oC_o_MFγ^_-𝐤 μ^C_o_MF - 2 J_c ρ^_C_e ],
        [1]
        
    𝐤'
           =
        H^0_12,𝐤
        + S^-1[         Γ^A_e_MF  g^A_eB_e_MFγ^_𝐤 g^C_eA_e_MFγ^_-𝐤; g^A_eB_e_MFγ^_-𝐤         Γ^B_e_MF  g^B_eC_e_MFγ^_𝐤;  g^C_eA_e_MFγ^_𝐤 g^B_eC_e_MFγ^_-𝐤           Γ^C_MF ],
        [1]
        
    𝐤'
           =
        H̅^0_12,𝐤
        + S^-1[         Γ^A_o_MF  g^A_oB_o_MFγ^_𝐤 g^C_oA_o_MFγ^_-𝐤; g^A_oB_o_MFγ^_-𝐤         Γ^B_o_MF  g^B_oC_o_MFγ^_𝐤;  g^C_oA_o_MFγ^_𝐤 g^B_oC_o_MFγ^_-𝐤         Γ^C_o_MF ],
        [1]
        
    𝐤^   =
        cos k_3
        [ J_c (1 + Δ)           0           0;           0 J_c (1 - Δ)           0;           0           0  J_c(1 - Δ) ]
        + S^-1cos k_3 
        [ (t^A_eA_o_MF)^∗               0               0;               0 (t^B_eB_o_MF)^∗               0;               0               0 (t^C_eC_o_MF)^∗ ],
        [1]
        
    𝐤' 
           = 
        cos k_3
        [ -J_c (1 - Δ)            0            0;            0 -J_c (1 + Δ)            0;            0            0  -J_c(1 + Δ) ]
        + S^-1cos k_3
        [ g_MF^A_eA_o           0           0;           0 g_MF^B_eB_o           0;           0           0 g_MF^C_eC_o ].

  Here the MF parameters are given as follows. First, those associated with the intralayer coupling are
  
    μ^A_e_MF   = 3J [
          ρ^_B_e - ρ^_C_e
          - 1 + Δ/2( ξ^_A_eB_e - ζ^_C_eA_e)
          - 1 - Δ/2( ξ^_C_eA_e - ζ^_A_eB_e)
          ],
        [1]
        
    μ^B_e_MF   = 3J [
          ρ^_A_e - ρ^_C_e
          - 1 + Δ/2( ξ^_A_eB_e - ζ^_B_eC_e)
          - 1 - Δ/2( ξ^_B_eC_e - ζ^_A_eB_e)
          ],
        [1]
        
    μ^C_e_MF   = 3J [
          -ρ^_A_e - ρ^_B_e 
          + 1 + Δ/2( ζ^_B_eC_e + ζ^_C_eA_e)
          - 1 - Δ/2( ξ^_B_eC_e + ξ^_C_eA_e)
          ],
        [1]
        
    
        t^A_eB_e_MF   = 3J [ ξ^_A_eB_e
          - 1 + Δ/4( ρ^_A_e + ρ^_B_e)
          + 1 - Δ/8( δ^∗_A_e + δ^_B_e) ],
        [1]
        
    
        t^B_eC_e_MF   = 3J [ -ξ^_B_eC_e
          + 1 + Δ/8( δ^∗_B_e + δ^_C_e)
          - 1 - Δ/4( ρ^_B_e + ρ^_C_e) ],
        [1]
        
    
        t^C_eA_e_MF   = 3J [ -ξ^_C_eA_e
          + 1 + Δ/8( δ^∗_C_e + δ^_A_e)
          - 1 - Δ/4( ρ^_C_e + ρ^_A_e) ],
        [1]
        
    Γ^A_e_MF   = 3J/2[
          1 + Δ/2(ξ^_C_eA_e - ζ^_A_eB_e)
          + 1 - Δ/2(ξ^∗_A_eB_e - ζ^_C_eA_e)
          ],
        [1]
        
    Γ^B_e_MF   = 3J/2[
          1 + Δ/2(ξ^∗_B_eC_e - ζ^_A_eB_e)
          + 1 - Δ/2(ξ^_A_eB_e - ζ^_B_eC_e)
          ],
        [1]
        
    Γ^C_e_MF   = 3J/2[
          1 + Δ/2(ξ^_B_eC_e - ξ^∗_C_eA_e)
          - 1 - Δ/2(ζ^_B_eC_e + ζ^_C_eA_e)
          ],
        [1]
        
    
        g^A_eB_e_MF   = 3J [ ζ^_A_eB_e
          - 1 + Δ/8( δ^_A_e + δ^_B_e)
          + 1 - Δ/4( ρ^_A_e + ρ^_B_e) ],
        [1]
        
    
        g^B_eC_e_MF   = 3J [ -ζ^_B_eC_e
          + 1 + Δ/4( ρ^_B + ρ^_C)
          - 1 - Δ/8( δ^_B + δ^_C)
          ],
        [1]
        
    
        g^C_eA_e_MF   = 3J [ -ζ^_C_eA_e
          + 1 + Δ/4( ρ^_C_e + ρ^_A_e)
          - 1 - Δ/8( δ^_C_e + δ^_A_e)
          ],

for even layers and

    μ_MF^A_o = μ_MF^B_e,  μ_MF^B_o = μ_MF^C_e,  μ_MF^C_o = μ_MF^A_e,  
       
      t_MF^A_oB_o = t_MF^B_eC_e,  
      t_MF^B_oC_o = t_MF^C_eA_e,  
      t_MF^C_oA_o = t_MF^A_eB_e,
      
       Γ_MF^A_o = Γ_MF^B_e,  Γ_MF^B_o = Γ_MF^C_e,  Γ_MF^C_o = Γ_MF^A_e,
      [1]
      
       
      g_MF^A_oB_o = g_MF^B_eC_e,  
      g_MF^B_oC_o = g_MF^C_eA_e,  
      g_MF^C_oA_o = g_MF^A_eB_e,

for odd layers. Similarly, the new MF parameters associated with the interlayer coupling are

    t^A_eA_o_MF   = J_c [
        - 1 + Δ/2( ρ^_A_e + ρ^_A_o)
        + 1 - Δ/4( δ^∗_A_e + δ^_A_o)
        ],
      
    
      t^B_eB_o_MF   = J_c [
        1 + Δ/4( δ^∗_B_e + δ^_B_o)
        - 1 - Δ/2( ρ^_B_e + ρ^_B_o)
        ],
      
    
      t^C_eC_o_MF   = J_c [
        1 + Δ/4( δ^∗_C_e + δ^_C_o)
        - 1 - Δ/2( ρ^_C_e + ρ^_C_o)
        ],
      
    
      g_MF^A_eA_o   =
      J_c [
        - 1 + Δ/4( δ^_A_e + δ^_A_o)
        + 1 - Δ/2( ρ^_A_e + ρ^_A_o)
        ],
      [1]    
      
    
      g_MF^B_eB_o   =
      J_c [
        1 + Δ/2( ρ^_B_e + ρ^_B_o)
        - 1 - Δ/4( δ^_B_e + δ^_B_o)
        ],
      
    
      g_MF^C_eC_o   =
      J_c [
        1 + Δ/2( ρ^_C_e + ρ^_C_o)
        - 1 - Δ/4( δ^_C_e + δ^_C_o)
        ].

Figure <ref> shows the Δ-dependence of these MF parameters.
Because these MF parameters are real valued, the coefficient matrix in Eq. (<ref>) has the form,

    H_NLSW =
      [ P^_𝐤 Q^_𝐤; Q^_𝐤 P^_𝐤 ],

with

    P^_𝐤
      =
      [ 𝐤^ 𝐤^; 𝐤^ 𝐤^ ],
          
      Q^_𝐤
      =
      [ 𝐤' 𝐤'; 𝐤' 𝐤' ].

This form can be diagonalized by a Bogoliubov transformation,

    [   𝐚^ _𝐤;  𝐚̅^ _𝐤;  𝐚^†_-𝐤; 𝐚̅^†_-𝐤 ]
      =
      ([        ; 3c.5U_𝐤 3c.5V_𝐤;        ; 3c.5V_𝐤 3c.5U_𝐤 ])
       ([         ;  .5α^ _𝐤;         ; .5α^†_-𝐤 ]),

where α^ _𝐤 (α^†_-𝐤) is the 6-component vector comprising the annihilation (creation) operators of Bogoliubov bosons. The transformation matrices satisfy U_𝐤^μκ = (U_-𝐤^μκ)^∗ and V_𝐤^μκ = (V_-𝐤^μκ)^∗. The poles, ω^_κ,𝐤, are the square-roots of the eigenvalues of S^2(P_𝐤 + Q_𝐤) (P_𝐤 - Q_𝐤).





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, c^_A_e = c^_B_e = -c^_C_e = c^_A_o = -c^_B_o = c^_C_o = 1 (see Fig. <ref>a), we have

    ⟨S^x_𝐫⟩   = c^_μ( S - ⟨ a^†_μ,𝐫 a^ _μ,𝐫⟩)
      = c^_μ( S - 1/∑_𝐤∑_κ| V_𝐤^μκ|^2 ),

for site 𝐫 in sublattice μ.


The dynamical spin structure factor is defined by

    𝒮^αα(𝐪,ω)
         = ∫_-∞^∞dt/2π e^iω t1/N∑_𝐫,𝐫' e^-i𝐪· (𝐫 - 𝐫')⟨
        S_𝐫^α(t) S_𝐫'^α(0)
      ⟩
       =
      ∑_n δ(ω - ω_n) |⟨0| S_𝐪^α|n⟩|^2,

where
S_𝐪^α = N^-1/2∑_𝐫 S_𝐫^α e^-i𝐪·𝐫
and
|n⟩ and ω_n denote the nth excited state and its excitation energy, respectively. The longitudinal spin component is

    S_𝐪^x   =
      √(N)/3S (δ_𝐪,0 + δ_𝐪,𝐐 + δ_𝐪,-𝐐)
      + δ S_𝐪^x,

with 𝐐 = (1/3,1/3,1) and

    δ S_𝐪^x
      = - √(1/N)∑_μ,𝐤 c^_μ a^†_μ,𝐤 - 𝐪 a^ _μ,𝐤,

We truncate the expansions of the transverse spin components at the lowest order:

    S_𝐪^y   ≈
      -i √(S/12)∑_μ(
      a^ _μ,𝐪 - a^†_μ,-𝐪),
      
    
      S_𝐪^z   ≈√(S/12)∑_μ (-c^_μ)
      (
      a^ _μ,𝐪 + a^†_μ,-𝐪).

The transverse components of the dynamical structure factor, 𝒮_⊥(𝐪,ω) = 𝒮^yy(𝐪,ω) + 𝒮^zz(𝐪,ω), reveal the magnon dispersion,

    𝒮^yy(𝐪,ω)
         =
      ∑_n δ(ω - ω_n) |⟨0| S_𝐪^y|n⟩|^2
      
       ≈S/12∑_κδ(ω - ω^_κ,𝐪) |∑_μ(U_𝐪^μκ - V_𝐪^μκ)
      |^2,
      
    𝒮^zz(𝐪,ω)
         =
      ∑_n δ(ω - ω_n) |⟨0| S_𝐪^z|n⟩|^2
      
       ≈S/12∑_κδ(ω - ω^_κ,𝐪) |∑_μ c^_μ(U_𝐪^μκ + V_𝐪^μκ)
      |^2.

Meanwhile, 𝒮^xx(𝐪,ω) comprises the elastic contribution and the longitudinal fluctuations,

    𝒮_∥(𝐪,ω)
         =
      ∑_n δ(ω - ω_n) |⟨0|δ S_𝐪^x|n⟩|^2,

which can be evaluated by using Wick's theorem.
The result at T = 0 is

    𝒮_∥(𝐪,ω)
         = Θ(ω) 
      N^-1∑_𝐤∑_κ,λ,μ,ν  
      c^_μ c^_νA_μν;κλ(𝐤;𝐪)
       δ(ω - ω^_κ,-𝐤+𝐪 - ω^_λ,𝐤),

where

    A_μν;κλ(𝐤;𝐪)
         = 
      1/2[
        ( U_𝐤-𝐪^μκ)^∗ V_𝐤^μλ
        + ( V_𝐤-𝐪^μκ)^∗ U_𝐤^μλ]
       [
        U_𝐤-𝐪^νκ( V_𝐤^νλ)^∗
        + V_𝐤-𝐪^νκ(U_𝐤^νλ)^∗].





 §.§ Data Availability.

All relevant data are available from the corresponding authors upon reasonable request.


naturemag_noURL




§ END NOTES



 §.§ Acknowledgments

We thank H. Tanaka, A. Chernyshev, and O. Starykh for valuable discussions. J.M. acknowledges the support of the Ministry of Science and Technology of China (2016YFA0300500). Y.K. acknowledges the financial support by JSPS Grants-in-Aid for Scientific Research under Grant No. JP16H02206. The work at Georgia Tech was supported by ORAU's Ralph E. Powe Junior Faculty Enhancement Award (M. Mourigal). H.D.Z. acknowledges support from NSF-DMR-1350002. The work performed in NHMFL was supported by NSF-DMR-1157490 and the State of Florida. We are grateful for the access to the neutron beam time at the neutron facilities at NCNR, BER-II at Helmholtz-Zentrum Berlin and HFIR operated by ORNL. The research at HFIR at ORNL was sponsored by the Scientific User Facilities Division (T.H., H.B.C., and M. Matsuda), Office of Basic Energy Sciences, U.S. DOE. This work utilized MACS supported in part by the National Science Foundation under Agreement No. DMR-1508249.



 §.§ Author contributions

Y.K. and J.M. conceived the project.
H.D.Z. prepared the samples.
T.H., Y.Q., D.L.Q., Z.L., H.B.C., M. Matsuda, L.G., M. Mourigal, and J.M. performed the neutron scattering experiments.
E.S.C. measured the magnetization.
Y.K., L.G., C.D.B., and M. Mourigal performed the NLSW calculations.
Y.K., J.M., M. Mourigal, and C.D.B. wrote the manuscript with comments from all the authors.



 §.§ Competing financial interests

The authors declare no competing financial interests.



 §.§ Correspondence and materials

Correspondence and requests for materials should be addressed to Y.K. (yoshitomo.kamiya@riken.jp) and J.M. (jma3@sjtu.edu.cn).

