Massless $m=0$ 4D Fourier transform of $(p^2 + i \epsilon)^{-2}$ This question is related to this one. I'm assuming that we're in or on the the light-cone $s \leq 0$ in what follows.
Suppose I'm interested in computing the following Fourier transform, in the massless $m=0$ case:
$$
H(x,y;m) = \int \frac{d^{4}p}{(2\pi)^4} \frac{e^{-i p \cdot (x-y)}}{(p^2 + m^2 + i \epsilon)^2}
$$
The $\epsilon$ is there in the Feynman-prescription sense. We know that the massive $m \neq 0$ propagator is given by;
$$
G_0(x,y;m) \ = \ \int \frac{d^4p}{(2\pi)^4} \frac{e^{-i p \cdot (x-y)}}{p^2 - m^2 + i \epsilon} \ = \ - \frac{i}{4 \pi^2} \frac{m}{\sqrt{-s}} K_{1}\left( m\sqrt{-s} \right) - \frac{1}{4 \pi} \delta(s)
$$
Where $s = (x^0 - y^0)^2 - (\mathbf{x} - \mathbf{y})^2$. Setting $\mu = m^2$, we notice that;
$$
\frac{\partial G_0(x,y;\sqrt{\mu})}{\partial \mu} = \frac{\partial }{\partial \mu } \left( \int  \frac{d^4p}{(2\pi)^4} \frac{e^{-i p \cdot (x-y)}}{p^2 - \mu + i \epsilon} \right) = \int \frac{d^{4}p}{(2\pi)^4} \frac{e^{-i p \cdot (x-y)}}{(p^2 - \mu + i \epsilon)^2} = H(x,y;\sqrt{\mu})
$$
So differentiating in the above manner, I have been able to find:
$$
H(x,y;m) = \frac{i}{8 \pi^2} K_{0}\left( m \sqrt{-s} \right)
$$
This seems like a good answer, however things go wrong when I set $m=0$. I find that:
$$
H(x,y;0) \ = \ \lim_{m\to 0^+} \left[ \frac{i}{8 \pi^2} \log\left( \frac{m(-s)}{4} \right) \right]
$$
Meaning the Fourier transformed function diverges like a log as I take the mass to zero. Why is this happening! Surely we can examine massless scalar theories? How can a rectify the above?
 A: The problem is that the Green's function is not analytic on the light cone, so the order of limits matters. 
One approach would be to go back to the differential equation that this chain is meant to solve. You can show that if
$$\left[\frac{\partial^2}{\partial t^2} - \nabla^2\right] G(t,\mathbf{x};t',\mathbf{x}') = \delta(t-t')\,\delta(\mathbf{x}-\mathbf{x}')$$
and we want to find a function that satisfies
$$\left[\frac{\partial^2}{\partial t^2} - \nabla^2\right]^2 G_2(t,\mathbf{x};t',\mathbf{x}') = \delta(t-t')\,\delta(\mathbf{x}-\mathbf{x}')$$
then the following works
$$G_2(t,\mathbf{x};t',\mathbf{x}') = \int \operatorname{d}^4 x'' G(t,\mathbf{x};t'',\mathbf{x}'')\, G(t'',\mathbf{x}'';t',\mathbf{x}').$$
Though where to take this "back to the beginning" approach from here to solve the problem, I'm not sure. One idea is to Fourier transform the spatial coordinates, giving you a 4th order ODE you can solve per mode, though it isn't clear what boundary conditions you'd need to apply to that to get the Feynman propagator.
I think your best bet, then, would be to take what I'd describe as "the standard approach." You just take your original expression for $H(x,y;m)$ with $m=0$ and evaluate the $p^0$ integral using the residue theorem from complex analysis. The remaining integrals you can probably find on a table of Fourier transforms.
