Is the Feynman propagator integrable over Minkowski spacetime? Setup: Consider the Feynman propagator as it appears in Peskin and Schroeder eq. 2.59:
\begin{align}
D_F(x-y)\equiv i\int_{\mathbb{R}^4} \frac{d^4p}{(2\pi)^4}\frac{e^{ip\cdot(x-y)}}{p^2-m^2+i\epsilon},
\end{align}
where $p,x,y$ denote 4-vectors; $p^2\equiv \eta_{\mu\nu}p^\mu p^\nu$, $p\cdot (x-y)\equiv \eta_{\mu\nu}p^\mu (x-y)^\nu $;  the sign convention here is $\text{diag}(\eta)=(+1,-1,-1,-1)$; and the $i\epsilon$ has been inserted for convergence.
The following "formal" manipulations would seem to imply $i\int d^4x D_F(x-y)=1/m^2$.    
\begin{align}
i\int_{\mathbb{R}^4}d^4x D_F(x-y)&=-\int_{\mathbb{R}^4} d^4x \int_{\mathbb{R}^4}\frac{d^4p}{(2\pi)^4}\frac{e^{ip\cdot(x-y)}}{p^2-m^2+i\epsilon}\\
&=- \int_{\mathbb{R}^4}d^4p\frac{e^{-ip\cdot y}}{p^2-m^2+i\epsilon}\int_{\mathbb{R}^4} \frac{d^4x}{(2\pi)^4}e^{ip\cdot x}\\
&=- \int_{\mathbb{R}^4}d^4p\frac{e^{-ip\cdot y}}{p^2-m^2+i\epsilon} \delta^{(4)}(p)\\
&=\frac{1}{m^2}, \tag{*}\label{*}
\end{align}
where I used a standard $\delta$-identity in going from the 2nd to the 3rd line, and took the $\epsilon\to 0$ limit after performing the $p$-integral.
Question(s) These manipulations are "formal" in the sense that I have interchanged the order of integration without justification.  I suppose, in general, one would want to apply something akin to Fubini's theorem.  However, I am not sure to what extent the regularization implicit in the Feynman propagator might alter such an argument.  
Can one justify these formal manipulations?  To what extent is $\eqref{*}$ a true result and is there any physical intuition for it?
 A: All of this becomes mathematically rigorous within the theory of temperate Schwartz distributions and their Fourier transforms. Put $y=0$ for simplicity.
Then $D_{F}(x)$ is in $S'(\mathbb{R}^4)$. Now you are taking its Fourier transform. No problem. Then you are evaluating it at zero momentum. That needs care because the Fourier transform (which is also the propagator in momentum space) is a distribution too. Fortunately, zero is not in the singular support of this distribution, so pointwise evaluation is legit and gives the $1/m^2$. Your computation simply amounts to the Fourier inversion theorem for a temperate distribution.
To learn about distributions see the old book by Laurent Schwartz "Mathematics for the Physical Sciences", or the book by Strichartz "A Guide to Distribution Theory and Fourier Transforms".
A lower tech approach (but with the theory of distributions lurking in the background)
is to do Fubini etc. with the propagator
$$
\frac{e^{-\epsilon_1 p^2}}{p^2-m^2+i\epsilon_2}
$$
and then take the limits $\epsilon_1,\epsilon_2\rightarrow 0^{+}$.
Careful: the $p^2$ in the exponential is with Euclidean signature.
This brings you from $S'(\mathbb{R}^4)$ to $S(\mathbb{R}^4)$ and basically proving the Fourier inversion formula in the latter space. In fact this proof usually requires introducing another $e^{-\epsilon_3 x^2}$, again with Euclidean signature. Good luck with all the epsilons.
PS: I said the theory of distributions is lurking in the background because of the following two facts. The space $S(\mathbb{R}^d)$ of rapidly decaying functions (for which classical tools like Fubini etc. work well) is dense in the space of distributions $S'(\mathbb{R}^d)$. The second fact is that operations on distributions like taking the Fourier transform are defined as the unique continuous extensions of their classical versions on $S(\mathbb{R}^d)$. So putting all the epsilons amounts to doing distributional analysis classically, i.e., at the level of a first graduate course in analysis, measure theory etc.
A: Only formally, in the same sense that
$$\int_{-\infty}^\infty \operatorname{d}x\, \frac{\operatorname{e}^{ikx}}{2\pi} = \delta(k).$$
Why? You can check that if you integrate the Feynman propagator in position space over all space you'll get
$$\int \operatorname{d}^3y\, D_F(x-y) = \frac{1}{2m}\sin(m|x^0-y^0|) - \frac{i}{m} \cos(m[x^0-y^0]), $$
where I may have screwed up the scale in front of the cosine term, but the one in front of the sine term is right. Getting that sine term right is an essential part of satisfying the boundary conditions in the defining equation for the propagator
$$\Box D_F(x-y) + m^2 D_F(x-y) = \delta^4(x-y).$$
Just integrate both sides of the defining equation over all space to see what I'm talking about.
Doing the time integral on the sine and cosine terms, you'll find that it converges in the mean to
\begin{align}
    \int \operatorname{d}x^0 \left[\frac{1}{2m}\sin(m|x^0-y^0|) - \frac{i}{m} \cos(m[x^0-y^0])\right] &= \frac{1}{m^2}\left[-\cos(mu) - 2i \sin(mu)\right]_{u=0}^{u_\mathrm{max}} \\
     & = \frac{1}{m^2}  - \frac{\cos(mu_{\mathrm{max}}) +2i\sin(mu_{\mathrm{max}})}{m^2} \\
     &\rightarrow \frac{1}{m^2}\ \mathrm{in\ the\ mean.}
\end{align}
