any of the three models can be written as a sum of local Hamiltonians $G_{k}^{(p)}$ of the $p$-sided polygonal shape, in particular, $${\cal H} = \sum\limits_{{\langle k \rangle}_p}^{~} G_{k}^{(p)}\, ,
\label{Hm1}$$ where $k$ labels the polygons and the sum runs over the set of all indices of the lattice polygons ${\langle k \rangle}_p$. The local Hamiltonian takes the form $$\begin{aligned}
\nonumber
G_{k}^{(p)} & = & -\frac{1}{2}\sum\limits_{i=1}^{p} \bigg[
      J_{xy} \left( S_{k_i}^{x} S_{k_{i+1}}^{x} + S_{k_i}^{y} S_{k_{i+1}}^{y} \right) \\
      & + & J_{z} S_{k_i}^{z} S_{k_{i+1}}^{z}
      + \frac{h}{4} \left( S_{k_i}^{x} + S_{k_{i+1}}^{x}\right) \bigg] \, ,
\label{Hm2}\end{aligned}$$ where $k_1,k_2,\dots,k_p$ label the spin positions on the $k^{\rm th}$ $p$-sided polygon (noticing that $k_{p+1} \equiv k_1$), and $S_{k_i}^x$, $S_{k_i}^y$, $S_{k_i}^z$ denote the corresponding Pauli spin-$\frac{1}{2}$ operators. We consider constant nearest-neighbor couplings $J_{xy}, J_{z}$ and a uniform external magnetic field $h$. By setting $J_{xy} = 0$ and $J_{z} = 1$ we obtain the TFIM at the transverse magnetic field $h$, whereas the choice $J_{xy} = 1,J_{z} = h = 0$ gives the XY model and $J_{xy} = - J_{z} = 1, h = 0$ the Heisenberg model  [@DaniskaGendiar].

Our task is to calculate an approximate ground-state of the system in the thermodynamic limit in the product form $$| 