Density of states in a system of interacting electrons

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?

$$\mu$$ would be a proper generalization of a single partile's engery in many body systems. $$G=G(N)$$ $$\mu =\frac{\partial G}{\partial N} =G(N+1)-G(N)$$ $$\text{Density of states:} \quad \rho(\mu) := \frac{\partial N}{\partial \mu}$$