Landau-Lifshitz Equation of Motion for Triangular Heisenberg Antiferromagnet 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.

I see two possible problems in your consideration.

*

*You've investigated perturbations of ferromagnetic ground state. When spin variations $\delta m$ are zeros, spins on three sublattices are the same:
$$
S_i = (0, 0, 1),\quad \forall i.
$$


*The Landau-Lifshitz equation is a nonlinear one. Effective field ${\cal H}_{i,{\rm eff}}$ depends on neighboring spins. Hence you need to take into account variations of effective field:
$$
\frac{d \delta S_i}{dt} = -\delta S_i \times {\cal H}_{i,{\rm eff}} - S_i \times \delta {\cal H}_{i,{\rm eff}}.
$$
I didn't analyze your application of the Bloch theorem. I think there also could be problems. Neel state on triangular lattice is invariant under translation of states of triangular cells of spins, not of individual spins.
