cells are connected by ${\bf e}_1 \equiv \hat{x}$ and ${\bf e}_2 \equiv \hat{y}$.

For clarity, we introduce $ D^T_{\alpha = y,z} ({\bf k}) = \{d_\alpha({\bf k},1), d_\alpha({\bf k},2),\cdots,d_\alpha({\bf k},15),d_\alpha({\bf k},16)\}$, and $A_\alpha({\bf k}) = \left\{ D^T_\alpha({\bf k}), D^\dagger_\alpha(-{\bf k})\right\}$ for clarity below. Within the linear flavor wave theory calculation, we find that the Hamiltonian is $H_\mu = H_c + H^\mu_B$, where $\mu =$ BC or PL, and $ H_c = 32\sum_{\bf k} \left[ J_1 + J_2 - \frac{13}{4} J_3 + K_1 + K_2 + 3 K_3 \right]$, and can be determined straightforwardly independent of the boson fields and independent of BC or PL order. The $H^\mu_B \equiv \sum_{n=1,2,3}H^\mu_n$ represent the boson Hamiltonian related to the $n$th neighbor couplings that are different for $\mu=BC,PL$. After Fourier transform, we find that $H^\mu_n = \sum_{{\bf k}, \eta = y,z} (A^\mu_\eta)^\dagger \mathcal{H}^\mu_{n,\eta} A^\mu_\eta $, with $$\begin{aligned}
\mathcal{H}^\mu_{n,\eta} = \begin{pmatrix}
\alpha^\mu_{n,\eta} & \gamma^\mu_{n,\eta} \\
(\gamma^\mu_{n,\eta})^\dagger & \alpha^\mu_{n,\eta}
\end{pmatrix},\end{aligned}$$ where $\alpha^\mu_{n,\eta}$ and $\gamma^\mu_{n,\eta}$ are $16\times16$ Hermitian matrices. The $H^\mu_n$ can be straightforwardly diagonalized to get the energies of BC and PL, which can be used to estimate the finite stiffness regimes of BC and PL shown in the main texts. Below we list the matrix elements 