between the localised $S=1$ magnetic moments. The relevant Hamiltonian is well-known in this case and reads [@Kunkemoeller2015; @Kunkemoeller2017; @Jain2017], $$\label{model1}
  {\mathcal H}_J = J \sum\limits_{\langle {\bf i}, {\bf j} \rangle} \mathbf{S_{\bf i}} \cdot \mathbf{S_{\bf j}} + \epsilon \sum\limits_{\bf i} \left( S_{\bf i}^z \right)^2 + \gamma \sum\limits_{\bf i} \left( S_{\bf i}^x \right)^2,$$ where the summation runs over all nearest-neighbour pairs on a 2D square lattice, $J$ is the spin exchange constant, and $\mathbf{S_{\bf i}}$ are the spin $S=1$ operators. As already discussed in Refs.  the spin model is highly anisotropic, with the suggested values of the $\hat{z}$ ($\hat{x}$) axis anisotropy being equal to $\epsilon = 5.6 J$ ($\gamma = 0.25 J$), respectively, with $J =  5.6$ meV reproducing the spin wave dispersion observed in the inelastic neutron scattering experiment [@Kunkemoeller2015; @Kunkemoeller2017; @Jain2017]. We note that such a large spin anisotropy originates in the large spin-orbit coupling on the ruthenium ions which, however, is quite widely considered as not strong enough to stabilise the $S=0$ ground state [@Mizokawa2001; @Kunkemoeller2015; @Kunkemoeller2017; @Das2018; @Pincini2019]. Although the latter result can naively be understood as a consequence of the crystal field splitting (between the $xz, yz$ and the $xy$ orbitals) being about twice larger than the 