Decay and renormalization of a longitudinal mode in a quasi-two-dimensional antiferromagnet

Abstract

An ongoing challenge in the study of quantum materials, is to reveal and explain collective quantum effects in spin systems where interactions between different modes types are important. Here we approach this problem through a combined experimental and theoretical study of interacting transverse and longitudinal modes in an easy-plane quantum magnet near a continuous quantum phase transition. Our inelastic neutron scattering measurements of Ba₂FeSi₂O₇ reveal the emergence, decay, and renormalization of a longitudinal mode throughout the Brillouin zone. The decay of the longitudinal mode is particularly pronounced at the zone center. To account for the many-body effects of the interacting low-energy modes in anisotropic magnets, we generalize the standard spin-wave theory. The measured mode decay and renormalization is reproduced by including all one-loop corrections. The theoretical framework developed here is broadly applicable to quantum magnets with more than one type of low energy mode.

Introduction

One of the strongest signatures of collective quantum behavior is the spontaneous quasi-particle decay in interacting bosonic systems, as observed in superfluids^1,2,3 and quantum magnets^4,5,6,7,8. In the latter case, spontaneous magnon decay has been studied in a growing number of lattice geometries and model systems where large quantum fluctuations enhance this many-body effect^9,10. A key finding of these studies is that the strong decay process is accompanied by a significant renormalization of the overall spectrum^{11,12,13,14,15,16}. This spectral renormalization leads to measurable effects in the thermal dynamic and transport properties¹⁷, which are inexplicable without considering the renormalization of the quasi-particle mass. At the same time, the renormalization of the spectra opens an avenue to understand quantum systems since the renormalized single-magnon dispersion provides a stringent test for theories that attempt to describe magnon decay. In other words, approaches that do not fully incorporate these many-body effects will not yield correct values of the interaction parameters extracted from experimental studies.

An important question is how to understand quasi-particle decay in quantum magnets when there is more than one type of low-energy mode, i.e. when the parent particles are not of the same type as the daughter particles. Anisotropic magnets with spin $S\ge 1$ provide a common example of this situation. The additional fluctuations (quadrupolar for $S\ge 1$, octupolar for $S\ge 3/2$, etc.) can generate modes which are not captured by standard SU(2) approaches at the linear level. Rather, the physics is more conveniently described in terms of generalized SU(N) spin-wave theory, where the low-energy modes are described by $N-1$ distinct bosons¹⁸. For example, anisotropic $S=1$ systems where both transverse and longitudinal modes are expected, have been previously treated by linear SU(3) theories^{17,18,19,20,21,22,23}. While linear SU(N) approaches to capture the correct number of low-energy modes, they are unable to reproduce the quasi-particle decay and renormalization generated by the interaction between these modes. To capture these effects requires going beyond the linear level and thus an objective of this paper is to generalize the 1/S-expansion of the SU(2) treatment to SU(3) in order to account for the quasi-particle decay and renormalization produced by the interaction (nonlinear) terms using the quintessential example of interacting longitudinal and transverse modes for an $S=1$ easy-plane quantum magnet as a test case.

In easy-plane quantum magnets, phase transitions can be driven by either fluctuations of the phase or the amplitude of the order parameter²⁴. The phase fluctuations are the transverse modes of the order parameter (Goldstone modes in the long-wavelength limit), whereas amplitude fluctuations correspond to the longitudinal modes. Due to the gapless nature of the Goldstone transverse modes, the longitudinal or “Higgs” mode is kinematically allowed to decay into two transverse modes. This decay becomes more significant in low-dimensional systems. Indeed, the longitudinal mode in two-dimensional (2D) antiferromagnets was originally assumed to be overdamped due to an infrared divergence of the imaginary part of the longitudinal susceptibility^25,26. However, more recent theoretical work predicted that the longitudinal peak should remain visible even in 2D^{27,28,29,30,31,32,33}. One aspect of this problem, which has not been emphasized in previous works, is that the rather strong decay of the longitudinal mode is accompanied by a significant renormalization of the gap and the dispersion of the modes. As noted above, this additional many-body effect provides a hard test for theories that attempt to reproduce the measured decay of the Higgs mode.

As a starting point to understand the physics described above, we focus on the quasi-2D Heisenberg square lattice with effective $S=1$ with an antiferromagnetic exchange coupling ($\widetilde{J}$) and a strong easy-plane single-ion anisotropy ($\widetilde{D}$). In this case, $\alpha =\widetilde{J}/\widetilde{D}$ can be viewed as a tuning parameter that can be used to drive a system from a quantum paramagnet (QPM) to an antiferromagnet (AFM) with an intervening QCP as shown in Fig. 1. Near the QCP, spontaneous symmetry breaking produces two transverse modes (one of them is a Goldstone mode) and a longitudinal Higgs mode. The longitudinal mode is unstable with respect to decay into a pair of transverse modes resulting in an intrinsic line broadening^9,34.

Fig. 1: Schematic diagrams near the quantum critical point.

The schematic phase diagram illustrates the O(2) quantum critical point (QCP) between the antiferromagnetic (AFM) state and the quantum paramagnet (QPM) as a function of ${{{{{\rm{\alpha }}}}}}=\widetilde{J}/\widetilde{D}$ ($\widetilde{J}$ is a Heisenberg exchange and $\widetilde{D}$ is a easy-plane single-ion anisotropy of effective $S$ = 1). The low-energy excitations of the QPM are two degenerate ${S}^{z}$=${\pm}\!1$ modes (black line) with a gap, $\Delta$, which closes at the QCP. The spontaneous U(1) symmetry breaking leads to a gapless magnon or transverse mode ($T$-mode), indicated with a blue line, which is accompanied by a gapped longitudinal mode ($L$-mode) indicated with the orange line. Near the QCP, the energy and the lifetime of the $L$-mode are strongly renormalized (dashed orange line) due to the decay into the continuum of two transverse modes (shaded orange region).

In this paper, we use inelastic neutron scattering to study the spin excitation spectrum of Ba₂FeSiO₇. The high-quality neutron-scattering data reveals a complex spectrum where transverse modes are resolution limited, whereas a longitudinal mode displays significant ${{{{{\bf{Q}}}}}}$-dependent broadening throughout the Brillouin zone (BZ), demonstrating the importance of quasi-particle decay even away from the long-wavelength limit. The neutron-scattering results further show that the longitudinal mode has a very small gap clearly demonstrating that Ba₂FeSiO₇ is relatively close to a QCP. To understand the inelastic neutron-scattering data, we implement a generalized SU(3) spin-wave calculation^17,18,22 and compute the low-energy excitation spectrum of an effective low-energy spin $S=1$ model. After demonstrating that the generalization of the well-known $1/S$-expansion of the SU(2) spin-wave theory^{35,36,37,38,39,40,41} is simply a loop expansion⁴² of the SU(3) spin-wave theory, we show that the one-loop correction is enough to account for the broadening of the longitudinal mode and the large renormalization of the gap and the dispersion of this mode. We further show that not including the one-loop corrections results in Hamiltonian parameters that place the exact ground state of the spin Hamiltonian for Ba₂FeSi₂O₇ on the nonmagnetic side of the QCP—contrary to experimental fact. This provides a dramatic demonstration of the importance of including renormalization effects, where the linear spin-wave calculation overestimates the stability range of the magnetically ordered state. The fact that the one-loop correction can simultaneously account for the real and imaginary part of the self-energy of the longitudinal mode, as well as of the renormalization the transverse mode dispersion, confirms that the easy-plane quantum magnet Ba₂FeSi₂O₇ is an ideal platform for studying many-body effects in the proximity of the O(2) QCP.

Results

Model material

Figure 2a illustrates the crystal structure of Ba₂FeSi₂O₇ comprising layers of FeSi₂O₇ separated by Ba atoms. As shown in Fig. 2b, the FeO₄ tetrahedra of the FeSi₂O₇ layer are connected via SiO₄ polyhedra and the two adjacent Fe²⁺ atoms are coupled through the superexchange interaction, J, that is mediated by the two oxygen ligands (red dashed line in Fig. 1b). The resulting square lattice of magnetic moments is vertically stacked along the c axis, leading to a quasi-2D simple tetragonal spin lattice.

Fig. 2: Crystal and magnetic structure of Ba₂FeSi₂O₇.

a Crystal structure of Ba₂FeSi₂O₇. Ba atoms separate layers composed of FeSi₂O₇, rendering a quasi-two-dimensional structure. b In the FeSi₂O₇ layer, FeO₄ tetrahedra are connected via SiO₄ polyhedra, and the adjacent two Fe²⁺ atoms are exchange coupled by two oxygen ligands. The red dashed line indicates the exchange pathway $J$ within two-dimensional square spin network. The interlayer coupling $J^{\prime}$ is found here to be much weaker than $J$. Red arrows indicate the moment direction in the collinear AFM phase as determined in ref. ⁴⁴. The black solid line indicates the chemical unit cell. c ${HK}$-reciprocal space with $L$ = 0.5 in the tetragonal structure ($P\bar{4}{2}_{1}m$). The blue solid line and the black circle indicate the Brillouin zone and zone center, respectively. The coordinates ($H$, $K$, $L$) of the reciprocal lattice of the origin lattice are related to $({k}_{x},{k}_{y},{k}_{z})$ of the magnetic lattice formed by the Fe²⁺ atoms through ${k}_{x}=\pi (H-K)$, ${k}_{y}=\pi (H+K)$, and ${k}_{z}=2\pi L$. d Illustration of the spin fluctuation modes. ${T}_{1}$ and ${T}_{2}$ indicate transverse fluctuation in the ${ab}$-plane and out-of the plane, respectively. $L$ indicates longitudinal fluctuation of spin.

A detailed description of the single-ion state of the Fe²⁺ ion is given in Note 1 of the Supplementary Information. The combination of a relatively large spin–orbit coupling (λ ~20 meV) and a dominant tetrahedral crystal field (Δ_Td), splits the free-ion levels, ⁵D (L = 2, S = 2), into several multiplets. The lowest energy $S=2$ multiplet has a significant orbital character due to the finite spin–orbit coupling, that combined with the tetragonal distortion (${\delta }_{{{{{{{\mathrm{Tetra}}}}}}}}$) by large compression of the FeO₄ tetrahedra leads to a rather strong easy-plane single-ion anisotropy^43,44. The five $S=2$ energy levels are then split into a singlet ${S}^{z}=0$ ground state and two excited ${S}^{z}={\pm}\!1$ and ${S}^{z}={\pm}\! 2$ doublets with energies D and 4D, respectively (see Fig. S1a of the Supplementary Information). Because the gap D of the ${S}^{z}={\pm}\!1$ doublet is four times smaller than the gap of the ${S}^{z}={\pm}\!2$ doublet and the dominant superexchange interaction J is smaller than D/4 in Ba₂FeSi₂O₇, the low-energy spectrum is well captured by projecting the $S=2$ spin Hamiltonian into the ${S}^{z}=0$ and ${S}^{z}=\pm 1$ low-energy states.

The resultant $S=1$ effective spin Hamiltonian describes the competition between a QPM ($\widetilde{J}\ll \widetilde{D}$) with each spin of the lattice having dominant ${S}^{z}=0$ character, and a collinear AFM state ($\alpha =\widetilde{J}/\widetilde{D} \; > \; {\alpha }_{c}$) with staggered magnetization in the ab plane (see Fig. 2b). Ba₂FeSi₂O₇ turns out to be on the antiferromagnetic side with a Néel temperature T_N = 5.2 K⁴⁴. Below T_N, the spins order antiferromagnetically with propagation vector ${{{{{{\bf{Q}}}}}}}_{{{{{{\bf{m}}}}}}}$ = (1, 0, 0.5), corresponding to $(\pi ,\pi ,\pi )$ as shown in Fig. 2c. The magnetic moments are highly confined in the ab plane due to easy-plane anisotropy, giving rise to the magnetic structure shown in Fig. 2b. A neutron diffraction study on a powder sample revealed a significantly reduced ordered moment of 2.95 μ_B, which is only 63% of the full moment of 4.36 μ_B (${g}_{{ab}}=2.18$) expected for an $S=2$ spin⁴⁴, suggesting the proximity to the quantum critical point. In addition, as described in further detail below, our analysis confirms that $\alpha =\widetilde{J}/\widetilde{D}\;\sim0.184$ is close to the critical value, ${\alpha }_{{{{{{\rm{c}}}}}}}^{2{{{{{\rm{D}}}}}}}$ = 0.18 and ${\alpha }_{{{{{{\rm{c}}}}}}}^{3{{{{{\rm{D}}}}}}}$ = 0.1 for 2D and 3D, respectively, obtained from quantum Monte Carlo simulations²².

The spin excitations of Ba₂FeSi₂O₇ are generically described by an antiferromagnetic $S=2$ spin Hamiltonian on a simple tetragonal lattice:

$${\mathcal H}= \, J\mathop{\sum}\limits _{{\langle{{{{\rm{r}}}}}},{{{{{\rm{r}}}}}}{^\prime} \rangle}[{S}_{{{{{{\rm{r}}}}}}}^{x}{S}_{{{{{{\rm{r}}}}}}^{\prime} }^{x}+{S}_{{{{{{\rm{r}}}}}}}^{y}{S}_{{{{{{\rm{r}}}}}}^{\prime} }^{y}+\varDelta {S}_{{{{{{\rm{r}}}}}}}^{z}{S}_{{{{{{\rm{r}}}}}}^{\prime} }^{z}]\\ +\,J^{\prime} \mathop{\sum}\limits _{{\langle\langle{{{{\rm{r}}}}}},{{{{{\rm{r}}}}}}^{\prime} \rangle\rangle}[{S}_{{{{{{\rm{r}}}}}}}^{x}{S}_{{{{{{\rm{r}}}}}}^{\prime} }^{x}+{S}_{{{{{{\rm{r}}}}}}}^{y}{S}_{{{{{{\rm{r}}}}}}^{\prime} }^{y}+\varDelta^{\prime} {S}_{{{{{{\rm{r}}}}}}}^{z}{S}_{{{{{{\rm{r}}}}}}{^\prime} }^{z}]\\ +\,D\mathop{\sum}\limits _{{{{{{\rm{r}}}}}}}{({S}_{{{{{{\rm{r}}}}}}}^{z})}^{2}.$$

(1)

The bracket 〈r,r′〉(〈〈r,r′〉〉) indicates that the sum runs over intralayer (interlayer) nearest-neighbor spins with isotropic superexchange interaction $J(J^{\prime} )$. Δ(Δ′) is the intralayer (interlayer) uniaxial anisotropy and the last term represents the easy-plane single-ion anisotropy (D > 0).

In the large D/J limit, the ${S}^{z}={\pm}\! 2$ doublet is separated from the ${S}^{z}={\pm}\!1$ doublet by an energy gap 3D. The low-energy subspace of magnetic excitations can then be further reduced by projecting out the ${S}^{z}={\pm}\!2$ doublet. The reduced low-energy Hamiltonian ${\mathcal H}_{{{{{{\rm{e}}}}}}{ff}}$ results from projecting ${{{{{\mathcal{H}}}}}}$ onto the low-energy subspace S spanned by the triplet of states with ${S}^{z}=0,\pm 1$: ${{{{{{\mathcal{H}}}}}}}_{{{{{{\rm{e}}}}}}{ff}}={{{{{{\mathcal{P}}}}}}}_{S}{{{{{\mathcal{H}}}}}}{{{{{{\mathcal{P}}}}}}}_{S}$. The resulting effective spin $S=1$ Hamiltonian is

$${ {\mathcal H} }_{{{{{{\rm{e}}}}}}ff}= \;\tilde{J}\mathop{\sum}\limits _{\langle {{{{{\bf{r}}}}}},{{{{{\bf{r}}}}}}^\prime\rangle }[{s}_{{{{{{\bf{r}}}}}}}^{x}{s}_{{{{{{\bf{r}}}}}}^\prime}^{x}+{s}_{{{{{{\bf{r}}}}}}}^{y}{s}_{{{{{{\bf{r}}}}}}^\prime}^{y}+\tilde{\Delta }{s}_{{{{{{\bf{r}}}}}}}^{z}{s}_{{{{{{\bf{r}}}}}}^\prime}^{z}]\\ +\,\widetilde{J}^{\prime}\mathop{\sum}\limits _{{\langle\langle{{{{\bf{r}}}}}},{{{{{{\bf{r}}}}}}}^{\prime}\rangle\rangle}[{s}_{{{{{{\bf{r}}}}}}}^{x}{s}_{{{{{{\bf{r}}}}}}^\prime}^{x}+{s}_{{{{{{\bf{r}}}}}}}^{y}{s}_{{{{{{\bf{r}}}}}}{^\prime}}^{y}+\tilde{\Delta }^{\prime} {s}_{{{{{{\bf{r}}}}}}}^{z}{s}_{{{{{{\bf{r}}}}}}{^\prime}}^{z}]\\ +\tilde{D}\mathop{\sum} \limits_{{{{{{\bf{r}}}}}}}{({s}_{{{{{{\bf{r}}}}}}}^{z})}^{2}.$$

(2)

with $\widetilde{J}=3J$, $\widetilde{J}^{\prime} =3J^{\prime}$, $\widetilde{\Delta }=\Delta /3$, $\widetilde{\Delta }^{\prime} =\Delta^ {\prime} /3$, and $\widetilde{D}=D$. As we will see below, this simple effective Hamiltonian can explain not only the in-plane antiferromagnetic ordering observed in Ba₂FeSi₂O₇ (see Fig. 2b), but also the spectra of quasi-particle excitations, including rather strong renormalization effects due to proximity to the QCP.

Inelastic neutron scattering

To investigate the spin excitation spectrum in Ba₂FeSi₂O₇, we performed inelastic neutron scattering using two instruments; the cold neutron triple-axis spectrometer (CTAX) at the High Flux Isotope Reactor, and the time-of-flight hybrid spectrometer (HYSPEC) at the Spallation Neutron Source at Oak Ridge National Laboratory⁴⁵. An overview of the inelastic neutron-scattering results is presented in Fig. 3 through contour maps of the neutron-scattering intensity, $I({{{{{\bf{Q}}}}}},\omega )$, along [H, H, 0.5] and [H, 0, 0.5]. For both spectra, strongly dispersive spin excitations extending up to energy ~2.7 meV are observed. Whereas the dispersion along [0, 0, L] direction is weak with a bandwidth of ~0.5 meV (see Note 4 in the Supplementary Information), which is expected for spin excitations of a quasi-two-dimensional spin system.

Fig. 3: Inelastic neutron scattering of Ba₂FeSi₂O₇.

a Contour map of the inelastic neutron-scattering (INS) data as a function of energy and momentum transfer along the [$H$, 0, 0.5] direction measured at $T$ = 1.6 K ($< {T}_{{{{{{\rm{N}}}}}}}$) using the HYSPEC time-of-flight spectrometer at SNS. d Contour map of the INS data as a function of energy and momentum transfer along [$H$, $H$, 0.5] direction measured at $T$ = 1.4 K ($< {T}_{{{{{{\rm{N}}}}}}}$) using the cold Neutron Triple-Axis spectrometer (CTAX) at HFIR. The instrumental resolutions at energy = 2.5 meV for each instrument are indicated with blue bars along the $y$-axis in a and d. The two transverse modes and the longitudinal mode are labeled with ${T}_{1}$, ${T}_{2}$, and $L$, respectively. b, c, e, f INS intensities calculated by the generalized linear spin-wave theory (GLSWT) and GLSWT plus one-loop corrections (GLSWT+one-loop) with the parameter sets ${{{{{\mathcal{A}}}}}}$ and ${{{{{\mathcal{B}}}}}}$ given in Table 1, respectively. The instrumental resolution of HYSPEC and CTAX was modeled in the calculated spectra using a Lorentzian function.

There are several distinct features in the inelastic neutron-scattering data. An intense spin-wave excitation emanates from the magnetic zone center (ZC), ${{{{{\bf{Q}}}}}}$ = (1, 0, 0.5), which arises due to the in-phase oscillation between Fe²⁺ spins in the plane. We refer to this mode as T₁. Along the [H, 0, 0.5] direction toward the zone boundary (ZB) at ${{{{{\bf{Q}}}}}}$ = (0, 0, 0.5), the T₁-mode reaches its maximum energy of ~2.5 meV. Another weak, but sharp mode, is visible along [H, 0, 0.5] with an energy of 2.5 meV at the ZC. We refer to this mode as T₂. These two modes are expected for a strong easy-plane antiferromagnet, where transverse magnons split into gapless in-plane fluctuations (T₁-mode) and gapped out-of-plane fluctuations (T₂-mode). The finite value of the energy gap of the out-of-plane fluctuation at the ZC is associated with the strength of the easy-plane single-ion anisotropy⁴⁶.

The T₁ and T₂ transverse magnon modes are also observed along the [H, H, 0.5] direction in Fig. 3d. Noticeably, an additional sharp mode is observed at the top of the T₁-mode. This mode is visible along the entire Brillouin zone boundary. We refer to this additional mode as “L”-mode. The L-mode is visible in the spectra along [H, 0, 0.5] as well, however, it exhibits dramatic line broadening near the ZC. To demonstrate more clearly the ${{{{{\bf{Q}}}}}}$-dependence of the modes, Fig. 4 shows cuts at constant momentum transfers for multiple points along [H, 0, 0.5] and [H, H, 0.5]. Two pronounced peaks, corresponding to the T₁- and L-modes, remain sharp along the ZB (Fig. 4b). As already noted, the situation is very different near the ZC where the L-mode is significantly broadened (Fig. 4a). We note that the L-mode remains a broad peak near the ZC, rather than a featureless excitation. To investigate the extent of the broadening effect, Gaussian line shapes for the T₁-, T₂-, and L-modes were fit to the individual cuts in Fig. 4. The line widths obtained from the fits are displayed in Fig. 7a–d. These results reveal that the L-mode is three times broader than the instrumental resolution at the ZC (see Fig. 4a), whereas it has comparable line width to instrumental resolution near the ZB.

Fig. 4: Detailed line-cuts of INS spectra.

a Constant momentum cuts at points along the [$H$, 0, 0.5] direction measured using HYSPEC at SNS, integrated over $H$ = [$H$−0.05, $H$+0.05] at selected $H$, $K$ = [−0.1, 0.1], and $L$ = [0.4, 0.6]. b Constant momentum cuts at points along the [$H$, $H$, 0.5] direction measured using CTAX at HFIR. Blue bars at the bottom of the panels indicate the instrumental resolutions for HYSPEC and CTAX at the proximate energy transfers. The blue and orange shaded regions are the results of fitting Gaussian line shape to transverse (${T}_{1}$, ${T}_{2}$) and longitudinal ($L$) modes, respectively.

Generalized spin waves

In this section, we introduce a generalized SU(3) spin-wave calculation^17,18,22,47, which is required to capture the two low-energy (longitudinal and transverse) modes of Ba₂FeSi₂O₇. Clearly, a linear treatment is not enough to capture the decay of the longitudinal mode into two transverse modes. Consequently, the main aim of this section is to lay the groundwork for introducing the loop expansion⁴² (generalization of the 1/S-expansion^{35,36,37,38,39,40,41}) in the section describing the nonlinear corrections.

To account for the transverse and longitudinal modes revealed by the INS experiment, the usual SU(2) spin-wave theory (SWT) must be generalized to SU(3)¹⁸, by introducing the SU(3) Schwinger boson representation of the spin operators ${S}_{{{{{{\bf{r}}}}}}}^{\nu }={{{{{{\boldsymbol{b}}}}}}}_{{{{{{\bf{r}}}}}}}^{{{\dagger}} }{{{{{{\mathcal{S}}}}}}}^{\nu }{{{{{{\boldsymbol{b}}}}}}}_{{{{{{\bf{r}}}}}}}$, where ${{{{{{\boldsymbol{b}}}}}}}_{{{{{{\bf{r}}}}}}}={({b}_{{{{{{\bf{r}}}}}},+1},{b}_{{{{{{\bf{r}}}}}},-1},{b}_{{{{{{\bf{r}}}}}},0})}^{T}$,

$${{{{{{\mathcal{S}}}}}}}^{x}=\frac{1}{\sqrt{2}}({\lambda }_{4}+{\lambda }_{6}),{{{{{{\mathcal{S}}}}}}}^{y}=\frac{1}{\sqrt{2}}({\lambda }_{5}-{\lambda }_{7}),{{{{{{\mathcal{S}}}}}}}^{z}={\lambda }_{3},$$

(3)

${\lambda }_{i}$ are the Gel–Mann matrices and the Schwinger boson operators satisfy the local constraint

$$\mathop{\sum} \limits_{m=\pm 1,0}{b}_{{{{{{\bf{r}}}}}},m}^{{{\dagger}} }{b}_{{{{{{\bf{r}}}}}},m}=M=1.$$

(4)

We note that the SU(3) Schwinger boson representation of the spin operators should not be confused with the Schwinger boson approximation^36,48,49,50, which is qualitatively different from the semi-classical approach that we describe below. The magnetically ordered state of Ba₂FeSi₂O₇ can be approximated by a product (mean-field) state of normalized SU(3) coherent states

$${{{{{\rm{|}}}}}}{\psi }_{{{{{{\bf{r}}}}}}}\rangle ={{\cos }}\theta {{{{{\rm{|}}}}}}0\rangle +({{\sin }}\theta {{\cos }}\phi {{{{{\rm{|}}}}}}1\rangle +{{\sin }}\theta {{\sin }}\phi {{{{{\rm{|}}}}}}-1\rangle ){e}^{i{{{{{{\bf{Q}}}}}}}_{{{{{{\bf{m}}}}}}}\cdot {{{{{\bf{r}}}}}}},$$

(5)

where ${{{{{{\bf{Q}}}}}}}_{{{{{{\bf{m}}}}}}}=(\pi ,\pi ,\pi )$ ((1, 0, 0.5) in the chemical lattice) is the AFM ordering wave vector. Although a general SU(3) coherent state is parameterized by four independent parameters for degenerate representations⁵¹, the two independent parameters $\theta$ and $\phi$ are enough to describe the collinear order under consideration. The three basis states $\left|m\right\rangle (m=0,{\pm}\! 1)$ are represented by creating a boson with the quantum number $m$ from the vacuum: $|m\rangle ={b}_{{{{{{\bf{r}}}}}},m}^{{{\dagger}} }|\phi \rangle$.

As in the usual spin-wave theory, we introduce an SU(3) transformation that rotates the boson operators, ${\widetilde{{{{{{\boldsymbol{b}}}}}}}}_{{{{{{\bf{r}}}}}}}={U}_{{{{{{\bf{r}}}}}}}{{{{{{\boldsymbol{b}}}}}}}_{{{{{{\bf{r}}}}}}}$, to a local basis that includes the coherent SU(3) state (5) as one of its three elements. This local transformation allows us to align the quantization axis with the direction of the local SU(3) order parameter. The spatial dependence of U_r can be removed by working in a twisted frame, where the original AFM order becomes a FM one. This can be done by rotating the spin reference frame of one of the two sublattices of the tetragonal lattice by an angle π along the z axis: ${s}_{{{{{{\bf{r}}}}}}}^{z}\to {s}_{{{{{{\bf{r}}}}}}}^{z}$, and ${s}_{{{{{{\bf{r}}}}}}}^{x,y}\to -{s}_{{{{{{\bf{r}}}}}}}^{x,y}$. In the new reference frame, the effective Hamiltonian (2) becomes

$${\tilde{ {\mathcal H} }}_{{{{{{\rm{e}}}}}}ff}=\tilde{J}\,\mathop{\sum} \limits_{\langle {{{{{\bf{r}}}}}},{{{{{\bf{r}}}}}}^\prime\rangle ,\nu }\,{a}_{\nu }{s}_{{{{{{\bf{r}}}}}}}^{\nu }{s}_{{{{{{\bf{r}}}}}}{^\prime}}^{\nu }+\tilde{J}^{\prime}\mathop{\sum} \limits_{\langle {{{{{\bf{r}}}}}},{{{{{\bf{r}}}}}}{^\prime}\rangle ,\nu }{b}_{\nu }{s}_{{{{{{\bf{r}}}}}}}^{\nu }{s}_{{{{{{\bf{r}}}}}}{^\prime}}^{\nu }+\tilde{D}\mathop{\sum} \limits_{{{{{{\bf{r}}}}}}}\,{({s}_{{{{{{\bf{r}}}}}}}^{z})}^{2},$$

(6)

with ${a}_{x}={a}_{y}={b}_{x}={b}_{y}=-1$, ${a}_{z}=\widetilde{\Delta }$ and ${b}_{z}=\widetilde{\Delta }^{\prime}$, and the SU(3) transformation reads

$$U=\left(\begin{array}{ccc}-{{\sin }}\phi & {{\cos }}\phi & 0\\ {{\cos }}\theta {{\cos }}\phi & {{\cos }}\theta {{\sin }}\phi & -{{\sin }}\theta \\ {{\sin }}\theta {{\cos }}\phi & {{\sin }}\theta {{\sin }}\phi & {{\cos }}\theta \end{array}\right).$$

(7)

The bosonic representation of ${\widetilde{{{{{{\mathcal{H}}}}}}}}_{{{{{{\rm{e}}}}}}{ff}}$ is

$${\tilde{ {\mathcal H} }}_{{{{{{\rm{e}}}}}}ff}= \;\tilde{J}\mathop{\sum} \limits_{\langle {{{{{\bf{r}}}}}},{{{{{\bf{r}}}}}}{^\prime}\rangle ,\nu }{a}_{\nu }{\tilde{{{{{{\boldsymbol{b}}}}}}}}_{{{{{{\bf{r}}}}}}}^{{{\dagger}} }{\tilde{{{{{{\mathcal{S}}}}}}}}^{\nu }{\tilde{{{{{{\boldsymbol{b}}}}}}}}_{{{{{{\bf{r}}}}}}}{\tilde{{{{{{\boldsymbol{b}}}}}}}}_{{{{{{\bf{r}}}}}}{^\prime}}^{{{\dagger}} }{\tilde{S}}^{\nu }{\tilde{{{{{{\boldsymbol{b}}}}}}}}_{{{{{{\bf{r}}}}}}{^\prime}}\\ +\tilde{J}{^\prime}\mathop{\sum} \limits_{\langle {{{{{\bf{r}}}}}},{{{{{\bf{r}}}}}}{^\prime}\rangle ,\nu }{b}_{\nu }{\widetilde{{{{{{\boldsymbol{b}}}}}}}}_{{{{{{\bf{r}}}}}}}^{{{\dagger}} }{\widetilde{{{{{{\mathcal{S}}}}}}}}^{\nu }{\widetilde{{{{{{\boldsymbol{b}}}}}}}}_{{{{{{\bf{r}}}}}}}{\widetilde{{{{{{\boldsymbol{b}}}}}}}}_{{{{{{\bf{r}}}}}}{^\prime}}^{{{\dagger}} }{\widetilde{S}}^{\nu }{\widetilde{{{{{{\boldsymbol{b}}}}}}}}_{{{{{{\bf{r}}}}}}{^\prime}}\\ +\tilde{D}\mathop{\sum} \limits_{{{{{{\bf{r}}}}}}}(1-{\tilde{{{{{{\boldsymbol{b}}}}}}}}_{{{{{{\bf{r}}}}}}}^{{{\dagger}} }\tilde{{{{{{\mathcal{A}}}}}}}{\tilde{{{{{{\boldsymbol{b}}}}}}}}_{{{{{{\bf{r}}}}}}}),$$

(8)

where ${\widetilde{{{{{{\mathcal{S}}}}}}}}^{\nu }=U{{{{{{\mathcal{S}}}}}}}^{\nu }{U}^{{{\dagger}} }$, $\widetilde{{{{{{\mathcal{A}}}}}}}=U{{{{{\mathcal{A}}}}}}{U}^{{{\dagger}} }$, and ${{{{{{\mathcal{A}}}}}}}_{\alpha \beta }={\delta }_{\alpha ,0}{\delta }_{\beta ,0}.$ Note that the unitary transformation (7) is chosen in such a way that the ${\widetilde{b}}_{{{{{{\bf{r}}}}}},0}$ boson is macroscopically occupied, namely $\langle {\widetilde{b}}_{{{{{{\bf{r}}}}}},0}\rangle =\langle {\widetilde{b}}_{{{{{{\bf{r}}}}}},0}^{{{\dagger}} }\rangle \simeq \sqrt{M}$. According to the constraint (4), $M=1$ for the case of interest. However, we will keep using $M$ because $1/M$ is the parameter of the perturbative expansion that will be introduced below. Note that $M=2S$ for the usual SU(2) spin-wave theory. The main assumption behind the $1/M$ expansion is that $\langle {\widetilde{b}}_{{{{{{\bf{r}}}}}},-1}^{{{\dagger}} }{\widetilde{b}}_{{{{{{\bf{r}}}}}},-1}\rangle ,\langle {\widetilde{b}}_{{{{{{\bf{r}}}}}},+1}^{{{\dagger}} }{\widetilde{b}}_{{{{{{\bf{r}}}}}},+1}\rangle \ll M$. Under this assumption, we can expand the spin operators ${S}^{\mu }$ and the quadrupolar operator ${({S}^{z})}^{2}$ in powers of $1/M$ (see Note 5 in the Supplementary Information). The resulting expansion of ${\widetilde{{{{{{\mathcal{H}}}}}}}}_{{{{{{\rm{e}}}}}}{ff}}$ is

$${\widetilde{{{{{{\mathcal{H}}}}}}}}_{{{{{{\rm{e}}}}}}{ff}}={M}^{2}{{{{{{\mathcal{H}}}}}}}^{(0)}+M{H}^{(2)}+{M}^{1/2}{H}^{(3)}+{M}^{0}{H}^{(4)}+O\left({M}^{-1}\right),$$

(9)

where the linear term ${H}^{(1)}$ vanishes because the parameters $\theta$ and $\phi$ in Eq. (5) are determined by minimizing the mean-field energy

$${{{{{{\mathcal{H}}}}}}}^{(0)}=(2\tilde{J}\widetilde{\Delta }+\widetilde{J}^{\prime} \widetilde{\Delta }^{\prime} ){{{\sin }}}^{4}\theta {{{\cos }}}^{2}2\phi\quad\quad\quad$$

$$\qquad \quad\quad\quad\quad\quad-\frac{1}{2}(2\widetilde{J}+\widetilde{J}^{\prime} ){{{\sin }}}^{2}2\theta (1+{{\sin }}2\phi )+\widetilde{D}{{{\sin }}}^{2}\theta .$$

(10)

Since the AFM order is invariant under time reversal followed by one lattice translation, the states ${S}^{z}={\pm}\!1$ must have equal weight in the mean-field state (5), implying that $\phi =\pi /4$. By minimizing ${{{{{{\mathcal{H}}}}}}}^{(0)}$ with respect to $\theta$, we obtain

$$x\equiv {{{\sin }}}^{2}\theta =\frac{1}{2}-\frac{\widetilde{D}}{8(2\widetilde{J}+\widetilde{J}^{\prime} )}.$$

(11)

The quadratic term ${{{{{{\mathcal{H}}}}}}}^{(2)}$ represents the generalized linear spin-wave (GLSW) Hamiltonian. After Fourier transforming the bosonic operators,

$${\tilde{b}}_{{{{{{\bf{r}}}}}}\alpha }=\frac{1}{\sqrt{{N}_{s}}}\mathop{\sum}\limits _{{{{{{\bf{k}}}}}}}{\tilde{b}}_{{{{{{\bf{k}}}}}}\alpha }{e}^{i{{{{{\bf{k}}}}}}\cdot {{{{{\bf{r}}}}}}},$$

(12)

where ${N}_{s}$ is the number of sites, ${{{{{{\mathcal{H}}}}}}}^{(2)}$ can be brought into a compact form by introducing the Nambu spinor ${\vec{b}}_{{{{{{\bf{k}}}}}}}={({\widetilde{b}}_{{{{{{\bf{k}}}}}},+1},{\widetilde{b}}_{{{{{{\bf{k}}}}}},-1},{\widetilde{b}}_{-{{{{{\bf{k}}}}}},+1}^{{{\dagger}} },{\widetilde{b}}_{-{{{{{\bf{k}}}}}},-1}^{{{\dagger}} })}^{T}$,

$${ {\mathcal H} }^{(2)}=\mathop{\sum} \limits_{{{{{{\bf{k}}}}}}}\mathop{\sum} \limits_{\alpha ,\beta =\pm 1}{\overrightarrow{b}}_{\!\!{{{{{\bf{k}}}}}}}^{{{\dagger}} }{ {\mathcal H} }^{(2)}({{{{{\bf{k}}}}}}){\overrightarrow{b}}_{\!\!{{{{{\bf{k}}}}}}},$$

(13)

with

$${ {\mathcal H} }^{(2)}({{{{{\bf{k}}}}}})=\bigg(\begin{array}{ll}{\Delta }_{\alpha \beta }({{{{{\bf{k}}}}}}) & {\Lambda }_{\alpha \beta }({{{{{\bf{k}}}}}})\\ {\Lambda }_{\beta \alpha }({{{{{\bf{k}}}}}}) & {\Delta }_{\beta \alpha }({{{{{\bf{k}}}}}})\end{array}\bigg).$$

(14)

The matrix elements are

$${\Delta }_{\alpha \beta }({{{{{\bf{k}}}}}})= \mathop{\sum} \limits_{\nu }\big[(2{a}_{\nu }\tilde{J}+{b}_{\nu }\tilde{J}^\prime)({\widetilde{{{{{{\mathcal{S}}}}}}}}_{\alpha \beta }^{\nu }{\widetilde{{{{{{\mathcal{S}}}}}}}}_{00}^{\nu }-({\widetilde{{{{{{\mathcal{S}}}}}}}}_{00}^{\nu })^{2}{\delta }_{\alpha \beta })\\ +(\tilde{J}{a}_{\nu }\mathop{\sum} \limits_{{\nu }^\prime=x,y}\,\cos \,{k}_{{\nu }^\prime}+\tilde{J}^\prime{b}_{\nu }\,\cos \,{k}_{z}){\widetilde{{{{{{\mathcal{S}}}}}}}}_{\alpha 0}^{\nu }{\widetilde{{{{{{\mathcal{S}}}}}}}}_{0\beta }^{\nu }\big]\\ -\frac{\tilde{D}}{2}({\tilde{{{{{{\mathcal{A}}}}}}}}_{\alpha \beta }-{\tilde{{{{{{\mathcal{A}}}}}}}}_{00}{\delta }_{\alpha \beta }),$$

(15)

$${\Lambda }_{\alpha \beta }({{{{{\bf{k}}}}}})=\mathop{\sum} \limits_{\nu }{\tilde{{{{{{\mathcal{S}}}}}}}}_{\alpha 0}^{\nu }{\tilde{{{{{{\mathcal{S}}}}}}}}_{\beta 0}^{\nu }\left[\tilde{J}{a}_{\nu }\mathop{\sum}\limits _{\nu ^\prime=x,y}\,\cos \,{k}_{\nu ^\prime}+\tilde{J}{^\prime}{\rm b}_{\nu }\,\cos \,{k}_{\rm z}\right].$$

(16)

The collinear mean-field state (5) has a residual Z₂ symmetry associated with a $\pi$-rotation along the direction of the ordered moments (local $\widetilde{z}$ axis). The bosonic operator ${\widetilde{b}}_{+1}^{{{\dagger}} }$ picks up minus sign under this Z₂ symmetry because it creates the state with ${\widetilde{S}}^{z}=-1$. In contrast, the bosonic operator ${\widetilde{b}}_{-1}^{{{\dagger}} }$ remains invariant because it creates the state with ${\widetilde{S}}^{z}=0$. This symmetry analysis implies that the ${\widetilde{b}}_{+1}$ and ${\widetilde{b}}_{-1}$ bosons must be decoupled in ${{{{{{\mathcal{H}}}}}}}^{(2)}$ because a non-vanishing hybridization term would otherwise break this Z₂ symmetry:

$${ {\mathcal H} }^{(2)}=\mathop{\sum} \limits_{{{{{{\bf{k}}}}}},\alpha =\pm 1}[{A}_{{{{{{\bf{k}}}}}},\alpha }{\tilde{b}}_{{{{{{\bf{k}}}}}},\alpha }^{{{\dagger}} }{\tilde{b}}_{{{{{{\bf{k}}}}}},\alpha }-\frac{{B}_{{{{{{\bf{k}}}}}},\alpha }}{2}({\tilde{b}}_{-{{{{{\bf{k}}}}}},\alpha }{\tilde{b}}_{{{{{{\bf{k}}}}}},\alpha }+{\tilde{b}}_{{{{{{\bf{k}}}}}},\alpha }^{{{\dagger}} }{\tilde{b}}_{-{{{{{\bf{k}}}}}},\alpha }^{{{\dagger}} })]$$

(17)

with ${\gamma }_{{{{{{\bf{k}}}}}}}^{{xy}}={{\cos }}({k}_{x})+{{\cos }}({k}_{y})$, ${\gamma }_{{{{{{\bf{k}}}}}}}^{z}={{\cos }}({k}_{z})$ and the expressions for ${A}_{{{{{{\bf{k}}}}}},\alpha }$ and ${B}_{{{{{{\bf{k}}}}}},\alpha }$ are given in Note 5 of the Supplementary Information.

The diagonal form of ${{{{{{\mathcal{H}}}}}}}^{(2)}$,

$${ {\mathcal H} }^{(2)}=\mathop{\sum}\limits _{{{{{{\bf{k}}}}}},\alpha =\pm 1}{\omega }_{{{{{{\bf{k}}}}}},\alpha }\bigg({\beta }_{{{{{{\bf{k}}}}}},\alpha }^{{{\dagger}} }{\beta }_{{{{{{\bf{k}}}}}},\alpha }+\frac{1}{2}\bigg)$$

(18)

is then obtained by applying an independent Bogoliubov transformation for each bosonic flavor,

$${\widetilde{b}}_{{{{{{\bf{k}}}}}},\pm 1}={u}_{{{{{{\bf{k}}}}}},\pm 1}{\beta }_{{{{{{\bf{k}}}}}},\pm 1}+{v}_{{{{{{\bf{k}}}}}},\pm 1}{\beta }_{-{{{{{\bf{k}}}}}},\pm 1}^{{{\dagger}} },$$

(19)

with

$${u}_{{{{{{\bf{k}}}}}},\pm 1}=\sqrt{\frac{1}{2}\bigg(\frac{{{{{{\rm{|}}}}}}{A}_{{{{{{\bf{k}}}}}},\pm 1}{{{{{\rm{|}}}}}}}{{\omega }_{{{{{{\bf{k}}}}}},\pm 1}}+1\bigg)},$$

$${v}_{{{{{{\bf{k}}}}}},\pm 1}=\frac{{B}_{{{{{{\bf{k}}}}}},\pm }}{{{{{{\rm{|}}}}}}{B}_{{{{{{\bf{k}}}}}},\pm }{{{{{\rm{|}}}}}}}\sqrt{\frac{1}{2}\bigg(\frac{{{{{{\rm{|}}}}}}{A}_{{{{{{\bf{k}}}}}},\pm 1}{{{{{\rm{|}}}}}}}{{\omega }_{{{{{{\bf{k}}}}}},\pm 1}}-1\bigg)}.$$

(20)

The operators ${\beta }_{{{{{{\bf{k}}}}}},\pm 1}^{{{\dagger}} }$ create quasi-particles with energy

$${\omega }_{{{{{{\bf{k}}}}}},\pm 1}=\sqrt{{A}_{{{{{{\bf{k}}}}}},\pm 1}^{2}-{B}_{{{{{{\bf{k}}}}}},\pm 1}^{2}},$$

(21)

where ${\omega }_{{{{{{\bf{k}}}}}},+1}({\omega }_{{{{{{\bf{k}}}}}},-1})$ is the dispersion relation of the transverse (longitudinal) modes. The neutron-scattering intensity $I({{{{{\bf{Q}}}}}},\omega )$ is related to the spin-spin correlation function through

$$I({{{{{\bf{Q}}}}}},\omega )\propto {f}^{2}({{{{{\bf{Q}}}}}})\mathop{\sum}\limits_{\mu ,\nu }\left({\delta }_{\mu \nu }-\frac{{\hat{Q}}_{\mu }{\hat{Q}}_{\nu }}{{Q}^{2}}\right)\\ \times \frac{1}{2\pi {N}_{s}}\mathop{\sum }\limits_{i,j}^{{N}_{s}}{\int }_{-\infty }^{+\infty }dt{e}^{i\omega t-i{{{{{\bf{Q}}}}}}\cdot ({{{{{{\bf{r}}}}}}}_{i}-{{{{{{\bf{r}}}}}}}_{j})}\langle {s}_{i}^{\mu }(t){s}_{j}^{\nu }(0)\rangle ,$$

(22)

where ${{{{{\bf{Q}}}}}}$ is the momentum vector transfer, and $f({{{{{\bf{Q}}}}}})$ is the magnetic form factor of Fe²⁺. In the Discussion section, we will show that although the GLSW approach discussed in this section can reproduce the dispersion relations of all observed low-energy modes in Ba₂FeSi₂O₇, it cannot account for various interaction effects that are revealed by the INS experiments. To capture these effects, we must then include the next order terms in the $1/M$-expansion.

Nonlinear corrections

In this section, we demonstrate that the generalization of the $1/S$-expansion is simply a loop expansion. Based on this result, we compute the one-loop corrections to the linear theory presented in the previous section. As we explain in the next section, the one-loop correction accounts for both the broadening and the energy renormalization of the longitudinal mode near the zone center.

After Fourier transforming and applying a Bogoliubov transformation, the cubic contributions to the generalized spin-wave theory become

$${{{{{{\mathcal{H}}}}}}}^{(3)}={{{{{{\mathcal{H}}}}}}}_{c}^{(3)}+{{{{{{\mathcal{H}}}}}}}_{l}^{(3)},$$

(23)

with

$${ {\mathcal H} }_{c}^{(3)}= \frac{1}{\sqrt{{N}_{s}}}\mathop{\sum} \limits_{{{{{{{\bf{q}}}}}}}_{i}}\mathop{\sum}\limits _{{\alpha }_{i}=\pm 1}\delta ({{{{{{\bf{q}}}}}}}_{1}+{{{{{{\bf{q}}}}}}}_{2}+{{{{{{\bf{q}}}}}}}_{3})\\ \times \big[\frac{1}{3!}{V}_{s}^{(3)}({{{{{{\bf{q}}}}}}}_{1,2,3},{\alpha }_{1,2,3}){\beta }_{{{{{{{\bf{q}}}}}}}_{1},{\alpha }_{1}}{\beta }_{{{{{{{\bf{q}}}}}}}_{2},{\alpha }_{2}}{\beta }_{{{{{{{\bf{q}}}}}}}_{3},{\alpha }_{3}}\\ +\frac{1}{2!}{V}_{d}^{(3)}({{{{{{\bf{q}}}}}}}_{1,2,3},{\alpha }_{1,2,3}){\beta }_{{\bar{{{{{{\bf{q}}}}}}}}_{1},{\alpha }_{1}}^{{{\dagger}} }{\beta }_{{\bar{{{{{{\bf{q}}}}}}}}_{2},{\alpha }_{2}}^{{{\dagger}} }{\beta }_{{{{{{{\bf{q}}}}}}}_{3},{\alpha }_{3}}+h.c.\big],$$

(24)

and

$${ {\mathcal H} }_{l}^{(3)} = \frac{1}{\sqrt{{N}_{s}}}\mathop{\sum} \limits_{{{{{{\bf{q}}}}}}}\mathop{\sum}\limits _{\alpha =\pm 1}[{V}_{l}^{(3)}({{{{{\bf{q}}}}}},0,{{{{{\bf{q}}}}}};\alpha ,-1,\alpha ){\beta }_{0,-1}^{{{\dagger}} }+h.c.]\\ = \sqrt{{N}_{s}}\mathop{\sum} _{\alpha =\pm 1}[{V}_{L,\alpha }{\beta }_{0,-1}^{{{\dagger}} }+h.c.].$$

(25)

Here ${V}_{d}^{(3)}$ and ${V}_{s}^{(3)}$ are the decay and sink vertices, respectively. The symmetry-allowed cubic vertices are depicted in the second and third lines of Fig. 4. Note that, unlike the SU(2) case, collinear magnetic ordering does not preclude cubic terms in the expansion (9) of the generalized SU($N$) spin-wave theory with $N \; > \; 2$. For the SU(3) case under consideration, the residual Z₂ symmetry ($\pi$-rotation along the local $\widetilde{z}$-axis) only requires that the ${\widetilde{b}}_{+1}$ boson must appear an even number of times (e.g., ${\widetilde{b}}_{+1}{\widetilde{b}}_{+1}$ or ${\widetilde{b}}_{+1}^{{{\dagger}} }{\widetilde{b}}_{+1}^{{{\dagger}} }$) in the cubic terms. ${{{{{{\mathcal{H}}}}}}}_{l}^{(3)}$ is a linear term that originates from the normal ordering of the cubic vertices. This term renormalizes the optimal value $\theta$ that was obtained from the minimization of ${{{{{{\mathcal{H}}}}}}}^{(0)}$. The integral of ${V}_{l}^{(3)}({{{{{\bf{q}}}}}}{{{{{\rm{;}}}}}}\alpha ,-1)$ over the entire Brillouin zone is the so-called cubic–linear vertex, which is nonzero only for the longitudinal boson at the ordering wave vector ${{{{{\bf{q}}}}}}$ = 0 (in the twisted frame). The explicit forms of ${V}_{d,s}^{(3)}$ and ${V}_{l}^{(3)}$ are derived in Note 7 of the Supplementary Information.

We will now describe the construction of a systematic perturbative field theory that is controlled by $1/M$. This scheme can be applied to study anharmonicities starting from any generalized spin-wave theory based on a Schwinger boson representation of the generators of SU($N$). The well-known $1/S$-expansion will be recovered for the particular case $N=2$ and $M=2S$. As we will demonstrate below, the $1/M$-expansion is just a particular example of the loop expansion that is commonly used to describe spontaneous symmetry breaking in particle theory⁴². The connection is more evident after noticing that $M$ becomes an overall prefactor of the rescaled Hamiltonian (Eq. (9)), $H={H}_{{{{{{\rm{e}}}}}}{ff}}/M$, once we also rescale the bosonic fields according to ${\bar{b}}_{{{{{{\bf{r}}}}}},\nu }={\widetilde{b}}_{{{{{{\bf{r}}}}}},\nu }/\sqrt{M}$. Since the original interaction vertices ${V}^{(n)}(n\ge 3)$ scale as ${V}^{(n)}\sim {(M)}^{2-\frac{n}{2}}$, all vertices of the rescaled Hamiltonian $H({\bar{b}}_{{{{{{\bf{r}}}}}},\nu },{\bar{b}}_{{{{{{\bf{r}}}}}},\nu }^{{{\dagger}} }{{{{{\rm{\}}}}}}})$ become of order $M$, while the propagator is still of order $1/M$. Thus, the order $p$ of a particular one-particle irreducible diagram is $V-I$, where $V$ is the number of vertices and $I$ is the number of internal lines (note that the frequency $\omega$ is of order ${M}^{0}$ because the quadratic contribution $\langle {H}^{(2)}\rangle$ is independent of $M$). Since the number of loops is $L=I-V+1$ (Every vertex introduces a delta function that reduces the number of independent momenta by one, except for one delta function that is left over for overall energy-momentum conservation), we obtain the desired result: $p=1-L$.

Let us rederive this result without rescaling the fields and the Hamiltonian. As we already mentioned, Eq. (9) tells us that the interaction vertices ${V}^{(n)}(n\ge 3)$ scale as ${V}^{(n)}\sim {(M)}^{2-\frac{n}{2}}$. The quasi-particle propagator

$${{{{{{\mathcal{G}}}}}}}_{0,\alpha }({{{{{\bf{k}}}}}},i\omega )={(-i\omega +{\omega }_{{{{{{\bf{k}}}}}},\alpha })}^{-1},\alpha ={\pm}\! 1$$

(26)

where $\omega$ is the Matsubara frequency, scales as ${{{{{{\mathcal{G}}}}}}}_{0,\alpha }(k)\sim {M}^{-1}$ because ${\omega }_{{{{{{\bf{k}}}}}},\alpha }$ is of order $M$ (see Eq. (9)). The dressed single-particle propagator is obtained from the Dyson equation,

$${{{{{{\mathcal{G}}}}}}}^{-1}({{{{{\bf{k}}}}}},i\omega )={{{{{{\mathcal{G}}}}}}}_{0}^{-1}({{{{{\bf{k}}}}}},i\omega )-\Sigma ({{{{{\bf{k}}}}}},i\omega ),$$

(27)

where $\Sigma ({{{{{\bf{k}}}}}},i\omega )$ is the single-particle self-energy. At a given order in $M$, the dressed propagator includes two external legs, $L$ independent loops, $I$ internal lines (bare propagators ${{{{{{\mathcal{G}}}}}}}_{0}$) and ${V}_{n}$ interaction vertices of the type ${V}^{(n)}$. After a summation over the Matsubara frequency $\omega \sim {M}^{1}$, each loop gives a contribution of order ${M}^{1}$. Hence, the order $p$ of a particular one-particle irreducible diagram contributing to $\Sigma ({{{{{\bf{k}}}}}},i\omega )$ is

$$p=L-I+\mathop{\sum}\limits _{n\ge 3}{V}_{n}\bigg[2-\frac{n}{2}\bigg].$$

(28)

Since each internal line connects a pair of vertices, we have

$$\mathop{\sum}\limits _{n\ge 3}\,n{V}_{n}=2I+2,$$

(29)

where ${\sum }_{n\ge 3}n{V}_{n}$ is the total number of lines. Furthermore, the number of loops is equal to the number of independent momentum integrals. From the conservation of momentum at each vertex, we have

$$L=I-\bigg[\mathop{\sum} \limits_{n\ge 3}\,{V}_{n}-1\bigg].$$

(30)

By combining the above results, we obtain

$$p=1+\mathop{\sum} \limits_{n\ge 3}{V}_{n}-\mathop{\sum}\limits _{n\ge 3}\frac{n{V}_{n}}{2}=1-L,$$

(31)

implying that the order of a given diagram is determined by the number of loops.

The lowest-order ${{{{{\mathcal{O}}}}}}$ (M⁰) Feynman diagrams are shown in Fig. 6. Since the inverse of the bare boson propagator is of order ${{{{{\mathcal{O}}}}}}$ (M¹), the remaining diagrams of order ${{{{{\mathcal{O}}}}}}$ (M⁰) give a relative 1/M-correction to the poles of the bare propagators. The real part of the new poles corresponds to the renormalized single-particle energy, whereas the imaginary part corresponds to the decay rate, which is responsible for the broadening of the quasi-particle peaks measured with INS.

Fig. 5: Basic ingredients of the perturbative field theory in 1/M for Ba₂F_eSi₂O₇.

Solid (dash) lines represent the bare propagator of the transverse (longitudinal) boson. The symmetry-allowed cubic vertices are shown on the second and third lines. The red (blue) dot represents a decay (sink) vertex. The cubic–linear vertices are listed on the fourth line. The last line represents the normal vertex ${V}_{\alpha \alpha }^{(4,N)}$ from ${{{{{{\mathcal{H}}}}}}}^{(4)}$.

The contributions to the self-energy from the decay and from the source diagrams shown in Fig. 6 are

$${\Sigma }_{\alpha }^{d}({{{{{\bf{q}}}}}},i\omega )=\frac{1}{2{N}_{s}}\mathop{\sum}\limits _{{{{{{\bf{k}}}}}},{\alpha }_{1},{\alpha }_{2}=\pm 1}\frac{|{V}_{d}^{(3)}(\bar{{{{{{\bf{k}}}}}}},{{{{{\bf{k}}}}}}+\bar{{{{{{\bf{q}}}}}}},{{{{{\bf{q}}}}}};{\alpha }_{1},{\alpha }_{2},\alpha ){|}^{2}}{i\omega -{\omega }_{{{{{{\bf{k}}}}}},{\alpha }_{1}}-{\omega }_{{{{{{\bf{q}}}}}}+\bar{{{{{{\bf{k}}}}}}},{\alpha }_{2}}},$$

(32)

and

$${\Sigma }_{\alpha }^{s}({{{{{\bf{q}}}}}},i\omega )=-\frac{1}{2{N}_{s}}\mathop{\sum}\limits _{{{{{{\bf{k}}}}}},{\alpha }_{1},{\alpha }_{2}=\pm 1}\frac{|{V}_{s}^{(3)}({{{{{\bf{k}}}}}},\bar{{{{{{\bf{k}}}}}}}+\bar{{{{{{\bf{q}}}}}}},{{{{{\bf{q}}}}}};{\alpha }_{1},{\alpha }_{2},\alpha ){|}^{2}}{i\omega +{\omega }_{{{{{{\bf{k}}}}}},{\alpha }_{1}}+{\omega }_{{{{{{\bf{q}}}}}}+\bar{{{{{{\bf{k}}}}}}},{\alpha }_{2}}},$$

(33)

Fig. 6: Diagrammatic representation of the Dyson equation.

a One-loop diagrams that contribute up to the order ${M}^{0}$ for the transverse boson. b One-loop diagrams that contribute up to the order ${M}^{0}$ for the longitudinal boson. The dressed propagator is denoted by a thick line, whereas the bare propagator is denoted by a thin line.

respectively.

Finally, the diagrams that appear in the last line for both panels of Fig. 6 arise from the normal ordering of the quartic term ${{{{{{\mathcal{H}}}}}}}^{(4)}$ in Eq. (9). These contributions simply renormalize the quadratic Hamiltonian:

$${ {\mathcal H} }_{NO}^{(4)}=\mathop{\sum}\limits _{{{{{{\bf{q}}}}}},\alpha ,\alpha {^\prime}}[{V}_{\alpha \alpha {^\prime}}^{(4,N)}{\beta }_{{{{{{\bf{q}}}}}},\alpha }^{{{\dagger}} }{\beta }_{{{{{{\bf{q}}}}}},\alpha {^\prime}}+({V}_{\alpha \alpha {^\prime}}^{(4,A)}{\beta }_{-{{{{{\bf{q}}}}}},\alpha }{\beta }_{{{{{{\bf{q}}}}}},\alpha {^\prime}}+h.c.)],$$

(34)

where ${V}_{\alpha \alpha {\prime} }^{(4,N)}({V}_{\alpha \alpha {\prime} }^{(4,A)})$ represents the normal (anomalous) contributions. Since ${{{{{{\mathcal{H}}}}}}}_{{NO}}^{(4)}$ is of order M⁰, only the diagonal normal contribution arising from the normal vertex ${V}_{\alpha \alpha {\prime} }^{(4,N)}{\delta }_{\alpha ,\alpha {\prime} }$ gives a relative correction of order $1/M$ to the bare single-particle energy given in Eq. (21) (the anomalous terms in Eq. (34) give a relative correction contribution order $1/{M}^{2}$). The derivation of ${V}_{\alpha \alpha }^{(4,N)}$ is included in Note 7 of the Supplementary Information.

We note the parallel between the decay, sink, and quartic diagrams that give the 1/M-correction to the single-particle self-energy and the ones that appear in the 1/S-expansion of the standard SU(2) spin-wave theory of non-collinear Heisenberg magnets¹¹. The main difference is that the SU(3) theory includes an extra bosonic flavor that enables more symmetry-allowed decay channels. In addition, the cubic–linear diagram exists even in absence of magnetic field because the magnitude of the ordered magnetic moment can be renormalized by changing the variational parameter θ. These diagrams, shown in the third line of Fig. 6a and the fourth line of Fig. 6b, are obtained by contracting one of the legs of the decay vertex with the cubic–linear vertex shown in Fig. 5. By using the Feynman rules, the cubic–linear diagrams are calculated as

$${\Sigma }_{\alpha }^{{cl}}({{{{{\bf{q}}}}}})=-\frac{1}{{\omega }_{{{{{{\bf{0}}}}}},-1}}({[{V}_{d}^{(3)}({{{{{\bf{0}}}}}},\bar{{{{{{\bf{q}}}}}}},{{{{{\bf{q}}}}}}{{{{{\rm{;}}}}}}\alpha ,-1,\alpha )]}^{\ast }{V}_{L,\alpha }+h.c.).$$

(35)

By applying the analytic continuation $\omega \pm i{\delta }^{+}\to i\omega$ and adopting the so-called on-shell approximation $\omega ={\omega }_{{{{{{\bf{q}}}}}}}$ for Eq. (32) and Eq. (33), the renormalized pole of the dressed propagator ${{{{{\mathcal{G}}}}}}$ is calculated as ${\widetilde{\omega }}_{{{{{{\bf{q}}}}}},\alpha }-i{\widetilde{\Gamma }}_{{{{{{\bf{q}}}}}},\alpha }={\omega }_{{{{{{\bf{q}}}}}},\alpha }+{V}_{{{{{{\bf{q}}}}}},\alpha \alpha }^{(4,N)}+{\Sigma }_{\alpha }^{{cl}}({{{{{\bf{q}}}}}})+{\Sigma }_{\alpha }^{s}({{{{{\bf{q}}}}}},{\omega }_{{{{{{\bf{q}}}}}},\alpha })+{\Sigma }_{\alpha }^{d}({{{{{\bf{q}}}}}},{\omega }_{{{{{{\bf{q}}}}}},\alpha })$, where the imaginary part of the pole ${\widetilde{\Gamma }}_{{{{{{\bf{k}}}}}},\alpha }$ arises from the decay term ${\Sigma }_{\alpha }^{d}$, that accounts for the observed broadening of the longitudinal mode in most regions of the BZ (see Fig. 3c, f) (the calculations are summarized in Note 9 of the Supplementary Information and ref. ⁵²). Moreover, the shift in the real part of the pole implies a corresponding renormalization in the model parameters. By fitting the neutron-scattering data with the renormalized dispersion peaks ${\widetilde{\omega }}_{{{{{{\bf{q}}}}}},\alpha }$ at the ZC, we obtain the set of optimal Hamiltonian parameters listed as set ${{{{{\mathcal{B}}}}}}$ in Table 1 and discussed further below.

Table 1 Parameter sets of GLSWT and GLSWT+one-loop models.

Discussion

Comparison between experiment and theory

To understand the spin excitation spectrum of Ba₂FeSi₂O₇ and demonstrate the importance of using the one-loop corrections, we start the comparison between experiment and theory with the GLSWT (i.e. without one-loop corrections). Figure 3b and e shows contour plots of $I({{{{{\bf{Q}}}}}},\omega )$ (Eq. (22)) calculated with the GLSWT along the [$H$, 0, 0.5], and [$H$, $H$, 0.5] direction, respectively. The Hamiltonian parameters (see set ${{{{{\mathcal{A}}}}}}$ in Table 1) are extracted by fitting the measured positions of the quasi-particle peaks (Gaussian-fitted peak centers of the experimental data) at the ZC. The GLSWT reproduces the dispersion of the observed two transverse modes T₁ and T₂ along the [H, 0, 0.5] and [H, H, 0.5] directions (Fig. 3b, e). Noticeably, the calculated longitudinal mode closely reproduces the experimental dispersion of the “L”-mode, which demonstrates that the SU(3) spin-wave theory describes the quasi-particles in Ba₂FeSi₂O₇.

Notably, the GLSWT does not reproduce the broadening and renormalization of the longitudinal modes observed in the inelastic neutron-scattering data. This is because the effect arises from the decay of a longitudinal mode into two transverse modes that is induced by the cubic term ${{{{{{\mathcal{H}}}}}}}^{(3)}$ of the expansion (Eq. (9)). To capture this effect, the $1/M$-correction from the one-loop expansion (see “Nonlinear correction” section) must be included. The GLSWT+one-loop correction can then describe the broadened spectrum of the longitudinal mode. The new Hamiltonian parameters, which are determined via the same procedure that is described above (see set ${{{{{\mathcal{B}}}}}}$ in Table 1), allow us to reproduce the observed spectrum (see Fig. 3c, f).

A more in-depth comparison between theory and experiment is shown in Fig. 7a and b. These figures show the quasi-particle dispersions along the [$H$, 0, 0.5] direction calculated with the GLSWT and GLSWT plus one-loop corrections compared to the measured dispersion. Near the ZC, ${{{{{{\bf{Q}}}}}}}_{{{{{{\bf{m}}}}}}}$ = (1, 0, 0.5) the energy of longitudinal mode obtained from the GLSWT is noticeably higher than the peak center of the measured “$L$”-mode (orange dots). The discrepancy in the dispersion is resolved by introducing the one-loop corrections. The real part of the self-energy renormalizes the energy of the longitudinal mode, leading to a better agreement with the observed peak positions near the ZC. At the same time, the imaginary part of the self-energy obtained from the decay diagrams, ${\Sigma }_{\alpha }^{d},$ leads to an intrinsic line broadening of the longitudinal mode that is missing in the GLSWT. In Fig. 7b and d, the lower (upper) boundary of the red-shaded region is given by ${\widetilde{\omega }}_{{{{{{\bf{k}}}}}},-1}\left(\mp \right){\widetilde{\Gamma }}_{{{{{{\bf{k}}}}}},-1},$ representing theoretical line broadening of the longitudinal mode that is compared against the experimental FWHM (orange error bars). In particular, the above-mentioned effects are most striking at ${{{{{{\bf{Q}}}}}}}_{{{{{{\bf{m}}}}}}}$ = (1, 0, 0.5), therefore we present a comparison of the intensity line-cut at this momentum transfer in Fig. 7e. It is interesting to note that the energy shift of the transverse mode is also captured by the one-loop corrections.

Fig. 7: Comparison between measured and calculated spectrum.

Comparison of the measured and calculated dispersion along the [$H$, 0, 0.5] (a, b) and [$H$, $H$, 0.5] (c, d) directions. In all panels of this figure, the theoretical results are obtained for the parameter set ${{{{{\mathcal{B}}}}}}$ in Table 1. a–d Blue and orange filled circles indicate the measured transverse and longitudinal modes, obtained from the Gaussian fitting to the data shown in Fig. 4a. Dots and error bars indicate peak centers and full width at half maxima (FWHM) of the observed modes, respectively. Lines indicate the calculated dispersions obtained from the GLSWT and GLSWT+one-loop corrections. The red-shaded region in b and d depict the decay (line broadening) of the longitudinal mode given by the one-loop corrections. e–g Comparison between the measured (blue dots) and calculated (orange and black lines) INS intensities at three high-symmetric ${{{{{\bf{Q}}}}}}$-points at (1, 0, 0.5), (0, 0, 0.5), and (0.5, 0.5, 0.5). All the experimental data were measured using CTAX with fixed ${E}_{{{{{{\rm{f}}}}}}}$ = 3 meV. For GLSWT, two transverse and longitudinal modes are denoted with ${T}_{1}$, ${T}_{2}$, and L.

After verifying that the one-loop corrections can simultaneously capture the broadening of the longitudinal mode and the energy shift of both the transverse and the longitudinal modes at the magnetic ZC, it is natural to ask if this also holds true far away from the ZC. Figure 7f, g is the intensity cuts for two representative points on the ZB. At a first glance, the peak centers of both modes are reasonably reproduced by the one-loop corrections. A more detailed analysis reveals that the experimental FWHM of both peaks is equal to the instrumental resolution. However, as illustrated in Fig. 8a, since the longitudinal modes are still inside the two-magnon continuum, the one-loop correction predicts an intrinsic broadening (black curves) in Fig. 7f, g.

Fig. 8: Kinematic constraints for the decay of the longitudinal mode.

The blue (orange) curve shows the calculated transverse (longitudinal) band dispersions along [H, 0, 0.5] with the GLSWT (using parameters set ${{{{{\mathcal{B}}}}}}$ in Table 1). The light blue-shaded areas indicate the two-transverse mode continuum, whose lower edge is indicated with a black solid line (${E}_{2}^{{{\min }}}$). a Results of the effective $S=1$ model. b Same as a but for a gapped branch of transverse modes (an ad hoc gap has been added to Eq. (21)).

To understand the origin of this discrepancy, we trace back the decay channel of the longitudinal mode on the zone boundaries. The two-magnon continuum at the zone edge starts at an energy equal to the sum of the single-magnon energies at the zone center and the zone boundary. Due to the U(1) symmetry of the effective Hamiltonian, the magnons are gapless at the zone center, implying that the onset of the two-magnon continuum coincides with the single-magnon branch (see Fig. 8). In absence of U(1) symmetry, the magnon modes become gapped and the longitudinal mode does not need to lie inside the two-magnon continuum for arbitrary values of the wave vector (see Fig. 8b). A small magnon gap pushes the onset of the two-magnon continuum to be above the energy of the longitudinal mode at the zone boundaries. This modification of the two-magnon spectrum precludes the decay of the longitudinal mode near the zone boundary and explains the experimental observation. We then conjecture that the single-magnon dispersion is indeed gapped.

Unfortunately, it is difficult to extract the size of this gap from our INS data because of the large quasi-elastic scattering. Nevertheless, the analysis presented in Note 2 of the Supplementary Information indicates that our data are indeed consistent with a finite spin gap. We note that the gap can be captured by working with the original spin $S=2$ Hamiltonian (Eq. (1)). The tetragonal symmetry allows for a single-ion anisotropy term of the form ${{{{{{\mathcal{H}}}}}}}_{A}=A{\sum }_{i}[{({S}_{i}^{x})}^{4}+{({S}_{i}^{y})}^{4}]$, which breaks the global U(1) symmetry, generating a finite gap for the transverse mode. However, when we project the original $S=2$ Hamiltonian onto the low-energy space to obtain the effective spin $S=1$ Hamiltonian (Eq. (2)), the term ${{{{{{\mathcal{H}}}}}}}_{A}$ simply renormalizes the single-ion anisotropy, implying that the low-energy model acquires an “emergent” U(1) symmetry that is absent in the original high-energy model. Lastly, we note that the energies of the longitudinal mode on the zone boundaries after the one-loop corrections are slightly lower than the measured values. This level of discrepancy can be attributed to the missing second-order corrections ${{{{{\rm{O}}}}}}\left(\frac{{J}^{2}}{3D}\right)$ to the low-energy model (2) or to missing terms in the original Hamiltonian (1). A simple analysis shows that a second nearest-neighbor AFM interaction with ${\widetilde{J}}_{2}\sim 0.2\widetilde{J}$ can account for this discrepancy. For simplicity, ${\widetilde{J}}_{2}$ is not included in our calculation. Except for the discrepancy near the zone boundaries, the effective $S=1$ model with one-loop corrections successfully captures most features of the INS data inside the BZ.

Finally, we emphasize that the loop expansion preserves the Goldstone mode that results from the spontaneous breaking of the emergent U(1) symmetry group of ${\widetilde{{{{{{\mathcal{H}}}}}}}}_{{{{{{\rm{e}}}}}}{ff}}$. More specifically, the ${{{{{\mathscr{O}}}}}}\left({M}^{0}\right)$ correction to the real part of the self-energy vanishes for the Goldstone mode (see Note 8 in the Supplementary Information), although the individual contributions from the diagrams shown in Fig. 5 diverge as 1/q in the long-wavelength limit. We note that previous attempts of computing the decay of the longitudinal mode³⁴ have not accounted for the renormalization of the single-particle dispersion arising from the 1/M-correction to the real part of the self-energy. This correction leads to a significant change in the extracted ratio $\alpha =\widetilde{J}/\widetilde{D}$ of Ba₂FeSi₂O₇, cf. ${\alpha }_{{{{{{\rm{GLSWT}}}}}}}=0.152,$ and ${\alpha }_{{{{{{\rm{GLSWT}}}}}}+{{{{{\rm{one}}}}}}-{{{{{\rm{loop}}}}}}}=0.187.$ This change is a direct consequence of the substantial renormalization of the energy ${\omega }_{L}\left({{{{{{\bf{Q}}}}}}}_{{{{{{\bf{m}}}}}}}\right)$ of the longitudinal mode at the ZC. In fact, an accurate calculation that goes beyond the one-loop approximation estimates that the critical ${{{{{{\rm{\alpha }}}}}}}_{{{{{{\rm{c}}}}}}}$ required to close the gap ${\omega }_{L}({{{{{{\bf{Q}}}}}}}_{{{{{{\bf{m}}}}}}})$ for $\widetilde{J}^{\prime} =0.1\widetilde{J}$, and $\widetilde{\Delta }=\widetilde{\Delta }^{\prime} =1/3$ is around 0.158. In other words, the Hamiltonian parameters extracted from fitting the experiment with the GLSWT place Ba₂FeSi₂O₇ on the quantum paramagnetic side of the phase diagram shown in Fig. 1, which obviously contradicts the experimental evidence. In contrast, the set of parameters obtained from the GLSWT+one-loop correction (${\alpha }_{{{{{{\rm{GLSWT}}}}}}+{{{{{\rm{one}}}}}}-{{{{{\rm{loop}}}}}}}$) place the material at the magnetically ordered phase of the exact phase diagram. Furthermore, the calculated ordered moment is very close to the measured value 2.95 ${\mu }_{{{{{{\rm{B}}}}}}}$ (see Note 12 of the Supplementary Information for discussion of the reduction of the ordered moment). In general, nonlinear corrections become increasingly important upon approaching the QCP and logarithmic corrections due to multi-loop vertex renormalizations become relevant very close to this point^28,31,53,54. The fact that one-loop correction is enough to reproduce the spectrum of Ba₂FeSi₂O₇ indicates that this material is still far enough from that critical regime.

In summary, Ba₂FeSi₂O₇ provides a natural realization of a quasi-2D easy-plane antiferromagnet in the proximity of the QCP that signals the transition into the QPM phase. Previous examples of low-dimensional easy-plane quantum magnets in the proximity of this QCP were typically located on the quantum paramagnetic side of the quantum phase transition^17,20,21,23. Ba₂FeSi₂O₇ then allows us to explain the strong decay and renormalization effects of the low-energy transverse and longitudinal modes of the AFM state. Furthermore, the distance to the O(2) QCP could be in principle controlled by chemical substitution, while the application of an in-plane magnetic field, that gaps out the transverse modes, can be used to control the decay rate of the longitudinal mode.

Here, we have used the INS data of Ba₂FeSi₂O₇ as a platform to test a loop expansion based on an SU(3) spin-wave theory^17,18,20,55, that captures the longitudinal and the transverse modes at the linear level. This loop expansion, which generalizes the well-known $1/S$-expansion of the SU(2) spin-wave theory, allows us to reproduce the measured width and renormalization of the longitudinal and transverse modes near the zone center by just including a one-loop correction. Small discrepancies near the zone boundary are attributed to limitations of the effective low-energy $S=1$ model that we adopted for this work.

The loop expansion that we have described in this manuscript provides a general scheme for treating quantum magnets with more than one type of low-energy mode. In general, quantum magnets that exhibit low-energy modes with $N-1$ different “flavors” can be treated semi-classically using an SU(N) spin-wave theory. The parameter of the semi-classical expansion is the number of loops in the Feynman diagrams that contribute to the single-particle propagator.

Methods

Sample preparation

A single crystal of Ba₂FeSi₂O₇ was grown using an optical floating zone melting method⁴⁴. Polycrystalline Ba₂FeSi₂O₇ feed-rods were prepared using the solid-state reaction method. The stoichiometric powders of BaCO₃ and Fe₂O₃, and SiO₂ were mixed, ground, pelletized, and sintered with intermediate heating in a reduced gas atmosphere. Ba₂FeSi₂O₇ single crystal was grown using a floating zone furnace in the same gas environment.

Inelastic neutron-scattering measurement

Inelastic neutron-scattering measurements were performed 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) at Oak Ridge National Laboratory⁴⁵. A 2.15-g single crystal was aligned with the ($H,H,L$) and ($H,0,L$) in the horizontal scattering plane for CTAX and HYSPEC experiments. A liquid helium cryostat was used to control temperature. At CTAX, the initial neutron energy was selected using a PG (002) monochromator, and the final neutron energy was fixed to ${E}_{{{{{{\rm{f}}}}}}}$ = 3.0 meV by a PG (002) analyzer. The horizontal collimation was guide-open-40′−120′, which provides an energy resolution with a full width half maximum (FWHM) = 0.1 and 0.18 meV for $\Delta E$ = 0 and 2.5 meV, respectively. For the HYSPEC experiment, ${E}_{{{{{{\rm{i}}}}}}}$ = 9 meV and a Fermi chopper frequency of 300 Hz were used, which provides an energy resolution of FWHM = 0.28 meV and 0.19 meV at $\Delta E$ = 0 and 2.5 meV, respectively. Measurements were performed at $T$ = 1.6 K and 90 K by rotating the sample from −50 to 170° with 1° steps. Data were symmetrized over positive and negative $H$ and integrated over $K$ = [−0.1, 0.1] and $L$ = [0.4, 0.6]. In Fig. 3a, there appears to be quasi-elastic scattering below 0.5 meV in low ${{{{{\bf{Q}}}}}}$-region. This scattering arises from the incompletely blocked direct beam due to the oscillating collimator. All of the datasets were reduced and analyzed using MANTID⁵⁶ and DAVE⁵⁷.