How is Green function in many-body theory introduced? Normally, for a (linear) operator $L$ and a DE 
$$ Lu(x) = f(x) $$
the Green function is defined as
$$ LG(x,s) = \delta(x-s) $$
and it is found that
$$ u(x) = \int G(x,s) f(s) ds $$
is the general solution of the DE.
Now, I've read some texts about Green functions in many-body theory (example), but the form is unfamiliar to me. 
Can you explain how those Green functions are introduced? I.e., why are objects of the form $\frac{1}{E - H_0 \pm i\eta}$ called Green functions (examples here and here)?
 A: Because these are actually Fourier transform of the usual Green functions.
Consider the Schrödinger equation :
$$
\hat{\mathcal{H}}|\Psi(t)\rangle=\mathrm{i}\partial_t|\Psi(t)\rangle
$$
The general solution $|\Psi(t)\rangle$ of such equation for a time-independant hamiltonian $\hat{\mathcal{H}}$ can be expressed in terms of Green function $G(x',x,t)$ :
$$
\Psi(x,t)=\langle x|\Psi(t)\rangle=\langle x|e^{-\mathrm{i}t\hat{\mathcal{H}}}|\Psi(t=0)\rangle=\int\mathrm{d}x'\,G(x',x,t)\,\Psi(x',t=0)
$$
where $G(x',x,t)=\langle x'|e^{-\mathrm{i}t\hat{\mathcal{H}}}|x\rangle$. The last equality is obtained by introducing the closure identity :
$$
\int\mathrm{d}x'\,|x'\rangle\langle x'|=\hat{1}
$$
One can then define a Green operator :
$$
\hat{G}(t)=-\mathrm{i}\,\Theta(t)\,e^{-\mathrm{i}t\hat{\mathcal{H}}}
$$
where $\Theta$ stands for the Heaviside step function which is here to ensure the causality of the solution $\Psi(x,t)$.
Then, one can compute the Fourier transform of such operator, which sometimes is called resolvent operator :
$$
\hat{G}(\epsilon)=\int\mathrm{d}t\,\hat{G}(t)\,e^{\mathrm{i}\epsilon t}=-\mathrm{i}\int\mathrm{d}t\,\Theta(t)\,e^{\mathrm{i}t(\epsilon-\hat{\mathcal{H}})}
$$
Then one can express the $\Theta$ function in terms of its Fourier transform :
$$
\Theta(t)=\int\frac{\mathrm{d}\omega}{2\pi\mathrm{i}}\,\frac{e^{-i\omega t}}{\omega-\mathrm{i}\eta}
$$
where $\eta$ is a positive infinitesimal parameter.
Taking all of this together, you will find that :
$$
\hat{G}(\epsilon)=-\frac{1}{2\pi}\int\mathrm{d}\omega\,\frac{1}{\omega-\mathrm{i}\eta}\int\mathrm{d}t\,e^{\mathrm{i}t(\epsilon-\omega-\hat{\mathcal{H}})}
$$
It is possible to recognize with the Fourier transform of the Dirac distribution :
$$
\delta(\omega)=\frac{1}{2\pi}\int\mathrm{d}t\,e^{\mathrm{i}\omega t}
$$
that :
$$
\hat{G}(\epsilon)=-\int\mathrm{d}\omega\,\frac{1}{\omega-\mathrm{i}\eta}\;\delta(\omega+\epsilon-\hat{\mathcal{H}})=\frac{1}{\epsilon+\mathrm{i}\eta-\hat{\mathcal{H}}}
$$
which is what you are looking for.
