\fl \hat{H} = - N \left( \Omega S + \frac{J_\mathrm{an}}{4} \right) + \sum_{{\mathbf{k}}} \Gamma_{{\mathbf{k}}} \hat{b}^{\dag}_{{\mathbf{k}}} \hat{b}_{{\mathbf{k}}} + \sum_{{\mathbf{k}}} \frac{\Delta_{{\mathbf{k}}}}{2} \left( \hat{b}_{{\mathbf{k}}} \hat{b}_{-{\mathbf{k}}} + \hat{b}^{\dag}_{{\mathbf{k}}} \hat{b}^{\dag}_{-{\mathbf{k}}} \right),\end{aligned}$$ where $\gamma_{{\mathbf{k}}} = 2 \cos (k_{x}a) + 2 \cos (k_{y}a) + 2 \cos (k_{z} a)$, $N$ is the total number of lattice sites, and $$\begin{aligned}
\Delta_{{\mathbf{k}}} & = \frac{J_{1}}{2} \gamma_{{\mathbf{k}}} + \frac{J_\mathrm{an}}{2}, \\
\Gamma_{{\mathbf{k}}} & = \Omega - \Delta_{{\mathbf{k}}}.\end{aligned}$$ Note that we have set $E = 0$ here, since the phonon spectrum is measured at zero field.

It is convenient to formulate the dynamics of the pseudomagnons using Green’s functions. The Green’s functions are correlation functions between the pseudomagnon creation and annihilation operators, and the equations of motion of the Green’s functions therefore include the equations of motion of $\hat{b}_{{\mathbf{k}}}$ and $\hat{b}^{\dag}_{{\mathbf{k}}}$. The spin-wave excitation spectrum can then be obtained from the poles of the Green’s function.

The Green’s function and its equation of motion are, respectively, $$\label{G1}
D_{1}({\mathbf{k}}, t) = -i \left\langle \left[ \hat{b}_{{\mathbf{k}}}(t), \hat{b}^{\dag}_{{\mathbf{k}}}(0) \right] \right\rangle \theta (t),$$ $$\label{G1_eom}
\frac{dD_{1}({\mathbf{k}}, t)}{dt} = -i \delta (t) - i \Gamma_{{\mathbf{k}}} D_{1}({\mathbf{k}}, t) - i \Delta_{{\mathbf{k}}} D_{2}({\mathbf{k}}, t),$$ where $\theta(t)$ is the step function. The second Green’s function that appears in equation  and its equation 