When we are introduced to the density of states in typical band-theory problems we neglect interaction between electrons, and thus we define the density of states of a sigle particle as: $D(E)=2\int_{1st BZ}\delta (E-\epsilon_\mathbf{k})d\mathbf{k}$ since the possible states available to occupy for an electron lie within a band (I assume here there is just one band described by $\epsilon_\mathbf{k}$).
Now, when we switch on the interactions the concept of single-electron density of states seems ill-defined to me, since there is not such a thing as "possible energy of an electron". Take the 1D Hubbard model Hamiltonian within mean-field approximation:
$H=\sum_{k}(\epsilon_{k\uparrow}n_{k\uparrow}+ \epsilon_{k\downarrow}n_{k\downarrow})-U N\langle n_\uparrow \rangle \langle n_\downarrow \rangle$
where $\epsilon_{k\sigma}=-2t\cos k+\langle n_{-\sigma}\rangle U$. $U$ is the on-site Coulomb repulsion and $t$ the hopping term and $\sigma$ the spin (=$\pm$).
In this case the eigenstates are multiparticle states, and thus we can ask only for the energy of the system, or average of one particle. How can I then calculate the DOS in this context? How do I define my bands in first place if I don't know what's the dispersion relation?