Off the top of my head, a derivation might go like, as \delta goes to zero, for x = \omega + \mu - H = 0 non zero, the denominator stays non-zero and the numerator goes to zero, so you have a zero. For x = 0, then it goes like 1/\delta as \delta goes to zero, so this is divergent. The only trick from there would be to show that the integral of \delta/(x^2 + \delta^2) over all x is a constant value (I think it might be pi) independent of \delta.
My experience with Green's functions is fairly limited. So the only ways I know of to evaluate Green's functions require knowledge of the full spectrum of the Hamiltonian. But I will list the advantages I can think of, all of which are condensed matter topics.
The density of states you calculate from a non-interacting (or mean-field) theory doesn't correspond very well with experimental photo absorption/emission, because these experiments involve exciting a particle to a state with a finite lifetime that the whole system can have a response to. If the physical system were non-interacting it makes sense that the density of states would be an adequate description, because the excited particle wouldn't affect any of the other particles. But in a system with strong interactions, this isn't the case. The Green's function calculated via G(t) = \langle c(t) c^\dagger \rangle is a direct calculation of particle excitation, so from this argument, one would expect the Green's function to give a more physically realistic picture.
If you calculate the band structure (ie. eigenvalues of H_k) via LDA, you will have a set of energy bands which are sharp delta-function peaks at a specific energy values, all having the weight of one (ignoring degeneracy). This doesn't reflect experiment well because a variety of factors will broaden these peaks, or shift their weight to other bands. Some of these factors are: impurities, disorder, temperature, correlation effects. (In many cases these can be safely neglected, but there are certainly interesting cases where they are relevant). Impurities can be treated with the Coherent Potential Approximation, which uses a Green's function formalism (and that is the extent of my knowledge on CPA).
My knowledge of Green's functions comes from doing Dynamical Mean Field Theory (DMFT) calculations. In the DMFT formalism, we take the Green's function at each k-value:
G_{k}(\omega) = \frac{1}{\omega + \mu - H_k - \Sigma(\omega)}
where you can see a frequency dependent potential has been added to the Hamiltonian (the so-called self-energy). The spectrum of these Green's functions gives the usual k-space band structure, if the self-energy is zero. The self-energy is calculated via some other means, often with some self-consistency condition with the Green's function. The self-energy is non-Hermitian and may have a substantial imaginary part at some frequency. This can broaden the spectrum of H, create additional peaks, shift weight around, etc. to create the kinds of effects that appear in experiment. In DMFT the self-energy is used to include correlation effects. I'd imagine it's difficult / impossible to include any frequency dependent potential without using Green's functions.
