JCH model, we will follow the usual route. Since the total number of excitations in the system $$\hat N=\sum_{   j} \left( \hat a_{   j}^\dagger\hat a_{   j}+ \hat\sigma^+_{   j} \hat\sigma^-_{   j}\right)$$ commutes with the full Hamiltonian (\[eq:JCHreal\]), it is enough to treat the system for a fixed number of excitations. The boundary of the $n$th Mott lobe can then be determined by calculating the total energy $E(N)$ for $N=nL-1$, $N=nL$ and $N=nL+1$ excitations in a system with $L$ sites. The chemical potential then reads $$\mu^\pm_n = \pm\Bigl[ E({n L\pm1})- E({n L})\Bigr] ,\label{eq:chemicalPotential}$$ where the plus sign belongs to the upper boundary of the Mott lobe and the minus sign to the lower one. For $t_d\equiv0$ , $\mu^ \pm_n$ can be calculated straight forwardly. Starting with the energy for $N=nL$ excitations with $n$ being an integer, i.e. for a commensurate number of excitations, it can be seen that due to the nonlinear dependence of the single-site energy $E^-_n$ on $n$, the excitations will distribute equally over the whole lattice. The ground state is therefore given by $\vec n=\{n,n,\dots,n\}$. Now, when adding (removing) a single excitation from the whole system, the ground state is given by $\{n\pm1,n,\dots,n\}$, where we have ignored the degeneracy 