Introduction

The computation of dynamical correlation functions at finite temperature is one of the important open problems of modern quantum many-body physics. These functions are not only crucial to test models against different spectroscopic techniques, but are also critical to the development of fast machine learning tools to accelerate and enhance understanding of problems at the forefront of condensed matter physics. For instance, the inelastic neutron scattering (INS) cross-section of quantum magnets is proportional to the dynamical spin structure factor, S(Q, E). A full calculation of this dynamical correlation function is complex because of the exponential numerical cost of computing the exact Hamiltonian eigenstates as a function of the number of spins Ns. To surmount this challenge, classical approximations such as Landau-Lifshitz dynamics (LLD) have been extensively adopted and applied1,2,3,4,5,6 because their numerical cost becomes linear in Ns. Furthermore, while semi-classical approaches7,8,9 are only applicable at the lowest temperatures (because fluctuations around the classical ground state are assumed to be small), the classical approach, based on the LLD, can be implemented at any temperature by sampling initial states via the Metropolis algorithm (classical Monte Carlo simulation) and solving the classical equations of motion. Therefore, the LLD and its generalizations can be used to compute dynamical spin structure factors over the full temperature range. Furthermore, the possibility to determine a Hamiltonian in the more tractable classical limit is important since a classical description is expected to become a good approximation at high enough temperatures. An approach of this type can even be applied to materials that exhibit strong quantum fluctuations at low temperatures, such as spin liquid candidates3.

LLD was originally introduced to describe the precession of the magnetization in a solid10. This dynamics can be derived as a classical limit of quantum spin systems, whose quantum mechanical state becomes a direct product of SU(2) coherent states. The time evolution of this product state is dictated by the LLD equations. It is known, however, that N-level quantum mechanical systems (N = 2S + 1 for spin systems) admit more than one classical limit11,12,13,14,15,16,17,18. As was pointed out in a recent work19, there are large classes of low-entangled quantum magnets, such as materials with strong single-ion anisotropy, for which it is necessary to use a generalized spin dynamics (GSD) which accounts for non-dipolar components of spin states (non-zero expectation values of other multipoles). This GSD is also necessary to describe magnets with significant biquadratic interactions7,20 and those comprising weakly-coupled entangled units, such as dimers21,22, trimers23, and tetramers24,25,26,27. The hypothesis that the present work seeks to test is that these systems are better described by direct products of SU(N) coherent states at any temperature. This implies that the traditional LLD must be extended to encompass these more general cases19.

The main goal of this work is to test the aforementioned hypothesis by modeling the INS cross-section of the effective S = 1 quasi-2D easy-plane antiferromagnet (AFM), Ba2FeSi2O728,29. In Ba2FeSi2O7, a significant single-ion anisotropy (D ~ 1.42 meV) induced by the large tetragonal distortion of FeO4 tetrahedron in conjunction with spin-orbit coupling of Fe2+(3d6) results in an effective low-energy three-level manifold generated by the spin states, \(\left\vert {S}^{z}=0,\pm 1\right\rangle\). The competition with a relatively weak Heisenberg exchange interaction, places the ground state of this material (α=J/D ~ 0.187) near the quantum critical point at αc ~ 0.158 that separates easy-planar AFM order from a quantum paramagnetic (QPM) phase (see Fig. 1). Since the AFM to QPM transition is driven by an enhancement of the local quadrupolar moment at the expense of the magnitude of the local dipolar moment, a proper classical description must allow for the coexistence of local dipolar and quadrupolar fluctuations, leading to transverse and longitudinal collective modes28. SU(3) coherent states fulfill this condition because an SU(3) spin has 8 = 3 + 5 components that include the three components of the dipole moment and the five components of the quadupolar moment (trace-less symmetric tensor)7,19,30,31,32.

Fig. 1: Schematic phase diagram of a weakly coupled S = 1 easy-plane antiferromagnet.
figure 1

A quantum phase transition occurs between the QPM and AFM phases as a function of α = J/D. The Bloch spheres depict the mean-field spin state: the arrow and ellipsoid indicate the dipole moment and its fluctuations. As α → 0, the spin is described by a three-level SU(3) representation with \(\left\vert {S}^{z}=0,\pm 1\right\rangle\). The \(\left\vert {S}^{z}=0\right\rangle\) level of the QPM (no net dipole moment) is represented by a small disk-shape. As the temperature increases, the \(\left\vert {S}^{z}=\pm 1\right\rangle\) levels (green sphere) become occupied and the system crosses-over into an SU(3) classical paramagnet. Near the QCP, the easy-plane AFM ordering is described by a mean-field product of SU(3) coherent states with small dipole moment and large quadrupolar component (ellipsoid). Easy-plane antiferromagnet Ba2FeSi2O7 (TN = 5.2 K) is located near the QCP (yellow star). For α ≫ 1, the system is better approximated by SU(2) coherent states.

Full size image

In this article, we use INS to explore the temperature dependence of S(Q, E) for Ba2FeSi2O7. The spin excitation spectra exhibit significant transfer of spectral weight across the Néel temperature (TN = 5.2 K) with a continuous evolution from well-defined SU(3) spin-wave modes for T < TN to a diffusive resonant excitation for T > TN, indicating temperature dependent spin dynamics. To account for this behavior, we generalize LLD to SU(3) coherent states. This GSD combined with Monte Carlo simulation for temperature dependent classical SU(3) coherent spin states19 describes the observed temperature dependence of S(Q, E) both above and below TN. This result verifies the main hypothesis of this work, namely that many anisotropic magnets such as Ba2FeSi2O7 can be described by direct products of SU(N) coherent states at any temperature.

Results

Temperature dependent inelastic neutron scattering spectra

Figure 2a shows the temperature dependence of the unpolarized neutron cross-section I(Q,E) of Ba2FeSi2O7, which is proportional to S(Q, E). Below TN, the spectrum exhibits sharp spin-waves corresponding to acoustic (T1) and optical (T2) transverse modes. The longitudinal (L) mode is observed as a broad continuum above the T1-mode throughout the entire Brillouin zone (BZ) due to decay into a pair of transverse modes as described in ref. 28. While the sharp spin-waves disappear above TN, a broad dispersion with a finite gap emerges at the magnetic zone center (ZC), Qm = (1, 0, 0.5). With increasing temperature, the gap size increases and the bandwidth becomes narrower.

Fig. 2: Temperature dependent INS spectra of Ba2FeSi2O7.
figure 2

a Contour maps of the INS data as a function of energy and momentum transfers along [H, 0, L] (L = 0.5, 0.58) measured at T = 1.6 K, 6 K, 10 K, and 40 K, corresponding to the scaled temperatures T/TN = 0.31, 1.15, 1.92, and 7.69 where TN = 5.2 K. Labels T1, T2, and L indicate the acoustic and optical transverse modes, and the longitudinal mode in the spectra, respectively. b Resolution convoluted INS intensities calculated by the GSD method. As described in more detail in the main text, the calculations are performed for the temperature ratio, \({T}^{{{{\rm{cl}}}}}/{T}_{{{{\rm{N}}}}}^{{{{\rm{cl}}}}}\), indicated in each panel. c Polarized neutron scattering at Q = (1, 0, 0) measured at 10 K (T > TN). The non-spin-flip (spin-flip) channel for the neutron spin polarization (∥(0, 1, 0)) response corresponds to the Sxx + Syy (Szz) component of I(Q,E). The green-shaded region indicates the magnetic INS response extracted by fitting a DLDHO function to data (see text). d The green (red) shaded region indicates the Sxx + Syy (Szz) component of I(Q,E) calculated by GSD.

Full size image

To understand the diffusive spectra above TN, we performed a polarized neutron scattering experiment. Figure 2c shows the neutron spin polarization dependence of I(E) at Q = (1, 0, 0). The spin-flip and nonspin-flip scattering cross-sections are coupled to the sample magnetization and the wave-vector, allowing us to extract the directional dependence of S(Q,E). The neutron spin was polarized along [0, 1, 0], which provides separate in-plane (S⊥ = Sxx(Q, E) + Syy(Q, E)) and out-of plane (S∥ = Szz(Q, E)) components of S(Q, E), for nonspin-flip and spin-flip channels, respectively. The nonspin-flip channel is by far the most intense, indicating the diffusive spectra at 10 K mainly comes from the in-plane components of S(Q, E).

Calculation of generalized spin dynamics for SU(3) spin

To account for the measured spectra at finite temperatures, we performed GSD calculations. The low-energy effective Hamiltonian for Ba2FeSi2O7 is \({{{\mathcal{H}}}}={\sum }_{{{{\boldsymbol{r}}}},{{{\boldsymbol{\delta }}}}}{S}_{{{{\boldsymbol{r}}}}}^{\mu }{{{{\mathcal{J}}}}}_{{{{\boldsymbol{\delta }}}}}^{\mu \nu }{S}_{{{{\boldsymbol{r}}}}+{{{\boldsymbol{\delta }}}}}^{\nu }+D{\sum }_{{{{\boldsymbol{r}}}}}{({S}_{{{{\boldsymbol{r}}}}}^{z})}^{2}\)28, with the convention of summation over repeated indices μ, ν = {x, y, z} and δ runs over the neighboring bonds with finite exchange interaction. This Hamiltonian can be recast in terms of SU(3) generators \({\hat{O}}_{{{{\boldsymbol{r}}}}}^{\eta }\)19

$$\hat{{{{\mathcal{H}}}}}=\frac{1}{2}\mathop{\sum}\limits_{{{{\boldsymbol{r}}}},{{{\boldsymbol{\delta }}}}}{\hat{O}}_{{{{\boldsymbol{r}}}}}^{\eta }{{{{\mathcal{J}}}}}_{{{{\boldsymbol{\delta }}}}}^{\eta \gamma }{\hat{O}}_{{{{\boldsymbol{r}}}}+{{{\boldsymbol{\delta }}}}}^{\gamma }+\frac{D}{\sqrt{3}}\mathop{\sum}\limits_{{{{\boldsymbol{r}}}}}\left({\hat{O}}_{{{{\boldsymbol{r}}}}}^{8}+\frac{2}{\sqrt{3}}\right),$$
(1)

where the exchange tensor \({{{{\mathcal{J}}}}}_{{{{\boldsymbol{\delta }}}}}^{\eta \gamma }={\delta }_{\eta \gamma }{J}_{{{{\boldsymbol{\delta }}}}}({\delta }_{\eta 1}+{\delta }_{\eta 2}+{\delta }_{\eta 3}{{{\Delta }}}_{{{{\boldsymbol{\delta }}}}})\) with \({J}_{{{{\boldsymbol{\delta }}}}}=J(J^{\prime} )\), \({{{\Delta }}}_{{{{\boldsymbol{\delta }}}}}={{\Delta }}({{\Delta }}^{\prime} )\) for nearest-neighbor intralayer (interlayer) bonds and 1 ≤ η, γ ≤ 8. The generators \({\hat{O}}_{{{{\boldsymbol{r}}}}}^{1-3}\) correspond to dipolar operators \(({S}_{{{{\boldsymbol{r}}}}}^{x},{S}_{{{{\boldsymbol{r}}}}}^{y},{S}_{{{{\boldsymbol{r}}}}}^{z})\) and \({\hat{O}}^{4-8}\) are the quadrupolar operators (bilinear traceless forms of the dipolar operators). See ref. 19 for the matrix representations. We note that the single-ion anisotropy becomes an external (quadrupolar) field that is linearly coupled to the SU(3) spins.

After taking the classical limit of Eq. (1) using SU(3) coherent states19, \({\hat{O}}_{{{{\boldsymbol{r}}}}}^{\eta }\to {o}_{{{{\boldsymbol{r}}}}}^{\eta }\equiv \left\langle {Z}_{{{{\boldsymbol{r}}}}}\right\vert {\hat{O}}_{{{{\boldsymbol{r}}}}}^{\eta }\left\vert {Z}_{{{{\boldsymbol{r}}}}}\right\rangle\), we obtain the classical equation of motion (EOM) of the SU(3) spins

$$\frac{d{o}_{{{{\boldsymbol{r}}}}}^{\eta }}{dt}=\mathop{\sum}\limits_{{{{\boldsymbol{\delta }}}}}{{{{\mathcal{J}}}}}_{{{{\boldsymbol{\delta }}}}}^{\gamma \alpha }{f}_{\eta \gamma \lambda }{o}_{{{{\boldsymbol{r}}}}}^{\lambda }{o}_{{{{\boldsymbol{r}}}}+{{{\boldsymbol{\delta }}}}}^{\alpha }+\frac{D}{\sqrt{3}}{f}_{\eta 8\lambda }{o}_{{{{\boldsymbol{r}}}}}^{\lambda },$$
(2)

where fηγλ are the SU(3) structure constants: \([{\hat{O}}_{{{{\boldsymbol{r}}}}}^{\eta },{\hat{O}}_{{{{\boldsymbol{r}}}}}^{\gamma }]=i{f}_{\eta \gamma \lambda }{\hat{O}}_{{{{\boldsymbol{r}}}}}^{\lambda }\). To compute S(Q, E) at finite temperature, the initial conditions of Eq. (2) are sampled with the standard Metropolis-Hastings Monte Carlo (MC) algorithm from the CP2 manifold (classical phase space) of SU(3) coherent states (see Supplementary Note 2 for detailed information). The numerical integration methods for the classical EOM (2) are explained in Supplementary Note 3 and ref. 33. The INS intensity I(Q, E) is obtained from the Fourier transform of the classical dipolar operators \({o}_{{{{\boldsymbol{r}}}}}^{\mu }(t)\,\mu =1,2,3\) (see Supplementary Note 2). For the calculation, we used J = 0.266 meV, \(J^{\prime} =0.1J^{\prime}\), and D = 1.42 meV from ref. 28 and finite lattices consisting of 24 × 24 × 12 sites. Additionally, for low values of the spin S, the Néel temperature of the classical spin Hamiltonian is significantly lower than the Néel temperature of the quantum mechanical Hamiltonian because quantum fluctuations further increase the energy gain of the ordered state relative to a disordered state. This is a well-known fact for the traditional classical limit based on SU(2) coherent states which remains true for the more general case that we are considering here. Hence the classical approximation used here underestimates the value of Néel temperature, \({T}_{{{{\rm{N}}}}}^{{{{\rm{cl}}}}}=1.38\)K, compared to experimental value TN by a factor of ~ 3.75. Therefore in Fig. 2a, b we compare the measured and calculated spectra at the same values of T/TN and \({T}^{{{{\rm{cl}}}}}/{T}_{{{{\rm{N}}}}}^{{{{\rm{cl}}}}}\), respectively.

Comparison of measured and calculated spectra

Below TN, the calculated spectrum exhibits T1-, T2- and L-modes (see Fig. 2), where the calculated intensities by the GSD are multiplied by the classical to quantum correspondence factor for a harmonic oscillator \({\beta }^{{{{\rm{cl}}}}}E/(1-{e}^{-{\beta }^{{{{\rm{cl}}}}}E})\) with βcl = 1/kBTcl 34. Since the GSD calculation at low-temperatures coincides with the generalized linear spin-wave calculation28, the decay and renormalization of the T2 and L-modes observed at 1.6 K (Fig. 2a) are not captured by this classical approximation. Capturing these feature requires the nonlinear approach described in ref. 28. Above TN, the GSD calculation reproduces the gapped nature of the spectrum representing a resonant excitation between \(\left\vert {S}^{z}=0\right\rangle\) and \(\left\vert \pm 1\right\rangle\) states with a finite dispersion due to the exchange interaction. In the classical description, this diffusive mode originates from the combined effect of the “external SU(3) field” D that induces a precession of each SU(3) moment with frequency D/ℏ (center of the peak) and the random molecular field due to the exchange interaction with the fluctuating neighboring moments that determines the width of the peak. When J ≪ T, the spectrum thus becomes a dispersion-less broad peak centered around an energy Δpara ≃ D. The computed spectra reproduces the main characteristics of the observed dispersions and bandwidth. Since the \(\left\vert {S}^{z}=0\right\rangle\) and \(\left\vert \pm 1\right\rangle\) states are connected by the components that are transverse to the z-axis, S± = Sx ± iSy, the corresponding intensity of S(Q, E) should appear in the channel S⊥ = Sxx(Q, E)+Syy(Q, E) (see Fig. 2d), which is qualitatively in good agreement with the polarization dependence of the measured S(Q, E).

The detailed spectral change across TN is shown in Fig. 3 which compares the measured and calculated constant momentum scans at the ZC with varying temperature. For T < TN, we consider a wave vector Q = (1, 0, 0.2) that is close but not exactly equal to the ZC Q = (1, 0, 0.5) in order to avoid the large tail of the magnetic Bragg reflection as well as the technical challenges associated with calculating the spectrum at the ZC. In this case, the T1 (Goldstone) mode becomes visible because of its finite energy at Q = (1, 0, 0.2) due to the non-zero [0, 0, L]-dispersion produced by the small inter-layer coupling Jinter. As a result, the three T1, T2, L modes are observed in the spectrum (see Fig. 3a). While the T1 and T2 transverse modes remain nearly unchanged with increasing temperature, the energy of the L-mode decreases and the mode becomes broader and indistinguishable from the quasielastic scattering near TN. Above TN, the quasielastic scattering continuously evolves into a broad peak centered at finite energy (Δpara), whose energy increases gradually with the temperature (see Fig. 3c). To extract the spectral weight of the resonant excitation above TN and the L-mode below TN, the data were fitted with a double Lorentzian function associated with a damped harmonic-oscillator (DLDHO),

$$S(E)=A(n(E)+1)\left[\frac{{{\Gamma }}}{{(E-{{\Delta }})}^{2}+{{{\Gamma }}}^{2}}-\frac{{{\Gamma }}}{{(E+{{\Delta }})}^{2}+{{{\Gamma }}}^{2}}\right],$$
(3)

that provides a simplified description of the contribution of an over-damped mode35,36. The n(E) + 1 is the Bose factor, and A, Δ, and Γ indicate the intensity, energy, and line-width of the peak, respectively. The extracted spectral weights and parameters of the L-mode and the resonant excitations are indicated by the shaded regions in Fig. 3a–c, and summarized in Fig. 3d, e, respectively.

Fig. 3: Detailed constant momentum scans of INS spectra.
figure 3

a Measured constant momentum scans at Q = (1, 0, 0.2) near the ZC with temperatures below TN and b corresponding calculated spectra by GSD method. The transverse and longitudinal modes are indicated with labels `T1', `T2', and `L'. The spectral weight for the L-mode was fitted with a DLDHO function, and the results are indicated by the orange (blue) shaded regions for experimental (calculated) spectra. The calculated intensities by the GSD are shown in b. c Comparison of the measured and calculated momentum scans at the ZC above TN. As in a, b, the spectral weights were determined by fitting with a DLDHO function and the results are indicated by the shaded regions. The extracted spectral shapes of the resonant modes are compared with the GSD calculations shown by the blue solid line. d, e Comparison of contour plots of the constant momentum scans between INS data and GSD calculation across TN and \({T}_{{{{\rm{N}}}}}^{{{{\rm{cl}}}}}\). The position (Δ) and line-width (Γ) of the L-mode and resonant excitation were quantified by the DLDHO function, and are exhibited as points with error-bars, respectively.

Full size image

Figure 3b, c show a comparison of the GSD calculation with the INS data across TN. Remarkably, the spectral weight for the L-mode is enhanced and shifts to low-energy with increasing temperature, which is consistent with the data. Above TN, the GSD calculation gives a diffusive resonant peak-shaped spectrum centered at the energy Δpara that approaches D for T ≫ TN. We note that the traditional LLD based on SU(2) coherent states cannot explain this gapped diffusive mode as well as the L-mode, leading to incorrect results in the high-temperature limit19. However, we notice that the calculated spectrum underestimates the width of the mode at temperatures T ≳ TN. This discrepancy arises from an inadequate normalization of the SU(3) spins at T ≳ TN19, similar to the issue raised in traditional SU(2) LLD37.

This discrepancy in the line-width can be removed in the high-temperature limit by applying an adequate renormalization of the SU(3) spins \({o}_{{{{\boldsymbol{r}}}}}^{\eta }\to \kappa {o}_{{{{\boldsymbol{r}}}}}^{\eta }\), with κ = 2 in the high-temperature (T ≫ TN) limit, as described in Ref. 19. This renormalization guarantees that the SU(3) S(Q, E) satisfies the exact sum-rule in the high-T limit. Properly renormalizing the spin has the additional virtue of bringing the theoretical Néel temperature to \({T}_{{{{\rm{N}}}}}^{{{{\rm{cl}}}}}\,\simeq 7.5\)K (κ = 2) closer to the experimental value TN = 5.2 K. Figure 4a shows the comparison between the new calculations including the renormalization factor and the data measured at the same temperature Tcl = T (experiment). This comparison reveals a better agreement in spectral shape at TN ≪ T than those without renormalization factor (see Fig. 3c). As expected, significant deviations in low-energy spectrum are observed at T = 10 K and 12 K because they are relatively close to the \({T}_{{{{\rm{N}}}}}^{{{{\rm{cl}}}}}\) which deviates from the experimental TN. In other words, since the energy Δpara of the diffusive peak goes to zero at TN, its position is shifted to the lower energy relative to the measured peak at T = 10 K. The emergence of a quasielastic peak (centered at E = 0) in the theoretical calculation indicates proximity to the \({T}_{{{{\rm{N}}}}}^{{{{\rm{cl}}}}}\). Figure 4b, c show the computed INS cross-section, I(Q, E), at T = 10 K and 40 K along the same direction in reciprocal space that is presented in Fig. 2a. The T = 10 K and 40 K spectra with renormalized spin provide a good description of the measured spectra at these temperatures.

Fig. 4: Simulated INS spectra with the renormalized SU(3) spins.
figure 4

a Comparison of the measured constant momentum scans (black dots) at the ZC above TN with the results obtained from GSD simulations (green solid lines) with renormalized SU(3) spins (\({o}_{{{{\boldsymbol{r}}}}}^{\eta }\to \kappa {o}_{{{{\boldsymbol{r}}}}}^{\eta }\) with κ = 2), where the recalculated theoretical Néel temperature \({T}_{{{{\rm{N}}}}}^{{{{\rm{cl}}}}}\) for κ = 2 is \({T}_{{{{\rm{N}}}}}^{{{{\rm{cl}}}}}=7.5\) K (see text). b, c Simulated I(Q, E) at T = 10 K and at T = 40 K with the renormalized SU(3) spins, respectively.

Full size image

Discussion

In summary, the GSD based on SU(3) spin provides a good approximation of the measured INS cross-section over a broad temperature range. The most quantitative deviations are observed at very low-temperatures T ≪ TN and close to TN. The former case is due to the requirement of a one-loop quantum correction to account for the decay of the L-mode28, and the latter is due to the expected discrepancy between the experimental and re-scaled values of TN originated from the renormalization factor κ = 2 to the classical SU(3) spins. This renormalization factor arises from enforcing the sum-rule in the infinite T-limit. Similarly, the GSD of unrenormalized classical SU(3) (κ = 1) leads to the correct sum-rule in the zero temperature limit after quantizing the normal modes. Therefore, the correct scaling factor should be defined as a function κ(T) that monotonically interpolates between the two limiting cases κ(0) = 1 and κ(∞) = 2.

The verification of the main hypothesis of this work has very important consequences for the characterization of quantum magnets. For instance, while INS is an ideally suited technique for extracting models from data, the solution of this “inverse scattering problem" requires the development of fast solvers of the direct problem (inferring the INS cross-section of a given model). A crucial advantage of the GSD demonstrated here is that the cost of the simulations scales linearly in the system size, while the computation cost of exact dynamics grows exponentially33, making it an ideal solver for attacking the inverse scattering problem with machine-learning-based approaches4. Moreover, since the GSD can reproduce the INS data in the high-temperature regime, the method can still be applied to quantum magnets that exhibit long-range entanglement at low enough temperatures, but undergo a quantum to classical crossover above a certain temperature TQC3,38. Finally, we note that the SU(N) approach described here is also relevant to the broad class of materials comprising weakly coupled antiferromagnetic magnets including dimers, trimers, or tetramers as well as magnets with strong single-ion anisotropy, where similar effects may be anticipated.

Methods

Inelastic neutron scattering experiment

For the INS experiments, the same single crystal of Ba2FeSi2O7 (mass: 2.13 g) as was used in ref. 28 was aligned on an aluminum plate with an [H, 0, L] horizontal scattering plane. Unpolarized INS data were collected using the cold neutron triple-axis spectrometer (CTAX) at the High Flux Isotope Reactor (HFIR) and the hybrid spectrometer (HYSPEC) at the Spallation Neutron Source (SNS) located at Oak Ridge National Laboratory39. A liquid helium cryostat was used to control temperature for both experiments. At CTAX, the initial neutron energy was selected using a PG (002) monochromator, and the final neutron energy was fixed to Ef = 3.0 meV by a PG (002) analyzer. The horizontal collimation was guide − open\(-40^{\prime} -120^{\prime}\), which provides an energy resolution with full width half maximum (FWHM) = 0.1 and 0.18 meV for E = 0 and 2.5 meV, respectively. For the HYSPEC experiment, Ei = 9 meV and a Fermi chopper frequency of 300 Hz were used, which gives FWHM = 0.28 meV and 0.19 meV of energy resolution at E = 0 and 2.5 meV, respectively. Measurements were performed by rotating the sample from − 50° to 170° with 1° steps. Data were integrated over K = [ − 0.16, 0.16] and L = [L − 0.1, L + 0.1], and symmetrized over positive and negative H.

We also performed polarized neutron scattering measurement as part of the HYSPEC experiment using XYZ-polarization analysis which is the same configuration as the experiment in ref. 40. The X-axis is defined along Q = [1, 0, 0] for the scattering wave vector, the Z-axis is defined along [0, 1, 0] perpendicular to the scattering plane, and the Y-axis is defined along the direction [0, 0, 1] perpendicular to the X- and Z-axes. In the experiment, the neutron was polarized along the Z-direction, and the nonspin-flip and spin-flip scattering cross-sections provide \({I}_{{{{\rm{n}}}}}(E)+{I}_{{{{\rm{mag}}}}}^{{{{\rm{Z}}}}}(E)\) and \({I}_{{{{\rm{mag}}}}}^{{{{\rm{Y}}}}}(E)\), respectively40,41. The In(E) indicates non-magnetic structure factor and \({I}_{{{{\rm{mag}}}}}^{\alpha }(E)\), where α ∈ {X, Y, Z}, is the component dependent magnetic structure factor. The measured spin-flip and nonspin-flip cross-sections provide distinct \({I}_{{{{\rm{mag}}}}}^{{{{\rm{Z}}}}}(E)\) and \({I}_{{{{\rm{mag}}}}}^{{{{\rm{Y}}}}}(E)\), which correspond to Sxx+yy(E) and Szz(E), respectively, on the crystal axis where x∥a, y∥b, and z∥c in the tetragonal crystal structure. All of data sets were reduced and analyzed using the MANTID42 and DAVE43 software packages.