Hamiltonian (\[eq:JCHreal\]) in the grand-canonical ensemble uncouples in real space with a local Hamiltonian being $$\begin{gathered}
 \hat H^{\rm MF} = \left(\omega-\mu\right)  \hat a^\dagger \hat a+ \left(\Delta-\mu\right) \hat \sigma^+\hat\sigma^- + g\left(\hat a^\dagger \hat\sigma^- + \hat a\hat \sigma^+\right)\\
 - 2 \widetilde J \Psi  \left(\hat a^\dagger + \hat a\right) + 2\widetilde J \left|\Psi\right|^2.\label{eq:JCHMeanField}\end{gathered}$$ At this point, we omitted the spatial index because the problem is purely local. The modified hopping amplitude $\widetilde J=-\sum_d t_d$ gives the effective coupling within the mean field scheme. The phase diagram is now found by diagonalizing the mean field Hamiltonian (\[eq:JCHMeanField\]) either exactly by means of perturbation theory or numerically, setting an upper bound for the maximal number of bosonic excitations in the system. The ground state energy is then given by $\underset{\Psi}{\min}\  E[\Psi]$ and the MI is distinguished from the SF by the value of $\Psi$ for the minimal energy. For $\Psi\equiv0$, the system is in a MI state, for $\Psi>0$, the ground state is superfluid. This sets the point of the MI to SF transition. It should be mentioned at this point, that this method gives inadequate results in one dimension ($D=1$) but is exact for $D\to\infty$. Additionally, the effective hopping $\widetilde J$ must be larger 