# Variational Principle to find the Mean Field Bose-Hubbard Ground State

I'm trying to find the ground state energy of the Bose-Hubbard in the context of a mean-field approximation. Assuming small fluctuations, the mean field Hamiltonian can be decoupled into the sum of single site terms, each given by $$\hat H = \frac U2 \hat n (\hat n-1) - \mu \hat n - \psi(\hat a + \hat a^\dagger)+ \psi^2$$ where $U,\mu \geq 0$ are the interaction potential and chemical potential respectively and $\psi = \langle \hat a \rangle = \langle \hat a^\dagger\rangle$ is the parameter we will use to minimize the ground state energy.

I suppose there are a total of $N$ bosons, and the unnormalized single site ground state is given by $$| \Phi_0 \rangle = \sum_{j=0}^N c_j | j\rangle \qquad c_j\in \mathbb R$$ in the Fock basis. Then the ground state energy is given by $$E_0 = \frac{\langle \Phi_0 |\hat H|\Phi_0\rangle}{\langle \Phi_0 |\Phi_0\rangle} = \left( \frac{1}{\sum_{\ell=0}^N c_\ell^2}\right) \sum_{k=0}^N \left[ c_k^2 \left(\frac U2 k(k-1)-\mu k + \psi^2 \right) - 2\psi c_k c_{k+1} \sqrt{k+1}\right]$$ where $c_{N+1}:=0$. To find the ground state, we must minimize $E_0$ with respect to each of the parameters:

$$0 = \frac{\partial E_0}{\partial c_j}= \left( \frac{1}{\sum_{\ell=0}^N c_\ell^2}\right) \left[ 2 c_j \left(\frac U2 j(j-1)-\mu j + \psi^2 \right) - 2\psi c_{j+1} \sqrt{j+1} - 2\psi c_{j-1} \sqrt{j}-2 c_j E_0\right]$$ and so $$c_{j+1} = \frac{c_j}{\psi \sqrt{j+1}} \left(\frac U2 j(j-1)-\mu j + \psi^2 -E_0 \right) - \sqrt{\frac{j}{j+1}}c_{j-1}.$$ At this point, I'm unsure of how to proceed as the ground state expansion coefficients are given by this recurrence relation which is quite complicated. How then do I determine the ground state energy as a function of $\psi$?