There is a paper (PhysRevB.95.014435) in which the dispersion relation for some Heisenberg model on the honeycomb lattice is derived from the Landau-Lifshitz equation: \begin{align} \frac{d S_i}{dt} = - S_i \times \mathcal H_{\rm eff} \end{align} Their attempt from Eq. 2 to Eq.4 is pretty simple and I'll try the same for the 2D triangular Heisenberg antiferromagnet (THAF) (in xy-plane), which has a much simpler Hamiltonian: \begin{align} \mathcal H = \sum_{\langle {ij}\rangle } J S_i S_j,\quad \mathcal H_{\rm eff} = J \sum_j S_j \end{align} where $\langle {ij}\rangle$ sums over all nearest neighbors. There are some papers out there (for example PhysRevB.74.180403) which have derived the dispersion to be \begin{align} \omega_{\bf k} = \sqrt{(1- \gamma_{\bf k} ) ( 1+ 2 \gamma_{\bf k} ) } \label{eq:thaf_disp} \end{align} with \begin{align} \gamma_{\bf k} = \frac{1}{z} \sum_{j} \mathrm{e}^{i \bf{k}( \bf{R}_i - \bf{R}_j )} = \frac{1}{3}\left(\cos k_{x}+2 \cos \frac{k_{x}}{2} \cos \frac{\sqrt{3}}{2} k_{y}\right) \, . \end{align} The ground-state of the THAF is the $120^{\circ}$-Neel order. My idea is similar to the derivation in Linear Spin Wave Theory and I'm starting by some rotation of spin vectors \begin{align} S_{i \in A} &= (\delta m_i^{x}, \delta m_i^{y}, 1) \\ S_{i \in B } &= ( \sqrt{3}/2 \delta m_i^{y} - 1/2 \delta m_i^{x}, -\sqrt{3}/2 \delta m_i^{x} - 1/2 \delta m_i^{y}, 1) \\ S_{i \in C} &= ( -\sqrt{3}/2 \delta m_i^{y} - 1/2 \delta m_i^{x}, \sqrt{3}/2 \delta m_i^{x} - 1/2 \delta m_i^{y}, 1) \end{align} where A,B,C are the three sublattices of the ground-state and $\delta m \ll 1$ . Then I tried to solve the Landau-Lifshitz equation: \begin{align*} \frac{d S_{i \in A}}{dt} &=- \begin{pmatrix} \delta m_i^{x} \\ \delta m_i^{y} \\ 1 \end{pmatrix} \times \left(\sum_j J S_{j\in B} + J S_{j \in C}\right) =- \sum_j J \begin{pmatrix} \delta m_i^{x} \\ \delta m_i^{y} \\ 1 \end{pmatrix} \times \begin{pmatrix} - \delta m_j^{x} \\ - \delta m_j^{y} \\ 2 \end{pmatrix} \approx - \sum_jJ \begin{pmatrix} \delta m_j^{y} + 2 \delta m_i^{y} \\ - \delta m_j^{x} - 2 \delta m_i^{x} \\ 0 \end{pmatrix} \\ \frac{d S_{i \in B}}{d t} &= -\begin{pmatrix} \frac{\sqrt{3}}{2} \delta m_i^{y} - \frac{1}{2}\delta m_i^{x} \\ -\frac{\sqrt{3}}{2} \delta m_i^{x} - \frac{1}{2} \delta m_i^{y} \\ 1 \end{pmatrix} \times \left(\sum_j J S_{j \in A} + J S_{j \in C} \right) \\ &= - \sum_j J \begin{pmatrix} \frac{\sqrt{3}}{2} \delta m_i^{y} - \frac{1}{2} \delta m_i^{x} \\ -\frac{\sqrt{3}}{2} \delta m_i^{x} - \frac{1}{2} \delta m_i^{y} \\ 1 \end{pmatrix} \times \begin{pmatrix} \frac{1}{2} \delta m_j^{x} - \frac{\sqrt{3}}{2} \delta m_j^{y} \\ \frac{\sqrt{3}}{2} \delta m_j^{x} + \frac{1}{2} \delta m_j^{y} \\ 2 \end{pmatrix} \approx - \sum_j J \begin{pmatrix} -(\sqrt{3} \delta m_i^{x} + \delta m_i^{y}) - ( \frac{\sqrt{3}}{2} \delta m_j^{x} + \frac{1}{2} \delta m_j^{y} ) \\ \frac{1}{2} \delta m_j^{x} - \frac{\sqrt{3}}{2} \delta m_j^{y} - (\sqrt{3} \delta m_i^{y} - \delta m_i^{x}) \\ 0 \end{pmatrix} \\ &=\sum_j J\begin{pmatrix} \frac{\sqrt{3}}{2} (2 \delta m_i^{x} + \delta m_j^{x} ) + \frac{1}{2}(2 \delta m_i^{y} +\delta m_j^{y} ) \\ \frac{\sqrt{3}}{2} (2\delta m_i^{y} + \delta m_j^{y} ) -\frac{1}{2} (2\delta m_i^{x} + \delta m_j^{x} ) \\ 0 \end{pmatrix} \\ \frac{d S_{i \in C}}{d t} &= - \sum_j \begin{pmatrix} -\frac{\sqrt{3}}{2} \delta m_i^{y} - \frac{1}{2} \delta m_i^{x} \\ \frac{\sqrt{3}}{2} \delta m_i^{x} - \frac{1}{2} \delta m_i^{y} \\ 1 \end{pmatrix} \times \begin{pmatrix} \frac{\sqrt{3}}{2} \delta m_j^{y} + \frac{1}{2} \delta m_j^{x} \\ -\frac{\sqrt{3}}{2} \delta m_j^{x} + \frac{1}{2} \delta m_j^{y} \\ 2 \end{pmatrix} \approx - \sum_j J \begin{pmatrix} \sqrt{3} \delta m_i^{x} - \delta m_i^{y} - (-\frac{\sqrt{3}}{2} \delta m_j^{x} + \frac{1}{2} \delta m_j^{y}) \\ (\frac{\sqrt{3}}{2} \delta m_j^{y} + \frac{1}{2} \delta m_j^{x}) + \sqrt{3} \delta m_i^{y} + \delta m_i^{x} \\ 0 \end{pmatrix} \\ &= \sum_j J \begin{pmatrix} \frac{1}{2} (2\delta m_i^{y} + \delta m_j^{y}) - \frac{\sqrt{3}}{2} (2 \delta m_i^{x} + \delta m_j^{x}) \\ - \frac{\sqrt{3}}{2} (2\delta m_i^{y} + \delta m_j^{y}) - \frac{1}{2} (2\delta m_i^{x} + \delta m_j^{x}) \\ 0 \end{pmatrix} \end{align*}
By using Bloch-Theorem: \begin{align} \delta m_i^{x} = X \exp(i \left( \bf{k} \bf{R}_i - \omega t \right) ), \quad \delta m_i^{y} = Y \exp(i \left( \bf{k} \bf{R}_i - \omega t \right) ) \end{align} Since I only have now one sublattice I don't need $X_A$, $X_B$ and $X_C$ etc. like in the paper. If you compare left-hand and right-hand side of the those equations of motions all do have the same structure. This structure looks like
\begin{align} i \omega \begin{pmatrix} X \\ Y \end{pmatrix} \mathrm{e}^{i (\bf{k} \bf{R}_i - \omega t)} = \sum_j J \begin{pmatrix} - 2 Y \\ 2X \end{pmatrix}\mathrm{e}^{i (\bf{k} \bf{R}_i - \omega t)} + \sum_j J\begin{pmatrix} -Y \\ X \end{pmatrix} \mathrm{e}^{i (\bf{k} \bf{R}_j - \omega t)} \end{align} where the Bloch theorem is already used. This would then lead to the following matrix \begin{align} i \omega \begin{pmatrix} X \\ Y \end{pmatrix} = J \begin{pmatrix} 0 & -2 - \gamma_k \\ 2 + \gamma_k & 0 \end{pmatrix} \begin{pmatrix} X \\ Y \end{pmatrix} = H \begin{pmatrix} X \\ Y \end{pmatrix} \end{align} The paper sugested using $\psi^{\pm} = (X\pm iY)/\sqrt{2}$. This can be achieved by the Matrix \begin{align} U = \begin{pmatrix} 1 & i \\ 1 & -i \end{pmatrix} \end{align} and by calculating $i/2 \sigma_z UHU^{-1}$ I ended up with an hermitian matrix which uses $\psi^{\pm}$ as the amplitudes like sugested in the paper above: \begin{align} \begin{pmatrix} - \gamma_k - 2 & 0 \\ 0 & \gamma_k + 2 \end{pmatrix} \end{align} which would lead to $\omega_k = \pm \sqrt{(\gamma_k + 2)^2}$ which is obviously wrong but I cannot figure out where my mistake is or where I'm thinking wrong.