following S $= 1$ Heisenberg Hamiltonian defined on a tetragonal lattice: $$\begin{aligned}
{\cal H} &=& - J \sum\limits_{ijk} \vec{S}_{ijk} \left( \vec{S}_{i+1 \, j k} + \vec{S}_{i \, j+1 \, k}
\right) - J' \sum\limits_{ijk} \vec{S}_{ijk}  \vec{S}_{ij \, k+1} \nonumber \\
&& - 
D \sum\limits_{ijk} \left( S^z_{ijk} \right)^2 - \sum\limits_{ijk} \vec{H} \vec{S}_{ijk} + 
{\cal H}_{nd}(J_{xz})\end{aligned}$$ with $\vec{S} = \left( S^x,S^y,S^z \right)$ and $\vec{H} = \left( H_x, 0, H_z \right)$. Here $i$, $j$, $k$ (corresponding to the a, b, c–directions respectively) label the lattice sites. The coupling in the layers, $J$, is ferromagnetic, whereas the coupling between adjacent layers, $J'$, may be antiferromagnetic ($J' < 0$) or ferromagnetic ($J' > 0$). The single-ion anisotropy, $D > 0$, is considered in order to obtain ground states with nonvanishing $z$–components of the spins. This may alternatively be achieved by considering instead an anisotropy in the exchange interactions (yielding three-dimensional XXZ models).

For the tetragonal models, I will consider two different couplings between planar and axial spin components: (1) a term involving only products of pairs of spins, being invariant under the symmetry operations of the point group $C_{4v}$ (this term can thus be considered to be analogous to the nondiagonal exchange term, see Equations 