Fourier transform of the free propagator squared - $\int d^{4}p\ \frac{e^{-i p\cdot x}}{p^{2}+m^{2}-i\epsilon}$ The point of the question is to ask what is the function given by the following integral:
$$
H(x,y) \ \equiv \ \int \frac{d^{4}p}{(2\pi)^{4}} \frac{e^{-i p \cdot (x-y)}}{(p^{2}+m^{2}-i\epsilon)^{2}}
$$
This is closely related to the propagator (for $(x-y)^{2} < 0$):
$$
G(x,y) \ = \ \int \frac{d^{4}p}{(2\pi)^{4}} \frac{e^{-i p \cdot (x-y)}}{p^{2}+m^{2}-i\epsilon} \ = \ - \frac{ i m }{ 4 \pi^{2} \sqrt{ - (x-y)^{2} }} K_{1}\left( m \sqrt{ - (x-y)^{2} } \right)
$$
The reason I ask this question, is because in calculations like for the following Feynman diagram:

I would have in position space the following:
$$
\propto \int d^{4}u\ G_{0}(x,u) G_{0}(u,u) G_{0}(u,y)\\
= \int d^{4}u \left[ \int \frac{d^{4}p}{(2\pi)^{4}} \frac{e^{-i p \cdot (x-u)}}{p^{2}+m^{2}-i\epsilon} \right] \left[ \int \frac{d^{4}k}{(2\pi)^{4}} \frac{1}{k^{2}+m^{2}-i\epsilon} \right] \left[ \int \frac{d^{4}q}{(2\pi)^{4}} \frac{e^{-i p \cdot (u-y)}}{q^{2}+m^{2}-i\epsilon} \right] \\
= \int \frac{d^{4}p}{(2\pi)^{4}} \frac{e^{-i p \cdot (x-y)}}{(p^{2}+m^{2}-i\epsilon)^{2}} \ \int \frac{d^{4}k}{(2\pi)^{4}} \frac{1}{k^{2}+m^{2}-i\epsilon} \\
= H(x,y) \int \frac{d^{4}k}{(2\pi)^{4}} \frac{1}{k^{2}+m^{2}-i\epsilon}
$$
You can regulate the integral over $k$ however you like, but what to do with $H(x,y)$ here?
It seems that $H(x,y)$ pops up a lot when you do calculations like this. Is there a way to evaluate $H(x,y)$? Surely this must have been done somewhere?
 A: Note that you don't need to go through the Fourier transform to get the result you're seeking. All you need is translation invariance of $G$ (i.e. $G(x,y) = G(x+w, y+w) = G(|x-y|)$) to say that the first integral gives:
\begin{align}
I & \equiv \int \operatorname{d}^4 u \ G_0(x,u)\, G_0(u,u)\, G_0(u,y) \\
& = G(0) \int \operatorname{d}^4 u \ G_0(x,u)\, G_0(u,y),
\end{align}
as long as the Green's function $G$ has been suitably regularized to make $G(0)$ finite. Note that the convolution chain can grow infinitely long. Consider an exactly solvable perturbation of the form $\Delta \mathcal{L} = - \frac{1}{2}(\Delta m^2) \phi^2$, you'll get just such an infinite chain in producing the net propagator with shifted mass.
One way to calculate $H$ is to note that it can be produced by the limit:
\begin{align}
H(x,y) &= \lim_{\mu^2 \rightarrow 0} \int \operatorname{d}^4 u \ G_0(x,u; m)\, G_0(u,y; \sqrt{m^2 + \mu^2}) \\
& = \lim_{\mu^2 \rightarrow 0} \int \frac{\operatorname{d}^4 p}{(2\pi)^4} \frac{\mathrm{e}^{-ip\cdot(x-y)}}{ \left(p^2 + m^2 - i\epsilon \right) \left(p^2 + m^2 + \mu^2 - i\epsilon \right) } \\
& = \lim_{\mu^2 \rightarrow 0} \int \frac{\operatorname{d}^4 p}{(2\pi)^4} \frac{\mathrm{e}^{-ip\cdot(x-y)}}{\mu^2} \left[\frac{1}{p^2 + m^2 - i\epsilon} - \frac{1}{p^2 + m^2 + \mu^2 - i\epsilon}\right] \\
& = \lim_{\mu^2 \rightarrow 0} \frac{1}{\mu^2} \left[G_0(m|x-y|) - G_0\left(\sqrt{m^2 + \mu^2}|x-y|\right)\right],
\end{align}
where I used partial fraction decomposition in the intermediate step. Working out the final result, which can likely be done using the Bessel function derivative identities, is left as an exercise for the reader.
