we pick the axis without an anisotropy term. Second, we perform a $\pi$ rotation of spins around the $\hat{x}$ axis on one of the two sublattices—such transformation maps the anticipated antiferromagnetic ground state onto a ferromagnetic state. The result is the rotated spin Hamiltonian $$\begin{aligned}
 \label{eq:1}
   \mathcal{\tilde{H}}_J &= \sum\limits_{\left< {\bf i,j} \right>} J \: \left[ -\: S_{\bf i}^z S_{\bf j}^z + \frac{1}{2} \left( S_{\bf i}^+ S_{\bf j}^+ + S_{\bf i}^- S_j^- \right) \right] \nonumber \\
    &+ \gamma \sum\limits_{\bf i}\left(S_{\bf i}^y\right)^2 + \epsilon \sum\limits_{\bf i}\left(S_{\bf i}^x\right)^2.\end{aligned}$$

The next step is to utilize the Holstein-Primakoff transformation and the linear spin wave approximation using the assumption that the ground state is ferromagnetic and dressed with magnons. The thus obtained Hamiltonian is then diagonalized using the successive Fourier and Bogoliubov transformations. The resulting spin Hamiltonian reads $$\begin{aligned}
 \label{eq:1}
   {H}_J &= \sum_{\bf q}   {\Omega}_{\bf q} \beta^\dag_{\bf q} \beta_{\bf q} \; ,\end{aligned}$$ where the magnon energies $ {\Omega}_{\bf q}$ are given by $${\Omega}_{\bf q} = Y \sqrt{1 - \left( \frac{Z_{\bf q}}{Y} \right)^2},$$ with $$\label{hj solved}
  \begin{array}{l}
    Y = 4J+{\epsilon}+{\gamma},\\
    Z_{\bf q} = 4J{\gamma}_{\bf q}+{\epsilon}-{\gamma},\\
  \end{array}$$ with $$\gamma_{\bf q} = \frac{1}{2} \left( \cos(q_x) + \cos(q_y) \right).$$ The Bogoliubov coefficients $u_{\bf q}$ and $v_{\bf q}$, which describe the 