Anharmonic Collective Oscillations in Isotropic Spin Systems and their Spectroscopic Signatures

TL;DR

This paper reveals that isotropic spin systems can exhibit soft quartic spin modes leading to a temperature-dependent spin-wave gap, expanding the understanding of spin dynamics and suggesting new spectroscopic signatures.

Contribution

It introduces a theoretical framework for soft spin modes in isotropic systems with no fine-tuning, supported by numerical simulations and phenomenological analysis.

Findings

Quartic potentials induce a temperature-dependent spin-wave gap.

Numerical models show gap growth with temperature in spiral spin states.

Predictions can be tested via inelastic neutron scattering experiments.

Abstract

Spin waves are the fundamental excitations in magnetically ordered spin systems and are ubiquitously observed in magnetic materials. However, the standard understanding of spin waves as collective spin oscillations in an effective harmonic potential does not consider the possibility of soft modes, such as those due to an effective quartic potential. In this work, we show that such quartic potentials arise under very general conditions in a broad class of isotropic spin systems without a fine-tuning of the interaction parameters. Considering models with spin spiral ground states in two and three spatial dimensions, we numerically demonstrate that quartic amplitude spin oscillations produce a fluctuation-induced spin-wave gap which grows with temperature according to a characteristic power-law. In conjunction with a phenomenological theory, the present work provides a general theoretical framework for describing soft spin modes, extending the previously discussed spin dynamics in the presence of order-by-disorder, and highlighting the important role of finite-size effects. Our predictions of a temperature-dependent gap in spiral spin systems could be tested in inelastic neutron scattering experiments, providing direct spectroscopic evidence for thermal effects arising from soft spin modes in magnetic materials.

Equations

$$H=J_1\sum _{{⟨i<j⟩}_1}{\mathbf{S}}_i\cdot {\mathbf{S}}_j+J_2\sum _{{⟨i<j⟩}_2}{\mathbf{S}}_i\cdot {\mathbf{S}}_j+J_3\sum _{{⟨i<j⟩}_3}{\mathbf{S}}_i\cdot {\mathbf{S}}_j.$$

$${\mathbf{S}}_i^{\text{GS}}=\left(\begin{array}{c}\cos (\mathbf{Q}\cdot {\mathbf{r}}_i+\varphi )\\ \sin (\mathbf{Q}\cdot {\mathbf{r}}_i+\varphi )\\ 0\end{array}\right),$$

$${\mathbf{S}}_i^{\text{quartic}}=\frac{1}{N_i}\left[{\mathbf{S}}_i^{\text{GS}}+\delta \left(\begin{array}{c}0\\ 0\\ \cos ({\mathbf{Q}}_{\mathrm{aB}}\cdot {\mathbf{r}}_i+{\varphi }^{\prime })\end{array}\right)\right].$$

$$H_4=\frac{p^2}{2m}+\lambda x^4,$$

$$H_{\mathrm{eff}}=H_4+\alpha T{x}^2=\frac{p^2}{2m}+\lambda x^4+\alpha T{x}^2.$$

$$\Delta =\sqrt{\frac{\alpha T}{m}}{\left(1+\sqrt{1+\frac{2T^{\ast }}{T}}\right)}^{1/2},$$

$${\mathbf{S}}_i$$

$${\hat{\mathbf{x}}}_i$$

$${\hat{\mathbf{y}}}_i$$

$${\hat{\mathbf{z}}}_i$$

$$H_2=\frac{1}{2}\sum _{i,\nu }{J}_{\nu }\left[\cos ({\mathbf{Q}}_{\mathrm{B}}\cdot \mathbf{\nu })\left(x_i{x}_{i+\nu }-x_i^2-p_i^2\right)+p_i{p}_{i+\nu }\right],$$

$$x_i=\frac{1}{\sqrt{N}}\sum _{\mathbf{k}}{e}^{i\mathbf{k}\cdot {\mathbf{r}}_i}{x}_{\mathbf{k}},p_i=\frac{1}{\sqrt{N}}\sum _{\mathbf{k}}{e}^{i\mathbf{k}\cdot {\mathbf{r}}_i}{p}_{\mathbf{k}},$$

$$H_2=\frac{1}{2}\sum _{\mathbf{k}}\left[L_{\mathbf{k}}{x}_{-\mathbf{k}}{x}_{\mathbf{k}}+M_{\mathbf{k}}{p}_{-\mathbf{k}}{p}_{\mathbf{k}}\right],$$

$$L_{\mathbf{k}}$$

$$M_{\mathbf{k}}$$

$$x_i=\frac{1}{\sqrt{N}}\left(e^{i{\mathbf{Q}}_{\mathrm{aB}}\cdot {\mathbf{r}}_i}{x}_{{\mathbf{Q}}_{\mathrm{aB}}}+e^{-i{\mathbf{Q}}_{\mathrm{aB}}\cdot {\mathbf{r}}_i}{x}_{-{\mathbf{Q}}_{\mathrm{aB}}}\right),$$

$$p_i=\frac{1}{\sqrt{N}}\left(e^{i{\mathbf{Q}}_{\mathrm{aB}}\cdot {\mathbf{r}}_i}{p}_{{\mathbf{Q}}_{\mathrm{aB}}}+e^{-i{\mathbf{Q}}_{\mathrm{aB}}\cdot {\mathbf{r}}_i}{p}_{-{\mathbf{Q}}_{\mathrm{aB}}}\right).$$

$$x_{{\mathbf{Q}}_{\mathrm{aB}}}$$

$$p_{{\mathbf{Q}}_{\mathrm{aB}}}$$

$$x_i=\sqrt{\frac{2}{N}}\left(X\cos ({\mathbf{Q}}_{\mathrm{aB}}\cdot {\mathbf{r}}_i)-\tilde{P}\sin ({\mathbf{Q}}_{\mathrm{aB}}\cdot {\mathbf{r}}_i)\right),$$

$$p_i=\sqrt{\frac{2}{N}}\left(P\cos ({\mathbf{Q}}_{\mathrm{aB}}\cdot {\mathbf{r}}_i)+\tilde{X}\sin ({\mathbf{Q}}_{\mathrm{aB}}\cdot {\mathbf{r}}_i)\right).$$

$${\mathbf{S}}_i^{\text{quartic}}=\frac{1}{\sqrt{1+{\delta }^2{\cos }^2\left(\frac{2\pi }{L}x\right)}}\left(\begin{array}{c}\cos \left(\frac{2\pi }{L}y\right)\\ \sin \left(\frac{2\pi }{L}y\right)\\ \delta \cos \left(\frac{2\pi }{L}x\right)\end{array}\right),$$

$$J_n\sum _i\sum _{j\in n^{\text{th}}\text{ n.n of }i}{\mathbf{S}}_i\cdot {\mathbf{S}}_j.$$

$$\sum _{x,y}\mathbf{S}(x,y)\cdot \mathbf{S}(x+a,y+b)+\mathbf{S}(x,y)\cdot \mathbf{S}(x-b,y+a).$$

$$\sum _{x,y}$$

$$-\frac{{\delta }^2}{2}\cos \left(\frac{2\pi }{L}b\right)\left[{\cos }^2\left(\frac{2\pi }{L}(x+a)\right)+{\cos }^2\left(\frac{2\pi }{L}x\right)\right]$$

$$-\frac{{\delta }^2}{2}\cos \left(\frac{2\pi }{L}a\right)\left[{\cos }^2\left(\frac{2\pi }{L}(x-b)\right)+{\cos }^2\left(\frac{2\pi }{L}x\right)\right].$$

$$\sum _{x,y}\frac{{\delta }^2}{2}\cos \left(\frac{4\pi }{L}x\right)\left[\cos \left(\frac{2\pi }{L}b\right){\sin }^2\left(\frac{2\pi }{L}a\right)+\sin \left(\frac{2\pi }{L}b\right){\cos }^2\left(\frac{2\pi }{L}a\right)\right]=0.$$

$$e(\delta )=\sum _{i=1}^3{e}_i(\delta )\approx -1.00729+0.03906{\delta }^4+O({\delta }^6),$$

$${\mathbf{S}}_i^{\text{quartic}}=\frac{1}{\sqrt{1+{\delta }^2({\cos }^2(\frac{2\pi }{L}x)+{\cos }^2(\frac{2\pi }{L}z))}}\left(\begin{array}{c}\cos (\frac{2\pi }{L}y)\\ \sin (\frac{2\pi }{L}y)\\ \delta (\cos (\frac{2\pi }{L}x)+\cos (\frac{2\pi }{L}z))\end{array}\right),$$

$$e(\delta )=e_1(\delta )+e_2(\delta )+e_4(\delta )\approx -1.48114+0.153863{\delta }^4+O({\delta }^6).$$

$$\frac{d{\mathbf{S}}_i}{dt}={\mathbf{S}}_i\times {\mathbf{h}}_i,$$

$$\mathcal{S}(\mathbf{q},\omega )=\frac{1}{2\pi N}\sum _{i,j=1}^N{\int }_{-\infty }^{\infty }{e}^{i\omega t}{e}^{-i\mathbf{q}({\mathbf{r}}_i-{\mathbf{r}}_j)}⟨\mathbf{S}{}_i(0)\cdot \mathbf{S}{}_j(t)⟩dt$$

Full text

Introduction.— The notion of a quasiparticle is widely used to describe collective excitations in many-body systems [1]. This term refers to a common situation in which, similarly to fundamental particles in high-energy physics, a collective excitation in a solid-state system can be described by characteristic properties such as its effective mass or, more generally, its dispersion relation. Well-known examples of quasiparticles include phonons, which arise from collective lattice vibrations, and magnons, which result from collective spin precessions [2]. In both cases, the quasiparticles can be classically understood as oscillations in an effective harmonic potential $V(x)=\alpha x^{2}$ with a frequency $\omega$, and with the quantum mechanical excitation energy, $E$, given by $E=\hbar\omega$.

In this Letter, we address the fundamental question of how the concept of a quasiparticle changes when the effective potential is softer than quadratic and how such anharmonicity can be exposed in experiments. For concreteness, we specifically consider the case of quartic potentials where $V(x)\propto x^{4}$. Such potentials and the associated collective excitations are of relevance in quite distinct physical systems, including trapped Bose-Einstein condensates [3, 4], chaotic systems [5, 6], soft phonons [7, 8], and, in high-energy physics, in the context of the Higgs self-coupling [9, 10, 11]. Here, we investigate quartic potentials in classical Heisenberg spin systems, which emerge under fairly general conditions and without the need to fine-tune model parameters to eliminate the harmonic $\alpha x^{2}$ component. Our starting point is a Heisenberg spin Hamiltonian that has several degenerate symmetry-related planar spin spiral ground states. We show that in a system exhibiting a planar spin-spiral long-range ordered state at one of several symmetry-related ground-state ordering wave vectors, an oscillation of the spin component perpendicular to the spiral ordering plane, taken at the wave vector of a different ground-state spiral, has a quartic dependence of the excitation energy on the amplitude. We refer to this oscillation as a quartic oscillation in the following.

Since a system with a purely quartic potential has a vanishing harmonic component ($\alpha=0$), one might naively expect it to display gapless excitations. While this is true in the zero-temperature limit, our numerical simulations reveal that thermal fluctuations dynamically generate a spin-wave gap $\Delta$. This implies a remarkable situation in which the elementary properties of a quasiparticle, such as its effective mass, are no longer determined by the system’s microscopic parameters, but instead become state-dependent, allowing them to be modified by external parameters such as the temperature $T$. Our numerical simulations of finite size spin systems in two and three spatial dimensions indicate a characteristic scaling of this fluctuation-induced spin wave gap of the form $\Delta(T)\sim T^{\mu}$ where $\mu=1/4$ ($\mu\approx 1/2$) at low (intermediate) temperatures. We explain these results using a phenomenological theory that takes into account the quartic oscillations and their coupling to other quadratic modes that inevitably exist in the system. Specifically, this theory reveals an intriguing finite-size effect resulting from the small number [$O(1)$] of quartic modes that couple to $O(N)$ quadratic modes (where $N$ is the total number of spins). Consequently, although the low-temperature gap scaling with $\mu=1/4$ is a direct consequence of the quartic oscillation, the low-temperature regime in which it occurs vanishes in the thermodynamic limit, leaving only the thermally generated gap with $\mu=1/2$.

Through our detailed understanding of the gap scaling associated with quartic spin oscillations, we extend previous theoretical work on a related type of fluctuation-induced spin-wave gap that has been predicted in the context of a thermal or quantum order-by-disorder effect due to an exact but accidental ground-state zero mode, giving rise to a so-called pseudo-Goldstone mode [12, 13, 14, 15]. In the thermal classical case, the spin-wave gap was found to scale with the same characteristic exponent $\mu=1/2$ [14]. This implies that pseudo-Goldstone modes do not require the existence of exact ground-state zero modes, which are rarely found in real materials. Instead, they can be observed under less stringent conditions in the presence of quartic potentials, which we demonstrate can occur under rather general assumptions for isotropic spin systems.

On the experimental front, the dynamical gap generation atop an anharmonic potential we study may have already been detected in neutron scattering experiments on magnetic materials such as the pyrochlore helimagnet ZnCr2Se4 [16, 17] and the half-Heusler compound GdPtBi [18]. Here, we provide the theoretical basis for these observations by relating them to quartic oscillations using a classical description to predict the ensuing temperature dependence of the spin-wave gap.

Spin Systems and Quartic Perturbations.— To set the stage for our study of anharmonic spin oscillations in isotropic spin systems, we first consider a minimal two-dimensional (2D) square lattice model that supports soft modes in a quartic potential. Later, we will extend this analysis to a three-dimensional (3D) cubic lattice to confirm that our results apply to a wider range of spin systems and the more common 3D situation. In the square lattice case, we consider a Heisenberg Hamiltonian with exchange couplings extending up to third-nearest neighbors, with ferromagnetic $J_{1}<0$ and antiferromagnetic $J_{2}$, $J_{3}>0$ couplings:

[TABLE]

Here, $\langle\ldots\rangle_{n}$ indicates the sum over the $n^{\textrm{th}}$-nearest neighbors, and the spins are treated classically as three-dimensional vectors with unitary norm $\lvert{\bm{S}}_{i}\rvert=1$. The ground-state spin configuration for Eq. (1) is a coplanar spiral [19],

[TABLE]

where $\bm{Q}=(Q_{x},Q_{y})$ is the spiral wave vector, $\bm{r}_{i}$ is the position of site $i$, and $\phi\in[0,2\pi)$ is a phase. Without loss of generality, we assume the spiral to be in the $x$-$y$ plane and also set $\phi=0$. The wave vector $\bm{Q}$ is then determined by substituting the spiral spin configuration into the Hamiltonian and identifying the $\bm{Q}$ that minimizes the system’s energy. We fix $J_{1}=-1$, $J_{2}=0.476393205$, and $J_{3}=\frac{1}{2}J_{2}-\frac{1}{5}$, such that the ground state is a one-dimensional spiral with two possible commensurate wave vectors $\bm{Q}=(0,\pm\frac{2\pi}{5})$ and $(\pm\frac{2\pi}{5},0)$, related by a $\pi/2$ lattice rotation (see Fig. 1). While the commensurability of the spiral is irrelevant to the physical properties that we study herein, it tremendously simplifies the numerical calculations as it enables periodic boundary conditions.

The continuous O(3) symmetry of the Hamiltonian cannot be broken at finite temperatures in two dimensions, excluding a transition into a long-range magnetically ordered phase 111The absence of a phase transition in 2D follows from the Mermin-Wagner-Hohenberg theorem [28, 29]. Below a certain temperature, finite systems appear magnetically ordered when the correlation length becomes comparable to the system size. However, true long-range order is absent in the thermodynamic limit.. However, the system undergoes a finite-temperature Ising transition associated with the discrete $\mathbb{Z}_{2}$ symmetry between the two (energetically degenerate) spiral wave vector directions [21]. Specifically, the one-dimensional spiral phase breaks the square lattice’s fourfold rotational symmetry down to a twofold symmetry, with the spiral oriented along either of the Cartesian directions (see Fig. 1). The selected ordering wave vector, $\bm{Q}_{\rm B}$, is responsible for the Bragg peaks of the ground-state order, and we refer to the wave vector that is not selected as $\bm{Q}_{\rm aB}$, where “aB” stands for anti-Bragg.

This system has two key properties that are essential in harboring the physics of anharmonic collective excitations that we aim to explore. First, the system supports spin oscillations in an effective quartic potential which correspond to spin fluctuations in the direction perpendicular to the planar spiral ground state in Eq. (2) and characterized by the wave vector $\bm{Q}_{\rm aB}$. Specifically, we consider $\bm{S}^{\text{quartic}}_{i}$ as a perturbed state about the long-range ordered ground state,

[TABLE]

Here $\delta$ is the perturbation amplitude, $N_{i}$ is the site ($i$)-dependent normalization to ensure $\lvert{\bm{S}}^{\text{quartic}}_{i}\rvert=1$, and $\phi^{\prime}\in[0,2\pi)$. We show in the Supplemental Material (SM) [22] that such a perturbed non-coplanar state always has a vanishing $\delta^{2}$ term in the energy, regardless of the range of the interactions, as long as the Hamiltonian conserves the $C_{4}$, reflection, and translation symmetries of the square lattice. For our particular set of interactions, the energy per site is $e=e_{\rm GS}+0.039\delta^{4}+O(\delta^{6})$ [22]. Second, aside from three zero modes arising from perturbations with the wave vector $\bm{Q}_{\rm B}$ — which correspond to global spin rotations that preserve the spiral order and would represent real Goldstone modes at $\bm{Q}_{\rm B}$ were the system magnetically long-range ordered at $T>0$ — the system does not exhibit any additional modes softer than the quartic ones discussed above. In particular, the system has no accidental ground-state degeneracies 222Accidental continuous ground-state degeneracies are known to occur in systems with fine-tuned interactions that display simpler ground-state wave vectors such as $\bm{Q}=(0,\pm\pi)$ and $\bm{Q}=(\pm\pi,0)$ [33]. By considering systems with spin spiral ground states, we avoid such accidental zero modes. leading to zero modes, such as occurs in the context of order-by-disorder phenomena [14]. Therefore, the system considered allows us to focus specifically on the effects of quartic modes, which have heretofore attracted limited attention.

To further establish the general scope and quantitative validity of our findings, we extend our analysis to a 3D classical Heisenberg model on a cubic lattice. As in the 2D case, we consider a model with general spiral ground states. This is achieved by using a Hamiltonian with the following interaction parameters: $J_{1}=-1$, $J_{2}=4/[4(4+1.8\cos{\gamma})]$, and $J_{4}=1.8/[4(4+1.8\cos{\gamma})]$, where $J_{n}$ is the coupling between $n$th-nearest neighbors [16]. In this case, the ground state is also given by a coplanar one-dimensional spiral given by Eq. (2), but with ordering wave vectors $\bm{Q}=(0,0,\pm\gamma)$, $(0,\pm\gamma,0)$ or $(\pm\gamma,0,0)$. We fix $\gamma=2\pi/5$, yielding again a five-site periodic spiral. Below a critical temperature $T_{c}$, one of these three spiral configurations is selected, and the perpendicular spin fluctuations with the not-selected wave vectors correspond to quartic oscillations [22].

There are two main differences between this 3D model and the 2D model discussed above. First, quartic oscillations can now involve contributions from two wave vectors, corresponding to the not-selected spirals in the ordered state. Second, in 3D, the SO(3) spin symmetry is explicitly broken in the ordered phase below $T_{c}$, giving rise to true magnetic long-range order (we find that both types of symmetry breaking, lattice rotation and spin rotation, occur at the same critical temperature $T_{c}$).

Numerical simulations.— To investigate the dynamic signatures of quartic modes, we perform classical Monte Carlo (cMC) [24] and spin molecular dynamics (sMD) [25] simulations on the models introduced above. Using equilibrated states generated from cMC as initial configurations for the sMD simulations, we compute the spin structure factor, $\mathcal{S}(\bm{Q},\omega)$, as a function of temperature [22]. Specifically, our focus is on the determination of the excitation gap $\Delta$ at the wave vector $\bm{q}=\bm{Q}_{\rm aB}$ of the quartic perturbation in Eq. (3), which can be extracted by fitting the spectra with a Lorentzian function [22].

In the 2D case, the cMC results show that at $T_{c}=0.071(3)$ the system breaks the fourfold lattice rotation symmetry and develops a one-dimensional spiral along the $\hat{\bm{x}}$ or $\hat{\bm{y}}$ direction. Below $T_{c}$, $\mathcal{S}(\bm{Q}_{\rm aB},\omega)$ displays a peak at a finite frequency that decreases with temperature [see Fig. 2(a)]. Extracting the value of the gap for each temperature results in Fig. 2(b), where we show the gap $\Delta$ as a function of the temperature in a log-log plot for several square systems with linear sizes $L=20$, 40, 60, and 80. This gap can be attributed to an anharmonic effect because linear spin wave theory would predict a gapless mode at this point as it treats quartic oscillations as zero modes [16].

In the 3D case, we use cubic systems with $N=L^{3}$ spins to compute $\mathcal{S}(\bm{Q}_{\rm aB},\omega)$. The cMC calculations show that at $T_{c}\simeq 0.17(1)$, the system orders in a one-dimensional spiral along $\hat{\bm{x}}$, $\hat{\bm{y}}$ or $\hat{\bm{z}}$. As expected, the dynamical structure factor at $\bm{q}=\bm{Q}_{\rm aB_{1}}$ and $\bm{q}=\bm{Q}_{\rm aB_{2}}$ displays a peak at the same frequency [22]. The results are shown in Fig. 2(c) for $L=10$ and 15. The similarities between the results for the 2D and 3D models show that the gap behavior is a characteristic property of a system that possesses quartic perturbations.

Overall, our results can be summarized by three key observations, which we will explain in detail below

At low temperatures, the gap is well described by the scaling $\Delta\sim T^{1/4}$ for both the 2D and 3D models, as can be seen in Fig. 2(b) and 2(c), respectively. 2. 2.

At higher temperatures, the gap increases faster than $\Delta\sim T^{1/4}$. Remarkably, and maybe counter-intuitively, the largest gap is found right below the phase transition, although the order nearly vanishes in this regime. This indicates that the gap at $\bm{Q}_{\rm aB}$ is a robust feature of the ordered phase. 3. 3.

Finite-size effects are more prominent at low $T$, where $\Delta\sim T^{1/4}$, as seen in Fig. 2(b) and 2(c).

The $T^{1/4}$ behavior.— To explain the origin of the exponent $1/4$ observed in Fig. 2 at low temperatures, we model the system’s dynamics by a simple classical oscillation in a quartic potential and in thermal equilibrium, described by the Hamiltonian

[TABLE]

where the first term is an effective kinetic term. Then, by a straightforward dimensional analysis, one finds that the frequency of the oscillation scales as $\omega\sim\lambda^{1/4}E^{1/4}$ for an excitation of energy $E$. Associating the frequency with the gap, $\Delta\sim\omega$, and taking into account $E\sim T$, we obtain $\Delta\sim T^{1/4}$, explaining our numerical observations (see End Matter for an alternative numerical approach to the $T^{1/4}$ gap scaling). The effective Hamiltonian $H_{4}$ accounts only for the energetic contribution of the quartic mode, and is therefore expected to be accurate only at very low temperatures.

Entropic effects.— At moderate finite temperatures, the entropic effects must be considered. These contributions can be effectively modeled by including an additional $\sim Tx^{2}$ term in Eq. (4), which arises from integrating out other spin-wave modes that interact with the quartic mode. This leads to a minimal modified effective Hamiltonian

[TABLE]

The additional term is negligible at very low temperatures, but it becomes comparable to, and ultimately dominates over the quartic term as temperature increases. To determine the temperature dependence of the frequency of the oscillator in Eq. (5), we perform a mean-field decoupling of the $x^{4}$ term [22], yielding

[TABLE]

where $T^{*}\equiv 6\lambda/\alpha^{2}$. This result has two interesting limits: i) for $T\gg T^{*}$, where the entropic contribution dominates over the quartic term, $\Delta\sim\sqrt{T}$, and ii) for $T\ll T^{*}$, where the quartic term remains dominant, $\Delta\sim T^{1/4}$. Thus, the temperature $T^{*}$ sets a rough crossover energy scale below which the quartic mode gap scales as $T^{1/4}$ and above which it follows $\sqrt{T}$. This result not only explains our observation that the gap increases more rapidly at higher temperatures, but also provides a phenomenological equation that can be fitted to our numerical simulations.

In Fig. 2(b) and 2(c), we fit our numerical results for the 2D and 3D models to the prediction of the phenomenological theory of Eq. (6) [22]. For both models, we find a good agreement between the simulations and the phenomenological theory, especially at low temperatures. The region where the phenomenological model predicts a $\sim\sqrt{T}$ behavior is harder to identify, but is still identifiable in our numerical data. Specifically, the red line in Fig. 2(b) corresponds to a plain $\sim\sqrt{T}$ function demonstrating that in a small intermediate temperature regime, our results indeed follow this behavior approximately. The vertical lines in Fig. 2(b) and 2(c) show the characteristic crossover temperature $T^{*}$ obtained from our fits.

Finally, for both the 2D and 3D models, significant deviations between the simulations and the phenomenological model occur near the critical temperature $T_{c}$. This is expected, since the order parameter fluctuations become strong and the fundamental underlying assumption of a state selection with wave vector $\bm{Q}_{\rm B}$ ceases to be fulfilled.

Finite-size effects.— As stated previously, finite-size effects become important at low temperatures, which is also evident from the $L$-dependent crossover temperatures $T^{*}$ highlighted in Fig. 2(b) and 2(c). To explain this, we formally derive here the phenomenological model of Eq. (5) presented above based on heuristic arguments. The first step consists of writing the spins in terms of canonically conjugate variables $\{x_{i},p_{j}\}=\delta_{ij}$,

[TABLE]

where $\{\boldsymbol{\hat{x}}_{i},\boldsymbol{\hat{y}}_{i},\boldsymbol{\hat{z}}_{i}\}$ is any right-handed local orthogonal frame. This is a faithful representation of spins since $\{x_{i},p_{j}\}=\delta_{ij}$ ensures the canonical Poisson bracket relation for the components of $\bm{S}_{i}$. Inserting Eq. (7) into the spin Hamiltonians above, with the local frame aligned with the ground-state spiral ordering at wave vector $\bm{Q}_{\rm B}$, and Taylor expanding in $x_{i}$ and $p_{i}$, allows us to isolate the quartic mode and derive the phenomenological model (see details in End Matter). The key finding is that $\lambda$ in Eq. (5) depends on the system size as $\lambda\sim 1/N$. This result can be traced back to a prefactor $1/\sqrt{N}$ in the Fourier transform, and implies that the crossover temperature behaves as $T^{*}\sim 1/N$, which is confirmed by our numerical results [inset of Fig. 2(b)]. This also affects the low-temperature regime, where the gap behaves as $\Delta\sim(T/N)^{1/4}$.

This result not only explains our observations in the numerical calculations, but also provides important non-trivial information about the approach to the thermodynamic limit $N\to\infty$. On one hand, the gap closes in the low-temperature regime, $\Delta\to 0$ for $N\to\infty$. At the same time, the crossover temperature $T^{*}$ vanishes for $N\to\infty$. This implies that only the $\Delta\sim\sqrt{T}$ regime arising from entropic effects survives for large enough systems. Interestingly, this is the same result that was previously found for order-by-thermal-disorder cases with exact zero modes arising from continuous accidental degeneracies [14].

Although the quartic term in the energy is suppressed by a factor of $1/N$ for individual magnon excitations, it can re-emerge in other settings. For example, in the presence of true or accidental degeneracies, one can form wide domain walls connecting states within the degenerate manifold at vanishing (or nearly vanishing) cost [26]. In our case, however, domain walls between the two symmetry-related incommensurate orders cost $O(J)$ because of the quartic potential. Relatedly, interactions among parametrically pumped magnons in ferromagnets can overcome the same $1/N$ suppression when the mode population is large, leading to a range of nonlinear dynamical phenomena [27].

Conclusion.— In this Letter, we demonstrated that quartic potentials which generically arise in isotropic spiral spin systems lead to a distinctive spin-wave gap of the form $\Delta\sim\sqrt{T}$ in the thermodynamic limit. Interestingly, this gap scaling does not result directly from spin oscillations in the quartic potential, but rather from entropic effects that consider the interaction between the quartic mode and the other (macroscopic number of) harmonic modes. The absence of direct signatures of quartic oscillation (which lead to a gap scaling $\Delta\sim T^{1/4}$) in the thermodynamic limit is a subtle finite-size effect, as we showed using analytical and numerical methods.

The firm understanding of the role of quartic oscillations developed here for classical spin systems paves the way for future investigations of similar quantum effects at zero and nonzero temperature [15]. By being controlled by the inverse spin length $1/S$, such effects are expected to conspire with thermal fluctuations to again generate dynamical spin wave gaps $\Delta$. The ultimate goal of such studies will be to investigate the impact of finite spin lengths $1/S>0$ on the temperature dependence $\Delta(T)$ in isotropic spiral spin systems and, in a parallel pursuit, seek to expose such behavior in real candidate materials such as ZnCr2Se4 [17].

Acknowledgments.— We thank D. Peets and D. Inosov for helpful discussions. A. F. and M. G. G. acknowledge the use of the JUWELS cluster at the Forschungszentrum Jülich and the HPC Service of ZEDAT, Freie Universität Berlin. S. K. acknowledges financial support from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy through the Würzburg-Dresden Cluster of Excellence on Complexity and Topology in Quantum Matter – ct.qmat (EXC 2147, project-id 390858490). The work at the U. of Waterloo and the U. of Windsor was supported by the NSERC of Canada (M. J. P. G. and J. G. R.). M. J. P. G. acknowledges support from the Canada Research Chair program (Tier I).

Appendix S1 Energy expansion of the quartic perturbation

In this section, we present the series expansion of the energy for a spiral state perturbed by a perpendicular perturbation with wave vector $\bm{Q}_{\rm aB}$ in both the square and cubic lattices, demonstrating the quartic dependence of the energy of such a state on the perturbation amplitude.

S1.1 Square lattice

To illustrate how quartic modes emerge in a spiral ground state perturbed by a perpendicular perturbation with a wave vector corresponding to the non-selected spiral, $\bm{Q}_{\rm aB}$, we consider a generic spiral with period $L$. First, it is useful to rewrite the perturbed state given in Eq. (3) of the main text as:

[TABLE]

where the position of site $i$ is given by $\bm{r}_{i}=(x,y)$ and we set $\bm{Q}_{\rm B}=(0,2\pi/L)$, $\bm{Q}_{\rm aB}=(2\pi/L,0)$. Furthermore, we have set the arbitrary phases to be $\phi=\phi^{\prime}=0$. Next, consider the energy term associated with the $n$th nearest-neighbor interactions:

[TABLE]

This term is proportional to the scalar product between the spin at site $\bm{r}_{i}=(x,y)$ and its $n^{\textrm{th}}$ nearest neighbors, given by $\bm{r}_{j}=(x\pm a,y\pm b)$ and $(x\pm b,y\pm a)$. To avoid double counting, we include only one term from each pair related by a global minus sign: $(x+a,y+b)$, $(x-a,y+b)$, $(x-b,y+a)$, and $(x-b,y-a)$. These four terms form two pairs related by a $90^{\circ}$ rotational symmetry around the $\bm{\hat{z}}$ axis. To show that the quadratic term vanishes, it is sufficient to consider one such pair, summed over the entire lattice:

[TABLE]

Expanding the scalar products in a series of $\delta$, we obtain for the second-order terms the following expression:

[TABLE]

Using algebraic manipulations and the fact that summing over the entire lattice causes terms proportional to $\cos{\left(\frac{2\pi}{L}x\right)}$ and $\sin{\left(\frac{2\pi}{L}x\right)}$, or equivalent expressions for $y$, to vanish, we arrive at:

[TABLE]

This final expression vanishes when summing over the entire lattice. Therefore, while the existence of a spiral ground state depends on the Hamiltonian, the quartic nature of the mode involving a perpendicular perturbation with the wave vector of the non-selected spiral is a general property of the spiral ground state, independent of the specific Hamiltonian. The same argument also applies to two-dimensional spirals.

The energy expansion up to the fourth order for the state in Eq. (S1), specific to the square lattice discussed in the main text, is given by

[TABLE]

where $e_{i}$ with $i=1$, 2, 3 are the energy contributions from the different lattice bonds with couplings $J_{1}$, $J_{2}$, $J_{3}$, respectively.

S1.2 Cubic lattice

To compute the energy expansion, we consider a state similar to the one in Eq. (S1), generalized to the equal superposition of the two quartic perturbations $S_{i}^{z}\sim\delta[\cos(\bm{Q}_{\rm aB_{1}}\cdot\bm{r}_{i})+\cos(\bm{Q}_{\rm aB_{2}}\cdot\bm{r}_{i})]$. We rewrite the state as:

[TABLE]

where $\bm{r}_{i}=(x,y,z)$, $\bm{Q}_{\rm B}=(0,\frac{2\pi}{L},0)$, $\bm{Q}_{\rm aB_{1}}=(\frac{2\pi}{L},0,0)$, and $\bm{Q}_{\rm aB_{2}}=(0,0,\frac{2\pi}{L})$. To show that the terms proportional to $\delta^{2}$ in the expression in Eq. (S2) vanish, the procedure is analogous to the one presented for the square lattice. The only difference in the cubic lattice case is that we must consider three contributions, as the perturbation consists of two terms. Specifically, the scalar product is taken between the spin at site $\bm{r}_{i}=(x,y,z)$ and three of its $n^{\rm th}$ nearest neighbors: (1) the spin at $(x+a,y+b,z+c)$ (2) the spin at $(x-b,y+a,z+c)$ obtained by a $90^{\circ}$ rotation around the $\hat{\bm{z}}$-axis (3) the spin at $(x+a,y-c,z+b)$ obtained by a $90^{\circ}$ rotation around the $\hat{\bm{x}}$-axis.

The energy expansion up to the fourth order for the state in Eq. (S7), specific to the cubic lattice discussed in the main text, is given by

[TABLE]

Appendix S2 Details on the Monte Carlo simulations

To generate the spin configurations needed for the spin Molecular Dynamics (sMD) simulations, we perform classical Monte Carlo (cMC) calculations on the square and cubic lattices. We use a particular version of the single-spin update algorithm called the adaptive Gaussian step algorithm, which allows for a controlled and fixed $50\%$ acceptance ratio at all temperatures [1]. Additionally, for each spin-update trial, we perform two over-relaxation steps. Then, at a given temperature, we perform $5\times 10^{5}$ cMC steps consisting each of $N$ single-spin update trials and $2N$ over-relaxation steps, where $N$ is the number of spins in the system. Data for measurements is collected during the second half of the cMC steps.

For both lattices, we perform two types of simulations: cooling down and heating, both with logarithmic temperature steps. For the cooling-down simulations, we do 120 temperature steps from $T/J=2$ down to $T/J=0.01$, where $J$ is the largest coupling. Then we perform an extra 31 steps down to $T/J=10^{-12}$. For the heating simulations, we depart from a ground-state configuration and perform 34 logarithmic steps from $T/J=10^{-12}$ up to $10^{-1}$ (16 steps in the cubic case). We save the last configuration at each temperature as a benchmark point.

From each configuration at a desired temperature, we generate 120 new spin configurations for the Molecular Dynamics simulations. This is done by performing $5\times 10^{5}$ cMC steps between spin configurations, implying 60 million cMC steps between the first and last spin configurations at a given temperature.

In the case of the square lattice, we do this for $N=L^{2}$ lattices with $L=20$, 40, 60, and 80. In the case of the cubic lattice, we do this for $N=L^{3}$ lattices with $L=5$, 10, 15, and 20. In both cases, we use periodical boundary conditions, which do not frustrate the 5-site spiral ground state.

For completeness, we show in Fig. S1 the specific heat $c_{v}(T)$ from the cMC calculations; in panel (a) for the square lattice model and in panel (b) for the cubic lattice. In both cases, the lines indicate the calculations for the cooling down simulations, while the dots are the results for the heating up simulations, noting an excellent agreement between both. Also, both calculations find evidence for only a single phase transition into a spiral phase.

Appendix S3 Details on the spin Molecular Dynamics simulations

In this section, we explain how the spin Molecular Dynamics (sMD) simulations are performed and how the dynamical structure factor is computed within this framework.

S3.1 Overview of the method and numerical details

The sMD simulations consist of numerically integrating the classical equations of motion for the spins:

[TABLE]

where $\bm{h}_{i}=-\frac{\partial H}{\partial\bm{S}_{i}}$ is the local effective field acting on $\bm{S}_{i}$, arising from interactions with neighboring spins. These equations correspond to the Landau-Lifshitz equations without the damping term [2]. The numerical integration is carried out using the Runge-Kutta Cash-Karp method, implemented through the Boost C++ Library [3].

S3.2 Computing the dynamical structure factor

The central quantity of our study is the dynamical structure factor

[TABLE]

evaluated at $\bm{q}=\bm{Q}_{\rm aB}$. To compute $\mathcal{S}(\bm{Q}_{\rm aB},\omega)$, we follow the approach of Ref. [4]. Specifically, the spin Fourier transform for each spin component is computed during the simulation. At the end of the time evolution, $\mathcal{S}(\bm{Q}_{\rm aB},\omega)$ is obtained by multiplying the Fourier transform by its complex conjugate and summing over the components. For the calculation of $\mathcal{S}(\bm{Q}_{\rm aB},\omega)$ at finite temperature, the results are obtained by averaging over 120 simulations. Each simulation starts from an independent, equilibrated MC configuration, which is obtained by heating up the system starting from zero temperature.

Appendix S4 Fitting procedures

In this section, we describe the fitting procedures used in this work. First, we explain how the gap value is extracted from the dynamical structure factor. Then, we discuss how these extracted values are fitted using the function in Eq. (6) of the main text, which is derived from the phenomenological model. The results of the second fit are shown in Figs. 2(b) and 2(c).

S4.1 Gap extraction from the dynamical structure factor

The value of the gap at $\bm{q}=\bm{Q}_{\rm aB}$ is obtained by fitting $\mathcal{S}(\bm{Q}_{\rm aB},\omega)$ with a Lorentzian function multiplied by a scaling factor $A$:

[TABLE]

The value of the gap corresponds to the position of the maximum, given by $x_{0}$. For each linear size and temperature, we compute the difference between the gap value extracted from the fit and the position of the maximum of $\mathcal{S}(\bm{Q}_{\rm aB},\omega)$. Since all observed differences lie within $3\gamma/4$ (square lattice) and $\gamma/2$ (cubic lattice), we adopt these values as estimates of the error.

S4.2 Fitting the simulated gaps

The gap values are fitted using the function in Eq. (6) of the main text, which describes the gap behavior derived from the phenomenological model. At low temperatures, this function predicts $\Delta\sim T^{1/4}$. This scaling is confirmed by the linear fit of $\log{T}$ versus $\log{\Delta}$, shown by the blue line in Fig. 2(b) for the square lattice and the blue line in Fig. 2(c) for the cubic lattice. Thus, we verify that the function in Eq. (6) provides a good fit to the data at low temperatures. Starting from the lowest temperatures, we gradually increase the number of data points used in the fit. For each step, we compute $\chi^{2}/\text{df}$, where $\text{df}=n-2$ is the number of degrees of freedom (with $n$ being the number of points used, minus two since the fit depends on two parameters). We then compare $\chi^{2}/\text{df}$ with the chi-squared critical value at a 95% confidence level, denoted as $\chi_{0.05}^{2}/\text{df}$. The computed $\chi^{2}/\text{df}$ values and corresponding critical values for the two lattices are shown in Fig. S2. The fits presented in Figs. 2(b) and 2(c) correspond to the maximum number of points where $\chi^{2}/\text{df}$ remains below the critical threshold across all linear system sizes. This results in $n=19$ for the square lattice and $n=8$ for the cubic lattice.

Appendix S5 Additional data for the cubic lattice

In this section, we present additional data for the cubic lattice. Specifically, the energy gap as a function of the perturbation strength, and the gap as a function of temperature at the two anti-Bragg wave vectors.

S5.1 Dynamical gap as a function of the perturbation strength

Similarly to the square lattice discussed in the End Matter, we compute $\mathcal{S}(\bm{Q}_{\rm aB},\omega)$, starting from a configuration in which the quartic modes are excited directly, using the expression in Eq. (S7). The sMD simulations show that $\mathcal{S}(\bm{Q}_{\rm aB_{1}},\omega)$ and $\mathcal{S}(\bm{Q}_{\rm aB_{2}},\omega)$ display a peak at the same frequency, that scales linearly with $\delta$, see Fig. S3.

S5.2 Equivalence of the gap at the anti-Bragg points

As explained in the main text, the cubic lattice exhibits quartic modes at the two not-selected spiral wave vectors denoted as $\bm{Q}_{\rm aB_{1}}$ and $\bm{Q}_{\rm aB_{2}}$. As shown in Fig. S4, the gap at $\bm{Q}_{\rm aB_{1}}$ and $\bm{Q}_{\rm aB_{2}}$ is the same, which agrees with the equivalence of these points. This agreement further confirms that the initial configurations are equilibrated correctly. In Fig. 2(c) of the main text, the gap values are obtained from the average of the gaps at the two anti-Bragg points.

Appendix S6 Frequency of the quartic oscillator in the presence of the entropic term

We consider a phenomenological model of a quartic oscillator with an additional entropic contribution, described by the effective Hamiltonian

[TABLE]

where the first two terms are energetic contributions and the last term arises from entropy (see the main text). To estimate the frequency of this oscillator, we apply a mean-field decoupling of the quartic term as $x^{4}\approx 6\,\langle x^{2}\rangle\,x^{2}$, where $\langle\cdots\rangle$ denotes thermal average, assuming thermal equilibrium at temperature $T$. The factor of $6$ comes from the number of ways to contract $x^{4}$ into a product involving $\langle x^{2}\rangle\,x^{2}$. This decoupling leads to an effective quadratic Hamiltonian

[TABLE]

For such a quadratic Hamiltonian, the thermal average

[TABLE]

This equation can be solved self-consistently for $\langle x^{2}\rangle$ to find the frequency $\omega$ or, alternatively, one may rewrite Eq. (S14) as

[TABLE]

This is a quadratic equation in $\omega$, giving

[TABLE]

where $T^{*}=6\lambda/\alpha^{2}$. This frequency is identified as the thermal gap $\Delta$ of the quartic mode at temperature $T$.