I. The Green's Function Method
The Green's function is immensely useful as a tool in Solid State Physics. Using a Green's function, one can compute all relevant data from a physical system. For example, the Green's function for the time-independent Schrodinger equation (TISE), $$G(E):=\frac{1}{H-E}$$ yields both the density of states, $$-\frac{1}{\pi}\lim_{\epsilon\to 0^+}\text{Im}\,\text{Tr}\,G(E+i\epsilon)=\sum_n\delta(E-E_n)~~~~~~~~~~~~~~~$$ that is, the eigenvalues, and, letting $\{\psi_n\}$ denote the associated eigenstates, the Green's function also yields the projected density of states, $$-\frac{1}{\pi}\lim_{\epsilon\to 0^+}\text{Im}\,(f_0,G(E+i\epsilon)f_0)=\sum_n|(f_0,\psi_n)|^2\delta(E-E_n)$$ which is equivalent to the eigenstate data. Moreover, Green's functions allow us to efficiently formulate effective field theory, perturbation theory, and the renormalization group in the Hamiltonian picture, which is indispensable.
II. Many-Body Green's Functions: Where Everything Falls Apart:
However, this is all misleading: when physicists mention, "the Green's function of a non-interacting Hamiltonian $H$", which is, explicitly, $$(\psi_0, T\{\Psi^\dagger(x,t)\Psi(x',t')\}\psi_0),$$ Where $\psi_0$ is the groundstate of $H$, a literalist would think that they mean the Green's function for the many-body, time-dependent Schrodinger equation (TDSE): $$\frac{1}{\frac{i}{\hbar}H-\partial_t},~~~~~~ H=\sum_{ij}A_{ij}c^\dagger_ic_j. $$ However, close calculation shows that this is instead the Green's function of the associated single-particle Hamiltonian: $$(\psi_0, T\{\Psi^\dagger(x,t)\Psi(x',t')\}\psi_0)=\frac{1}{\frac{i}{\hbar}\mathcal H-i\partial_t},~~~~~~\mathcal H= \sum_{ij}A_{ij}\,\left|f_i\right>\left<f^j\right|.$$ However, this does not generalize straightforwardly to interacting systems. In particular, for an interacting system, there is no such single-particle hamiltonian! So the above virtues of the Green's function method no longer hold. We do not have the density of states, and we don't have the projected density of states.
So here's my question: what use is this method if it only characterizes non-interacting systems, which we already know how to solve? (Also, this gives a very boring renormalization group flow).