Help with understanding Green's Functions 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).
 A: In mathematical physics the resolvent of the Hamiltonian $H$, the operator
associated with the Green's function, is quite important in studying various
aspects of $H$. See, for instance, the various volumes of  "Methods of
modern mathematical physics " by Reed and Simon. Consider the time evolution
\begin{equation*}
\psi (t)=\exp [-iHt]\psi (0)
\end{equation*}
Its complex Laplace transform is
\begin{equation*}
\hat{\psi}(z)=\int_{0}^{\infty }dt\exp [izt]\exp [-iHt]\psi (0),\;Im
z>0
\end{equation*}
Performing the integration
\begin{equation*}
\hat{\psi}(z)=i[z-H]^{-1}\psi (0)
\end{equation*}
Note that here the "energy " is non-real so the $E+i\varepsilon $ occurs
naturally. In the simplest case the Green's function is
\begin{equation*}
G(x,y,z)=<x|[z-H]^{-1}|y>
\end{equation*}
and we can define
\begin{equation*}
G(x,y,E)=\lim_{\varepsilon \downarrow 0}G(x,y,E+i\varepsilon )
\end{equation*}
if it exists which is usually not the case. The advantage of introducing $
[z-H]^{-1}$ instead of $\exp [-iHt]$ is that it is much easier to study its
properties and the perturbation theory of the spectrum of $H$.
In conclusion the Green's function is a useful tool but apparently not in
the case you are interested in.
A: In http://arxiv.org/abs/1602.07793, I compute the honest (in the sense of inverting a differential operator) Green's function of the time-dependent Schrodinger equation, for a system of non-interacting particles. 
In short, the full Green's function is a simple extrapolation of the single-particle Green's function, and may be computed rather easily in the non-interacting case (though no one is interested in such an object, for some reason).
