of the state since we are only interested in the energy and the system is homogeneous. With this, the energies at $t_d=0$ can be written as $$\begin{aligned}
 E(nL-1)  &= (L-1) E^-_n +  E^-_{n-1} \\
 E(nL)      &=L E^-_n \\
 E(nL+1) &=(L-1) E^-_n +  E^-_{n+1},\end{aligned}$$ and the chemical potentials evaluate to $$\begin{aligned}
\mu^+_n&=E^-_{n+1}-E^-_n \\
 &= \omega - \frac {\chi_{n+1}} 2 +(1-\delta_{n0}) \frac {\chi_{n}} 2+\delta_{n,0} \frac{\Delta-\omega}2,\\
 \intertext{for any $n$ and}
 \mu^-_n&= E^-_n-E^-_{n-1} \\
 &= \omega - \frac {\chi_{n}} 2 +(1-\delta_{n1}) \frac {\chi_{n-1}} 2+\delta_{n,1} \frac{\Delta-\omega}2,\end{aligned}$$ for $n>0$. Thus for a commensurate number of excitations the system displays particle-hole gaps. Since $\mu^-_{n+1}=\mu^+_n$ the chemical potential for non-commensurate total number of excitations between $N=nL$ and $N=(n+1)L$ is the same, corresponding to a critical point. For non-vanishing tunneling the critical points extend to critical regions.\
The simplest numerical method to obtain a qualitative phase diagram is the so called mean field theory. As described for instance in [@Greentree2006; @Hartmann2008a; @Lei2008; @Makin2008], mean field theory can be implemented by introducing an order parameter $\Psi$, which in our case is chosen to be homogeneous and real valued. Decoupling the hopping term by using $$\hat a_j^\dagger \hat a_l \mapsto \Psi  \left(\hat a_j^\dagger + \hat a_l\right) - \left|\Psi\right|^2,$$ the whole JCH 