How to get an imaginary self energy? The Lehman representation of the frequency-dependent single particle Green's function is
$$G(k,\omega) = \sum_n \frac{|c_k|^2}{\omega - E_n + i\eta}$$
where $n$ enumerates all the eigenstates of the system, $c_k$ is the overlap between 
$|\Psi_n\rangle$, the eigenstate to eigenenergy $E_n$, and the state $|k\rangle = 
c_k^\dagger | 0 \rangle$. And $\eta$ is a very small convergence factor, sometimes also written as $0^+$.
The corresponding spectral function is
$$A(k,\omega) = \sum_n |c_l|^2 \frac{\eta}{\pi} \cdot \frac{1}{(\omega - E_k)^2 + \eta^2)}$$
This means we have Lorenzian with width $\eta$ at the eigenergies of the system, i.e., for small $\eta$
we see that the spectral function has sharp peaks for all values of $\omega$ that correspond to the system's eigenergies. (Note that for $\eta \rightarrow 0$, the Lorenzian becomes a $\delta$-function).
So far, so good. Now we can also write the Green's function using the self energy. If $\varepsilon(k)$ is the dispersion of the non-interacting particle, then the self energy is defined via
$$G(k,\omega) = \frac{1}{\omega - \varepsilon(k) - \Sigma(k,\omega) + i\eta}$$
It is mentioned on many occasions that a self-energy with a non-zero imaginary part corresponds to a quasi-particle with a finite lifetime. In the spectral function, an imaginary self energy leads to a broad peak, and the width of that peak is related to the imaginary part of $\Sigma(k,\omega)$.
Question: At first, I have shown that the spectral function has infinitely sharp peaks at precisely the system's eigenergies. But later, I have shown that an imaginary self energy leads to broad peaks. How can I reconcile these two facts?
Related Question: If I find a "continuum" of states in my spectral function, how can I tell if these form a continuum of "true" eigenstates of the system or a bunch of decaying quasiparticles?
The motivation: In the Holstein model, we have a single electron on a lattice in the tight-binding approximation interacting with Einstein phonons of frequency $\Omega$. In the non-interacting limit, we have the non-interacting electron band $\varepsilon(k) = -2t \cos(ka)$ and then, beginning at an energy $\Omega$ above $-2t$, a continuum of a states corresponding to an electron of momentum $k - q$ and a phonon of momentum $q$. These are "true" eigenstates of the (non-interacting) system. But if I then turn on the interaction of the type
$$ \frac{g}{\sqrt{N}} \sum_q c_{k-q}^\dagger c_k (b_q^\dagger + b_{-q})$$
I still have a continuum an energy $\Omega$ above the ground state, but now these are not true eigenstates any more, since the additional phonon can decay via the interaction term. But when they are not true eigenstates any more, then they shouldn't show up in the Lehman representation and hence shouldn't contribute to the spectral function?
 A: Unstable particles are concepts of effective field theories (or few-particle systems) in reduced descriptions where the decay products are ignored. In these reduced descriptions, they appear as particles with complex masses, and their Green's functions have complex poles.
In an unreduced description, unstable particles appear as poles of the analytically continued resolvent to the second, unphysical sheet. The spectral decomposition of the resolvent is a sum over the discrete spectrum plus an integral over the continuous spectrum along a complex curve defined by the $+i\epsilon$ prescription. By embedding the Hilbert space in a bigger space, the so-called rigged Hilbert space (or Gelfand triple), one can deform the integration path to a complex detour, as long as no pole is crossed. 
If one crosses a pole (which can happen only by deforming the path into the nonphysical sheet), an additional discrete term proportional to $1/(E-E_0)$ appears, which is the contribution  of the corresponding resonance. If the imaginary part of $E_0$ is tiny, corresponding to a long-living particle, and the residue (the numerator of the appearing term) is not tiny, then undoing the deformation implies a huge spike in the spectral distribution of the continuous spectrum, nearly a delta fucntion, which is why these objects can be approximately treated as particles. (This can be seen in full detail for the so-called Wigner--Weisskopf atom.) If the imaginary part is larger, one gets less pronounced peaks (resonances), which is the way shortlived unstable particles are detected in practice. 
This is well explained in the book
Kukulin, V. I. ; Krasnopol'sky, V. M. ; Horáček, J.,
Theory of resonances : principles and applications,
Kluwer 1989
together with the computation aspects for few-particle systems.
If the resolvent can be analytically continued everywhere near the real continuous spectrum, one can move the integration contour into the unphysical sheet until all significant resonances have become part of the now complex discrete spectrum (of the Hamiltonian in an indefinite non-Hilbert space obtained by a corresponding deformation of the integral in the inner product). The residual continuous spectrum has no significant peaks left, and can usually be neglected, resulting in a much more tractable effective subsystem.
To explain why the self-energy is typically complex, solve the defining equation $G(E)=\frac{1}{E-\epsilon(E)-\Sigma(E)}$ (valid for nonreal $E$ only) for
$\Sigma(E)=E-\epsilon(E)-G(E)^{-1}$. Note that this expression is not Hermitian, as under conjugation, $E$ becomes $\overline E \ne E$! Thus, even in the limit where $E$ becomes real (which is now often possible because poles of $G(E)$ don't create harm), there is no reason for $\Sigma(E)$ to be Hermitian, as $G(E)$ looks different on both sides of the branch cut given by the real continuous spectrum.
A: This is a great question. The point is we've made a huge structural change in going from Lehmann representation to the, um, other one. In Lehmann representation, we've chosen to write $G$ as a sum over an infinite number of real poles, which is fine of course, but when you start adding infinite numbers of delta functions together, things can get tricky.
(Aside: consider for a moment if we integrated Im$G$ over all $k$ -- we're left with the density of states $\rho(\omega)=\sum_n \delta(\omega-E_n)$... and yet when you sit down and work it out for, say, free particles in continuous 3-space, you get a nice smooth function of $\omega$. Why? Well because the eigenvalues $E_n$ become arbitrarily close to one another, the sum becomes an integral, under which delta functions immediately switch from foe to friend.)
So with that in mind, the spectral function you got from the Lehmann representation will have structure at all of the eigenvalues of the full many-body system (which will merge into smooth functions in an infinite system), but this is silly since those eigenvalues don't generally correspond to excitations that we could ever think about in any intuitive way.
So, instead of ever thinking about the real eigenstates of many-body systems, we give up and introduce quasiparticles and work with their Green functions, which, as you showed, have one complex pole related to $\Sigma$, with the imaginary part giving the spread on the real-frequency axis that the real eigenvalues have.
(I think this should help with your related question too, but it's late -- let me know if I've been opaque.)
