b_{i\beta}\right],\end{aligned}$$ where $j = i + \delta_n,$ and $\delta_n$ (with $n = 1,~2,~3$) connects site $i$ to its $n$th nearest neighbor sites, and we ignore the constant terms $K_n - J_n$ in the above equation. For performing flavor wave theory calculation, we introduce different local rotations for different site $j$. For site $j \in | S^z = \pm1\ra$, we introduce $$\begin{aligned}
\begin{pmatrix}
b_{jx} \\
b_{jy} \\
b_{jz}
\end{pmatrix} =
\begin{pmatrix}
\mp \frac{i}{\sqrt{2}} & \frac{1}{\sqrt{2} }& 0\\
\frac{1}{\sqrt{2}} & \mp \frac{i}{\sqrt{2}} & 0 \\
0 & 0 & 1
\end{pmatrix}
\begin{pmatrix}
d_{jx}\\
d_{jy}\\
d_{jz}
\end{pmatrix},\end{aligned}$$ which still preserves the local constraint of Eq. (\[Eq:constraint\]) with $b_{i\alpha} \rightarrow d_{i\alpha}$. At each site only one flavor of bosons $d_{ix}$ condenses, and we replace $d^\dagger_{ix}$ and $d_{ix}$ by $(M - d_{iy}^\dagger d_{iy} - d^\dagger_{iz} d_{iz})^{1/2}$, with $M=1$ in the present case. A $1/M$ expansion up to the quadratic order of the Holstein-Primakoff bosons $d_y$ and $d_z$ followed by an appropriate transformation allows us to extract the ground state energy. We will replace the labeling $d_{i\alpha} = d_{\alpha}({\bf r},a)$, where ${\bf r}$ runs over the Bravais lattice of unit cells of the square network and $a=1,2,3,4,\cdots, 15,16$ runs over the sub lattices. The different unit 