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)?
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.
