single-site mean field Hamiltonian, the total Hamiltonian of the system is $$\hat{H}^{{\rm TBH}} = \sum_{p,q}\hat{h}_{p,q}^{{\rm TBH}}.$$ For the details of the derivations see [@bai_18]. To get the ground state we digonalize the Hamiltonian in Eq. (\[ham\_ss\_tbec\]) at each site. And, for this we use the Gutzwiller ansatz, based on which the ground state at the ($p,q$)th site is [@anufriiev_16] $$|\psi\rangle_{p,q} = \sum_{n_1, n_2 }c^{(p,q)}_{n_1, n_2} 
                       |n_1, n_2\rangle_{p,q}.
  \label{gw_2s}$$ Here $|n_1, n_2\rangle$ is a Fock state, which is direct product of the $n_1$ and $n_2$ occupation number states of the 1st and 2nd species, respectively. The occupation number states $n_k\in [0,N_b -1]$ where $N_b$ is the total number of local Fock states used in the computation and $c_{n_1, n_2}^{p,q}$ are complex co-efficients with $\sum_{n_1,n_2} |c^{(p,q)}_{n_1,n_2}|^2$ = 1. From the ground state, we can compute the new SF order parameter of the two species as

$$\begin{aligned}
    \phi_{p, q}^1& =& _{p,q}\langle\psi|\hat{b}_{p, q}^1|\psi\rangle_{p,q} 
            = \sum_{n_1, n_2}\sqrt{n_1} 
              {c^{(p,q)*}_{n_1-1, n_2}}c^{(p,q)}_{n_1,n_2},
              \label{gw_phi2s_1}\\
    \phi_{p, q}^2& =& _{p,q}\langle\psi|\hat{b}_{p, q}^2|\psi\rangle_{p,q} 
            = \sum_{n_1, n_2}\sqrt{n_2} 
              {c^{(p,q)*}_{n_1, n_2-1}}c^{(p,q)}_{n_1,n_2}.
    \label{gw_phi2s_2}              
  \end{aligned}$$

Similarly, corresponding the lattice occupancies are

$$\begin{aligned}
   \rho_{p, q}^1& =& _{p,q}\langle\psi|\hat{n}_{p, q}^1|\psi\rangle_{p,q} 
            = \sum_{n_1, n_2} n_1 |c^{(p,q)}_{n_1,n_2}|^2,
            \label{num2s_1}\\
   \rho_{p, q}^2& =& _{p,q}\langle\psi|\hat{n}_{p, q}^2|\psi\rangle_{p,q} 
            = \sum_{n_1, n_2}n_2|c^{(p,q)}_{n_1,n_2}|^2. 
   \label{num2s_2}              
  \end{aligned}$$

Using the new SF order 