Integral involving two energy Green's functions The problem
I am attempting to evaluate this integral:
\begin{equation}
I(\vec{k}) = \lim_{\epsilon\to 0} \int d^3q \, 
\frac{1}{E-E_{\vec{q}}+i\epsilon}
\frac{1}{E-E_{\vec{k}+\vec{q}}+i\epsilon}
\end{equation}
with $E_\vec{q}=q^2/2$ and $E_{\vec{k}+\vec{q}}=|\vec{k}+\vec{q}|^2/2$.
If possible, it would be nice to see multiple methods of solution.
My attempt
I chose to work in a spherical coordinate system with $\vec{k}$ aligned along the $z$-axis. Letting $\phi$ represent the polar angle between $\vec{k}$ and $\vec{q}$, I wrote $E_{\vec{k}+\vec{q}}$ as 
\begin{equation}
E_{\vec{k}+\vec{q}}
= \frac{1}{2} \left( k^2 + q^2 - 2kq \cos (\pi-\phi) \right)
= \frac{k^2}{2} + \frac{q^2}{2} + kq \cos \phi .
\end{equation}
It is then fairly straightforward to evaluate both angular integrals, leaving only the radial integral:
\begin{equation}
I(\vec{k}) =
\lim_{\epsilon\to 0} \frac{2\pi}{k}
\int_0^\infty dq \, 
\frac{q}{E-q^2/2+i\epsilon} 
\ln \left| \frac{E-(k-q)^2/2+i\epsilon}{E-(k+q)^2/2+i\epsilon}\right| .
\end{equation}
But this is where I get stuck.
 A: As suggested by Philip Cherian, you may want to check in the Math section. However, I want to point out one thing about how to handle the limit.
I will assume that $|\vec{k}|\ne 0$. The first thing you want to do it separate each propagator into its real and imaginary parts. In general, $$\lim_{\delta \rightarrow 0^+}\frac{1}{E-E(q)+i\delta} = \lim_{\delta \rightarrow 0^+}\left[\frac{1}{E-E(q)}-i\frac{\delta}{[E-E(q)]^2 + \delta^2}\right],$$
where I multiplied top and bottom by the complex conjugate of the denominator, and disregarded the $\delta^2$ in the denominator of the real part, since it'll just go away when you take the limit. Then you use the identity $$\lim_{\delta\rightarrow 0^+}\frac{\delta}{[E-E(q)]^2+\delta^2} = \pi\delta(E-E(q)),$$
that is, the Dirac delta. What you get is $$\begin{split}&\int d^3q\,\left[\frac{1}{E-E(q)}-i\pi\delta(E-E(q)) \right]\left[\frac{1}{E-E(|\vec{q}+\vec{k}|)}-i\pi\delta(E-E(\vec{q}+\vec{k})) \right]\\ &=\int d^3q\,\frac{1}{E-E(q)}\frac{1}{E-E(|\vec{q}+\vec{k}|)} + i\pi\int d^3q\,\left[ \frac{\delta(E-E(q))}{E-E(|\vec{q}+\vec{k}|)} + \frac{\delta(E-E(|\vec{q}+\vec{k}|))}{E-E(q)} \right]. \end{split}$$
Notice that, since I'm assuming $k\ne 0$, the product of two Dirac deltas will have to vanish, because their arguments are never simultaneously zero. You can easily evaluate the integral with the Dirac deltas. Recall that, since the dispersion is quadratic in $q$, you need to follow certain rules when you change variables. That gives you the imaginary part of the integral. I have no real advice on how to evaluate the real part $$\mathrm{Re}\{I(\vec{k})\}=\int d^3q\,\frac{1}{E-E(q)}\frac{1}{E-E(|\vec{q}+\vec{k}|)},$$ except to check that you actually need it. Depending on the problem you're solving, sometimes all you need is the imaginary part. On the other hand, I also recommend you check that both denominators in $I(\vec{k})$ have $+i \epsilon$, instead of one having $+i \epsilon$ and the other $-i\epsilon$. I say this because that expression looks like the skeleton approximation to a polarization bubble, which is the product of a retarded and an advanced propagator, in which case their imaginary parts have opposite signs. 
A: This is a modification of classic manipulations that are done when computing loop integrals. The tool to use is usually referred to as "Feynman-parameters''. Inserting in the definitions of $ E _i $ and defining $ m  ^2 \equiv 2 E $ to make connection with the usual computations I find (I drop the $\epsilon$'s since they are irrelevant here):
\begin{align}
  \frac{1}{2} I & =  \int d^3q \frac{1}{ ( q ^2 - m ^2  ) ( (  {\mathbf{q}} + {\mathbf{k}}  ) ^2 - m ^2 ) } \\
 & = \int _0 ^1 dx \int d^3q \frac{1}{ \left[ ( ( {\mathbf{q}} + {\mathbf{k}} )  ^2 - m ^2  ) x + ( q ^2 - m ^2 ) ( 1 - x ) \right] ^2 }
\end{align}
Expanding and rewriting the expression a bit gives,
\begin{equation}
\frac{1}{2} I = \int _0 ^1 d x  \int d^3q \frac{1}{ \left[ ( {\mathbf{q}} + {\mathbf{k}} x ) ^2 - {\mathbf{k}} ^2 x ^2 + {\mathbf{k}} ^2 x - m ^2  \right] ^2 } 
  \end{equation} 
  Now shifting the integral variable I find:
  \begin{equation}
\frac{1}{2} I = \int _0 ^1 d x  \int d^3q \frac{1}{ \left[ {\mathbf{q}} ^2  + \Delta   \right] ^2 } 
  \end{equation} 
  where $ \Delta \equiv - {\mathbf{k}} ^2 x ^2 + {\mathbf{k}} ^2 x - m ^2 $. The angular part of the integral is now trivial and the radial integral is straightforward. I find,
  \begin{equation}
\frac{1}{2} I = \pi ^2 \int _0 ^1 d x \Delta ^{ - 1/2} 
\end{equation}
The $ x $ integral is usually left untouched, but I suspect in this case you could try carrying that out too.
