In the ground state the spins are ordered: all parallel in a ferromagnet, antiparallel on alternating sublattices in an antiferromagnet (Néel). Depending on the anisotropy, the ordering direction is the \(\hat z\) axis or the \(xy\) plane. Tilt one spin slightly and the exchange coupling carries the deviation along the lattice as a wave whose precession phase advances from site to site. The simulation solves the full (nonlinear) Landau–Lifshitz equation and compares the result with linear spin-wave theory.
\(J>0\); upper sign ferromagnet (FM), lower sign antiferromagnet (AFM) · \(\lambda\): XXZ exchange anisotropy · \(D\): Dzyaloshinskii–Moriya vector along \(\hat z\) · \(h\): field along \(\hat z\) · \(K>0\) easy axis, \(K<0\) easy plane · \(|\mathbf S_n|=S\). The simulation uses \(\hbar=a=S=1\). Boundaries are periodic by default (a ring, a torus in 2D/3D); an absorbing boundary can be selected in the controls.
Every site has a frame attached to its ground state: \(\hat e_3\) is the ground-state spin direction, \(\hat e_2=\hat y\) (rotating with the spiral in a spiral state), \(\hat e_1=\hat e_2\times\hat e_3\). From the components of the deviation in this frame,
\[ \varphi_n=\arg\big(\delta\mathbf S_n\cdot\hat e_1+i\,\delta\mathbf S_n\cdot\hat e_2\big) \]is the hue. The wheel is drawn as seen from the tip of \(\hat e_3\): \(\varphi\) increases counter-clockwise from \(\hat e_1\) towards \(\hat e_2\).
For FM \(\parallel\hat z\), \(\hat e_1=\hat x\) and \(\hat e_2=\hat y\): the colour is the azimuthal angle of the spin in the laboratory.
Saturation and brightness grow with \(w=\sqrt{|\delta\mathbf S_\perp|/\max|\delta\mathbf S_\perp|}\). Undeviated spins are grey, the most deviated ones fully coloured; the square root makes small amplitudes visible. In the phase map the colour blends from the background to the full hue with \(w\): background colour means no deviation.
The scale is set every frame by the current largest amplitude, so saturation shows relative, not absolute, amplitude. Use Figure 2 for absolute values.
A mode with \(\omega>0\) rotates clockwise about \(\hat e_3\), so \(\varphi\) decreases: the colour of a site cycles red → magenta → blue → cyan → green → yellow → red. In a plane wave \(\varphi_n=\mathbf k\cdot\mathbf r_n-\omega t\), so the colours repeat once per wavelength and the drift of the stripes shows the phase velocity.
The single-spin tilt and a click tilt the spin towards \(+\hat e_1\); the excited spin is red at that instant. The whole pattern often rotates with a common factor \(e^{-i\Omega t}\); what carries meaning is the colour difference between sites.
Where the amplitude goes to zero (e.g. Bessel nodes, packet tails) the phase is undefined and the colour can jump; such sites are drawn faint anyway. With elliptical precession (easy plane, spin flop, spiral) \(\varphi\) does not advance uniformly, and some colours persist longer.
In the AFM the frame of sublattice B is reversed (\(\hat e_3=-\hat z\), \(\hat e_1=-\hat x\)); the same colour on the two sublattices means opposite \(x\) directions in the laboratory. In the figures, orange dots are the full LLG solution, the teal line is linear theory, red/orange shading marks unstable and \(\omega<0\) regions, and the orange edge band is the absorbing layer.
The bar right below the picture controls the simulation: play/pause, reset and speed. The "Numerics and display" group in the right panel (size, boundary condition, arrow thinning, display magnification) opens when its title is clicked.
Keyboard: Space play/pause · R reset · +/− speed · V 3D arrows ↔ phase map (2D, 3D) · F full page · Esc leave. While a control (list or button) has the focus, the corresponding keys go to that control; the ⌨ button in the bar shows the list.
For \(J>0\) each bond energy is minimal when the spins are parallel. The ground state is
\[ \mathbf S_n = S\hat z,\qquad E_0=-N\,(JS^2+hS+KS^2). \]At \(h=K=0\) every direction is equivalent by O(3) symmetry; choosing \(\hat z\) is spontaneous symmetry breaking, and the gapless Goldstone mode at \(k\to0\) follows from it. A field \(h\) or anisotropy \(K\) breaks the symmetry explicitly and opens the gap \(\Delta=h+2KS\).
With the classical Poisson bracket \(\{S^a,S^b\}=\epsilon_{abc}S^c\) each spin precesses about its local effective field:
\[ \dot{\mathbf S}_n=\mathbf S_n\times\mathbf H_n,\quad \mathbf H_n=-\frac{\partial H}{\partial\mathbf S_n}=J(\mathbf S_{n-1}+\mathbf S_{n+1})+(h+2KS^z_n)\hat z. \]For small deviations, \(S^z_n\simeq S\), and with \(\psi_n=S^x_n+iS^y_n\)
\[ i\dot\psi_n=(2JS+\Delta)\,\psi_n-JS\,(\psi_{n+1}+\psi_{n-1}). \]This is exactly the tight-binding Schrödinger equation on a lattice: the spin deviation behaves like a particle hopping to neighbouring sites with amplitude \(JS\).
Substituting a plane wave \(\psi_n=A\,e^{i(kn-\omega t)}\) gives
\[ \omega_k=\Delta+4JS\sin^2\!\tfrac{k}{2}\;\approx\;\Delta+JS\,k^2,\qquad v_g=2JS\sin k. \]A single tilted spin \(\psi_n(0)=\psi_0\,\delta_{n0}\) contains every \(k\) with equal weight, and its exact solution is
\[ \psi_n(t)=\psi_0\, i^{\,n}\, J_n(2JSt)\; e^{-i(2JS+\Delta)t}. \]\(J_n(x)\) decays exponentially for \(|n|>x\), so the front moves along \(|n|\approx 2JSt\), at the maximum group velocity, reached at \(k=\pi/2\).
For \(H=+J\sum_n\mathbf S_n\cdot\mathbf S_{n+1}\) each bond energy is minimal when the spins are antiparallel:
\[ \mathbf S_n=(-1)^n S\hat z,\qquad E_0=-N\,(JS^2+KS^2). \]This is the classical ground state; in quantum mechanics the Néel state is not even an eigenstate of the Hamiltonian. At \(K=0\) any small \(h\) already favours the spin-flop state; the Néel state is the ground state only for \(h
The effective field is \(\mathbf H_n=-J(\mathbf S_{n-1}+\mathbf S_{n+1})+(h+2KS^z_n)\hat z\). With \(S^z_n\simeq\sigma_n S\), \(\sigma_n=(-1)^n\):
\[ \begin{aligned} i\dot\psi_A&=(2JS+2KS+h)\,\psi_A+JS\,(\psi_{n-1}+\psi_{n+1})\\ i\dot\psi_B&=-(2JS+2KS-h)\,\psi_B-JS\,(\psi_{n-1}+\psi_{n+1}) \end{aligned} \]Because \(S^z<0\) on sublattice B, its precession sense is reversed; the two sublattices couple like two particles of opposite "mass".
With \(\psi_A=A\,e^{i(kn-\omega t)}\), \(\psi_B=B\,e^{i(kn-\omega t)}\) the determinant must vanish:
\[ \varepsilon_k=2S\sqrt{(J+K)^2-J^2\cos^2k},\qquad \omega_{\pm}=\varepsilon_k\pm h. \]For \(K=0\), \(\varepsilon_k=2JS|\sin k|\): the dispersion is linear as \(k\to0\), with spin-wave velocity \(c=2JS\). Anisotropy opens the gap \(\Delta_0=2S\sqrt{K(2J+K)}\), which scales as \(\sqrt{JK}\) and is much larger than the FM gap \(2KS\) (exchange enhancement).
The magnetic unit cell is doubled, so \(\varepsilon_k\) has period \(\pi\). \(h\sum S^z\) commutes with the total \(S^z\); the field therefore only adds a uniform rotation at frequency \(h\) and splits the two branches of opposite chirality by \(\pm h\). Here \(\omega\) is the magnon energy (positive-norm branch); the \(\omega_-\) branch reaches zero at \(k=0\) when \(h=h_{sf}\), signalling the spin-flop transition.
Because these terms can change the ground state, for each parameter set the simulation first finds the ground state by energy minimisation (one angle for the FM, one per sublattice for the AFM) and then diagonalises the \(2p\times2p\) dynamical matrix in local frames for each \(k\). The closed-form results below are special cases of these numerical results.
For a uniform FM the energy difference between \(\hat z\) and the \(xy\) plane is set by a single effective anisotropy:
\[ A=J(1-\lambda)-K. \]\(A<0\): easy axis, order along \(\hat z\), and
\[ \omega_k=h+2S(\lambda J+K-J\cos k). \]\(A>0\): easy plane. At \(h=0\) the spins lie in the \(xy\) plane; for \(h>0\) they cant towards \(\hat z\) with \(\cos\theta=h/2SA\). For the AFM the same condition selects the Néel axis.
With order along \(\hat x\), the deviations are in-plane \(\delta S^y\) and out-of-plane \(\delta S^z\). The second-order energy is \(\sum_k[A_k|\delta S^y_k|^2+B_k|\delta S^z_k|^2]\), and for \(h=0\)
\[ \omega_k=2S\sqrt{J(1-\cos k)\,\big[J(1-\lambda\cos k)-K\big]}\;\approx\;\sqrt2\,S\sqrt{JA}\;|k|. \]The quadratic dispersion of the easy-axis FM becomes linear: the broken symmetry is now a continuous U(1) rotation and the order parameter is no longer conserved. The ellipse ratio is \(|\delta S^z|/|\delta S^y|=\sqrt{A_k/B_k}\to0\) as \(k\to0\); at long wavelength the spins oscillate almost entirely within the plane.
The term \(D\,\hat z\cdot(\mathbf S_n\times\mathbf S_{n+1})\) breaks inversion symmetry. With the order along \(\hat z\) (\(\mathbf M\parallel\mathbf D\)):
\[ \text{FM: }\ \omega_k=\omega_k^{0}+2DS\sin k,\qquad \text{AFM: }\ \omega_\pm(k)=\varepsilon_{\pm k}\pm h,\quad \varepsilon_k=2S\sqrt{(\lambda J+K)^2-(J\cos k+D\sin k)^2}. \]In the FM \(\omega_k\ne\omega_{-k}\): waves travelling in opposite directions have different velocities. If \(\omega_k<0\) for some \(k\), the uniform state is not the ground state; the true ground state is a spiral with \(\tan q=D/J\). In the AFM, \(\omega\) becomes complex and the Néel state dynamically unstable when \(\sqrt{J^2+D^2}>\lambda J+K\). With the order in the \(xy\) plane (\(\mathbf M\perp\mathbf D\)), DMI does not enter the linear spectrum at all; the non-reciprocity is proportional to \(\mathbf M\cdot\mathbf D\).
Each site has \(z=4\) neighbours. For order along \(\hat z\), with \(\gamma_{\mathbf k}=\tfrac12(\cos k_x+\cos k_y)\),
\[ \omega_{\mathbf k}=h+2KS+4JS(\lambda-\gamma_{\mathbf k})+2DS(\sin k_x+\sin k_y). \]At long wavelength \(\omega\approx\Delta+JS\,k^2\) is isotropic; the fourfold symmetry of the lattice appears only at short wavelength. The group velocity \(\mathbf v=2JS(\sin k_x,\sin k_y)\) is largest along the diagonal, \(2\sqrt2\,JS\).
Because the linear equation separates in \(x\) and \(y\), the exact solution for a single tilted spin is a product of Bessel functions:
\[ \psi_{xy}(t)=\psi_0\,i^{\,x+y}J_x(2JSt)\,J_y(2JSt)\,e^{-i(4\lambda JS+\Delta)t}. \]The front is not circular but the square \(|x|,|y|\le2JSt\), with edges parallel to the axes and corners along the diagonals. The square is clearly visible in the phase map. The continuum (long-wavelength) approximation cannot capture this anisotropy.
With the two sublattices \((-1)^{x+y}\):
\[ \varepsilon_{\mathbf k}=4JS\sqrt{\big(\lambda+\tfrac{K}{2J}\big)^2-\gamma_{\mathbf k}^2},\qquad \omega_\pm=\varepsilon_{\mathbf k}\pm h. \]For \(K=0,\lambda=1\) there is a linear Goldstone mode both at \(\Gamma\) and at \(M=(\pi,\pi)\), with velocity \(c=2\sqrt2\,JS\). The magnetic Brillouin zone is the square \(|k_x|+|k_y|\le\pi\), which is why the \(X\)–\(M\) and \(M\)–\(\Gamma\) segments are symmetric.
\(\mathbf S_n=R_z(\mathbf Q\cdot\mathbf r_n)\,\mathbf S^0_{s(n)}\), \(\mathbf S^0_s=S(\sin\theta_s,0,\cos\theta_s)\). For the FM, per bond,
\[ E(q,\theta)=S^2\sin^2\theta\,(-J\cos q+D\sin q)-\lambda JS^2\cos^2\theta, \]and the in-plane part reaches \(-\tilde JS^2\) at \(|q|=\arctan(D/J)\), with \(\tilde J=\sqrt{J^2+D^2}\). The spins lie in the \(xy\) plane and close into a conical spiral for \(h>0\).
Periodic boundaries require \(q=2\pi m/L\). The simulation scans every \(m\); the panel shows the locked \(q\) together with the infinite-system value.
Since \(B\) commutes with rotations about \(z\), in the rotating frame the bond becomes \(\tilde B=B\,R_z(q)\) and translation symmetry is restored. For the FM the effective in-plane coupling is \(\tilde J\) while the \(z\) coupling stays \(\lambda J\):
\[ \omega_k=2S\sqrt{\tilde J(1-\cos k)\,\big[\tilde J-\lambda J\cos k-K\big]}. \]The spiral is equivalent to an easy-plane FM with effective anisotropy \(\lambda J/\tilde J<1\). The gapless linear mode is a shift of the spiral phase (the phason). Here \(k\) is measured in the rotating frame; in the laboratory frame (e.g. neutron scattering) the branches shift to \(k\pm Q\).
In the AFM the effective in-plane coupling becomes \(\tilde J\). The Néel state becomes unstable for \(\tilde J>\lambda J+K\); this is the complex-\(\omega\) region shown in red in the dispersion plot when the spiral is not allowed. The in-plane AFM spiral is then the stable ground state.
Here the DMI vector points along \(\hat z\) on every bond, and the spiral rotates in the \(xy\) plane along \(\mathbf Q\parallel(1,1[,1])\). This is a toy model. Real interfacial DMI has \(\mathbf D\perp\) bond and in-plane, producing cycloids and Néel skyrmions; B20 chiral magnets have \(\mathbf D\parallel\) bond, giving helices and Bloch skyrmions. The ansatz covers only single-\(Q\) conical spirals; multi-\(Q\) structures and skyrmion lattices are outside it.
With \(z=6\) and \(\gamma_{\mathbf k}=\tfrac13(\cos k_x+\cos k_y+\cos k_z)\):
\[ \text{FM: }\ \omega_{\mathbf k}=\Delta+6JS(\lambda-\gamma_{\mathbf k}),\qquad \text{AFM: }\ \varepsilon_{\mathbf k}=6JS\sqrt{\big(\lambda+\tfrac{K}{3J}\big)^2-\gamma_{\mathbf k}^2}. \]The AFM ground state is the G-type (NaCl-type) Néel state \((-1)^{x+y+z}\). The dispersion is plotted along Γ–X–M–Γ–R–X with \(R=(\pi,\pi,\pi)\). The simulation displays one \(z\) = constant slice of the lattice.
A white-noise thermal field is added to the effective field; the fluctuation–dissipation theorem ties the noise amplitude to the damping:
\[ \mathbf H_n\to\mathbf H_n+\boldsymbol\xi_n,\qquad \langle\xi^a_n(t)\,\xi^b_m(t')\rangle=2\alpha k_BT\,\delta_{ab}\delta_{nm}\delta(t-t'). \]The stationary state is the canonical Boltzmann distribution. At harmonic order every magnon mode is populated by classical equipartition:
\[ \langle|\psi_{\mathbf k}|^2\rangle=\frac{2k_BT}{\omega_{\mathbf k}},\qquad \Delta m=\frac{k_BT}{N}\sum_{\mathbf k}\frac{1}{\omega_{\mathbf k}}. \]Classical spin waves follow the Rayleigh–Jeans distribution, not the Bose distribution. Hence, for the isotropic 3D FM,
\[ \Delta m_{\rm cl}=\frac{W}{6}\,\frac{k_BT}{JS^2}\simeq0.253\,\frac{k_BT}{JS^2}\;\propto T, \]where \(W\simeq1.516\) is the Watson integral. Quantum spin waves give Bloch's \(T^{3/2}\) law instead; the difference comes from modes with \(\hbar\omega\gg k_BT\) being frozen out quantum mechanically. In 1D and 2D, in the gapless case, \(\sum 1/\omega_{\mathbf k}\) diverges as \(N\to\infty\): the spin-wave counterpart of the Mermin–Wagner theorem. On a finite lattice the \(k=0\) mode is excluded.
The helix \(\varphi_n=kn\) with fixed cone angle \(\theta\) is an exact solution of the full LL equation; nonlinearity only renormalises its frequency:
\[ \omega(\theta)=h+2KS\cos\theta+2JS\cos\theta\,(1-\cos k). \]Increase θ₀ in plane-wave mode: the measured ω departs from the linear value and settles on this expression. In packets, nonlinearity competes with dispersion (in the continuum limit the LL equation is integrable and supports solitons; the discrete Heisenberg chain is not integrable).
In the ferromagnet the one-magnon state \(|k\rangle=N^{-1/2}\sum_n e^{ikn}S^-_n|0\rangle\) is an exact eigenstate of the Hamiltonian with energy exactly \(\hbar\omega_k\). The classical dispersion therefore coincides with the quantum one. This no longer holds for the antiferromagnet: the Néel state is not an eigenstate, zero-point corrections appear, and in the \(S=1/2\) chain the excitations become a spinon continuum.
By the Mermin–Wagner theorem there is no long-range order in 1D at \(T>0\); in the classical chain the correlation length at low temperature is \(\xi\simeq 2JS^2/k_BT\) (Fisher, 1964). At \(T=0\) the simulation shows the deterministic motion about the ordered state, at \(T>0\) the thermal equilibrium reached by stochastic LLG; on a finite lattice \(|m|\) stays nonzero even in 1D, because \(\xi\) can exceed the system size.
RK4 integration at \(T=0\), Heun (Stratonovich) at \(T>0\); \(|\mathbf S_n|=S\) is restored after every step, with periodic boundaries. For \(\alpha=0,\ T=0\), \(E\) and \(\sum_n S^z_n\) (rotation symmetry about \(z\)) must be conserved; the readout line tracks both. On a finite lattice \(k=2\pi m/L\) is quantised. With periodic boundaries the single-spin front reaches the antipode at \(t_s=(L/2)/v_{x,\max}\) (\(L/4JS\) for the FM); the readout line shows this time and warns when it is passed. After that the numerical–linear comparison remains valid (both solve the same ring), but infinite-system results such as the Bessel solution no longer apply. The absorbing boundary raises the damping smoothly in that layer, \(\alpha_n=\alpha+0.8\,f_n^2\) (\(f_n\) is the relative depth inside the layer, layer width \(\approx0.15L\)); reflection is small but not zero, energy is not conserved, and the plane wave is no longer an exact solution.
Unlike the FM, linear spin-wave theory is not exact for the AFM. In 1D the zero-point correction to the sublattice magnetisation diverges logarithmically, like \(\int dk/|\sin k|\): even at \(T=0\) the quantum chain has no Néel order. The exact (Bethe ansatz) lower edge for the \(S=1/2\) chain is the des Cloizeaux–Pearson spectrum \(\tfrac{\pi}{2}J|\sin k|\), and the excitations form a two-spinon continuum; for integer \(S\) a Haldane gap opens.
In the large-\(S\) limit and with at least some easy-axis anisotropy (\(K>0\)), the classical Néel picture is a good starting point; semiclassical corrections are of order \(1/S\). The simulation solves the classical equations exactly; it shows the classical limit, not the dynamics of the quantum chain.
In 2D the \(T>0\) magnetisation correction of the ferromagnet, \(\int d^2k\,/(e^{\beta\omega_k}-1)\), diverges logarithmically because \(\omega\propto k^2\): the isotropic model has no long-range order. Any easy-axis anisotropy (\(K>0\) or \(\lambda>1\)) opens a gap and restores order at \(T>0\); this is an Ising-like transition. In the easy-plane case a Berezinskii–Kosterlitz–Thouless transition is expected.
Unlike 1D, the 2D square-lattice \(S=1/2\) antiferromagnet has Néel order at \(T=0\). Linear spin-wave theory gives the sublattice magnetisation \(m\approx S-0.197\), consistent with the quantum Monte Carlo value \(m\approx0.307\), and the spin-wave velocity is renormalised only by a factor \(\approx1.18\).
In the code the model is defined by a single bond matrix \(B\) and a single-ion term; the ground state, effective field and dynamical matrix are all derived from them.
| Model | Linear dispersion | New physics |
|---|---|---|
| 1D ferromagnet | \(\Delta+2JS(1-\cos k)\) | Quadratic Goldstone mode, Bessel propagation |
| 1D antiferromagnet | \(h\pm2S\sqrt{(J+K)^2-J^2\cos^2k}\) | Two sublattices, linear dispersion, two chiral branches |
| XXZ / easy plane | \(2S\sqrt{J(1-\cos k)[J(1-\lambda\cos k)-K]}\) | Elliptical precession, linear Goldstone mode, canted and spin-flop states |
| Dzyaloshinskii–Moriya | \(\omega_k^{0}+2DS\sin k\) (FM, \(\mathbf M\parallel\mathbf D\)) | Non-reciprocal propagation, spiral instability, AFM dynamical instability |
| 2D square lattice | \(\Delta+4JS(\lambda-\gamma_{\mathbf k})\) | Square wave front, Γ–X–M path, checkerboard Néel |
| 3D cubic lattice | \(\Delta+2JS\sum_\mu(1-\cos k_\mu)\) | Γ–X–M–Γ–R–X path, G-type Néel, slice view |
| Finite temperature | \(\langle|\psi_{\mathbf k}|^2\rangle=2k_BT/\omega_{\mathbf k}\) | Stochastic LLG, equipartition, classical \(\Delta m\propto T\) |
| DMI spiral ground state | \(2S\sqrt{\tilde J(1-\cos k)[\tilde J-\lambda J\cos k-K]}\) | Conical spiral, phason mode, AFM spiral, locking of \(q\) |