After substitution, changing to a rotating co-ordinate system and a bit of trigonometry we get:
$$ \begin{align} \mathcal{H} = & \sum_{i,j} J(\mathbf{d}) \left( S_i^\eta S_j^\eta + \sin (\mathbf{Q} \cdot \mathbf{d}_{i,j})(S_i^\mu S_j^\xi - S_i^\xi S_j^\eta \right) + \\ & \cos(\mathbf{Q}\cdot \mathbf{d}_{i,j}\left( S_i^\mu S_j^\mu + S_i^\xi S_j^\xi\right) \end{align}$$Spin state with a small field along the $z$ axis:
$$ \lvert n \rangle \equiv S, m = S -n $$Ladder operators for a harmonic oscillator:
$$ \begin{align} a^+ \lvert n \rangle & = \sqrt{n+1} \lvert n + 1 \rangle \\ a \lvert n \rangle & = \sqrt{n} \lvert n -1 \rangle \\ \left[ a, a^+ \right] & = 1 \end{align}$$Ladder operators for spin:
$$ \begin{align} S^- \lvert n \rangle & = \sqrt{2S}\sqrt{1 - \frac{n}{2S}}\sqrt{n+1} \lvert n + 1 \rangle \\ S^+ \lvert n \rangle & = \sqrt{2S}\sqrt{1 - \frac{n - 1}{2S}}\sqrt{n} \lvert n -1 \rangle \\ \end{align}$$Assuming $\langle n \rangle \ll S$ the coupled harmonic oscillator is a good model. The excited states of the harmonic oscillator are magnons.
Spin wave approximation
Ladder operators for spin:
$$ \begin{align} S^- & = \sqrt{2S} a^+ \hat{f} \\ S^+ & = \sqrt{2S} \hat{f} a \end{align}$$Where:
$$ \hat{f} = \sqrt{1 - \frac{n}{2S}} $$Applying the transformation, $S$ in boson creation and annihilation operators:
$$ \begin{align} S_i^\eta & = \frac{1}{2}(S^+ + S^- ) = \sqrt{S/2}\hat{f}(a^+_i + a_i) \\ S_i^\mu & = -\frac{i}{2}(S^+ - S^- ) = \sqrt{S/2}\hat{f}(a^+_i - a_i) \\ S_i^\xi & = S - n_i = S - a_i^+a_i \end{align}$$So:
$$ S_i^\eta S_j^\eta = S/2 \hat{f}_i \hat{f}_j (a^+_i + a_i)(a^+_j + a_j) $$etc...
We have a hamiltonian of boson operators which can be expanded in powers of $1/S$ to obtain:
$$ \mathcal{H} = E_0 + \mathcal{H}_1 + \mathcal{H}_2 + \mathcal{H}_3 + \mathcal{H}_4 ... $$Where $E_0$ is the classical result which is derived in many textbooks:
$$ E_0 = S^2 \sum_{i, \mathbf{d}} J(\mathbf{d})\cos (\mathbf{Q} \cdot \mathbf{d}) $$1 operator term:
$$ \mathcal{H}_1 = S^{3/2} \sum_{i,j} \frac{i}{2} J(\mathbf{d}_{i,j}) \sin (\mathbf{Q} \cdot \mathbf{d}_{i,j})(a_i^+ - a_i + a_j - a_j^+) $$2 operator term:
$$\begin{align} \mathcal{H}_2 = S \sum_{i,j} \frac{1}{2} J(\mathbf{d}_{i,j}) &( (1 - \cos (\mathbf{Q} \cdot \mathbf{d}_{i,j}))(a_ia_j + a_ia_j) \\ & + (1 + \cos (\mathbf{Q} \cdot \mathbf{d}_{i,j}))(a_ia_j^+ + a_i^+a_j) \\ & - 2\cos (\mathbf{Q} \cdot \mathbf{d}_{i,j})(a_i^+a_i + a_j^+a_j) ) \end{align}$$3 operator term:
4 operator term:
Luckily, $\mathcal{H}_2$ can be written in matrix form:
$$\mathcal{H}_2 = \sum_\mathbf{k} \mathbf{x}^\dagger H(\mathbf{k}) \mathbf{x} $$Where $\mathbf{x}$ is a vector of Boson operators:
$$ \mathbf{x} = \begin{bmatrix} a_\mathbf{k} \\ a_{-\mathbf{k}}^+\end{bmatrix}$$And the matrix of the Hamiltonian has the form:
$$ H = \begin{bmatrix} A & B \\ B & A \end{bmatrix} $$In Bogoliubov method, we define new operator $b$ with the following transformation:
$$ \begin{align} b & = ua + va^+ \\ b^+ & = ua^+ +va \end{align}$$The new operator has to fulfill the commutation relations:
$$ \begin{align} [b, b^+] & = 1 \\ u^2 + v^2 & = 1 \end{align}$$With the right parameter choice:
$$ \mathcal{H}_2 = \sum_\mathbf{k} \omega_\mathbf{k} \left( b^+b + \frac{1}{2} \right) $$And the spin wave dispersion:
$$ \omega_\mathbf{k} = \sqrt{A^2 - B^2} $$The above calculation is equivalent to solving the eigenvalue problem of $gH$, where $g = [x, x^\dagger]$ commutator matrix.
Now we have $\omega_\mathbf{k}$ and half the story. Remember for neutron scattering:
With the correlation function: