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_{jin 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 Jbegin{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} (2delta m_i^{y} + delta m_j^{y} ) -frac{1}{2} (2delta 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} (2delta m_i^{y} + delta m_j^{y}) – frac{sqrt{3}}{2} (2 delta m_i^{x} + delta m_j^{x})
– frac{sqrt{3}}{2} (2delta m_i^{y} + delta m_j^{y}) – frac{1}{2} (2delta 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 Jbegin{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} = (Xpm 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.