element. From the single-site wavefunction the SF order parameter and lattice occupancy can be calculated from the expressions in Eqns. (\[gw\_phi2s\_1\]-\[gw\_phi2s\_2\]) and (\[num2s\_1\]-\[num2s\_2\]).

Mean-field decoupling theory
----------------------------

To calculate the phase boundaries between MI and SF phases analytically we use the site decoupled mean-field theory [@oosten_01; @iskin_09; @bandyopadhyay_19]. For this, we adapt perturbative analysis of the mean-field Hamiltonian in Eq. (\[ham\_ss\_tbec\]). It is important to note that the SF order parameter $\phi^{k}_{p,q} $ is zero in the MI phase, but nonzero in the SF phase. So, the vanishing of the SF order parameter $\phi^{k}_{p,q} \rightarrow 0^{+}$ marks the MI-SF phase boundary in the phase diagram. With this consideration, in the site decoupled mean-field theory, the interaction and the chemical potential terms constitute the unperturbed Hamiltonian $\hat{h}_{p,q,0}^{\rm TBH}$. From Eq. (\[ham\_ss\_tbec\]), $$\begin{aligned}
\hat{h}_{p,q,0}^{{\rm TBH}} 
            &= & \sum_{k}\left[ \frac{U_{kk}}{2}\hat{n}_{p, q}^k
                  \left (\hat{n}_{p, q}^k-1\right) 
                  - \tilde{\mu}^k_{p,q}\hat{n}_{p, q}^k\right]
                  \nonumber\\
                  &&+ U_{12}\hat{n}_{p, q}^1 \hat{n}_{p, q}^2,
\label{ham_ss_h0}\end{aligned}$$ which is diagonal with respect to the Fock basis states.Then, the hopping terms in Eq. (\[ham\_ss\_tbec\]) as the perturbation, $$\begin{aligned}
\hat{h}_{p,q,1}^{{\rm TBH}} &= &- \sum_{k}\Bigl [J_x^k 
                  \left(\hat{b}_{p+1, q}^{\dagger k}\phi_{p,q}^k 
                  + \phi^{k*}_{p + 1, q}\hat{b}_{p, q}^k \right ) + {\rm H.c.}
                   \nonumber\\
                  && +  J_y^k\left(\hat{b}_{p, q+1}^{\dagger k} \phi_{p,q}^k 
                  +  \phi^{k*}_{p, 