= \left[ i |1\ra - i |\bar{1}\ra \right]/\sqrt{2}, ~ |y\ra = \left[ |1\ra + |\bar{1}\ra \right]/\sqrt{2},~~ |z\ra = -i |0\ra,$ where we abbreviate $|S^z = \pm1\ra \equiv |\pm1\ra$, $|S^z = 0\ra \equiv |0\ra$, and $|\bar{1}\ra \equiv |-1\ra$, we can introduce the complex site-factorized wavefunction vectors at each site to characterize any possible ordered state with short-ranged correlations as $ {\bf d}_i =  ( d^x_i, d^y_i, d^z_i  )$ in the $\{|x\ra, |y\ra, |z\ra\}$ basis (${\bf d}_i = {\bf u}_i + i {\bf v}_i$). The normalization of the wavefunction leads to the constraint ${\bf d}_i \cdot \bar{\bf d}_i = 1$, or equivalently, ${\bf u}^2_i + {\bf v}^2_i =1$, and the overall phase can be fixed by requiring ${\bf d}_i^2 = \bar{\bf d}_i^2$, i.e., ${\bf u}_i \cdot {\bf v}_i = 0$. In terms of ${\bf d}$, the Hamiltonian  is expressed as $$\begin{aligned}
&&H = \sum_{i, \delta_n} \left[J_n\left| {\bf d}_i \cdot \bar{\bf d}_j\right|^2
+  \left(K_n - J_n\right) \left|{\bf d}_i \cdot {\bf d}_j \right|^2 + K_n\right]~~~~ \nonumber \\ 
&& - \sum_{\triangpic}  \mathfrak{R}_3 (\bf{d}_1 \cdot \bar{\bf{d}}_3) ( \bf{d}_2 \cdot \bar{\bf{d}}_1) ( \bf{d}_3 \cdot \bar{\bf{d}}_2) + \Hc \nonumber \\ 
&&  + \sum_{\rhombpic} \mathfrak{R}_4 (\bf{d}_1 \cdot \bar{\bf{d}}_4) (\bf{d}_2 \cdot \bar{\bf{d}}_1)(\bf{d}_3 \cdot \bar{\bf{d}}_2) (\bf{d}_4 \cdot \bar{\bf{d}}_3) + \Hc. \label{Hd_4ring}~~~~~\end{aligned}$$ 