Let $x\ge 0 $ and $y \ge 0 $ and consider a following improper integral.
\begin{equation} I(x,y) := \frac{1}{2\pi} \int\limits_{{\mathbb R}} \frac{(-1-\imath k y)}{k^2} e^{\imath k (x-y)} dk \end{equation}
We calculate this integral in two different ways.
First way:
\begin{eqnarray} I &=& \frac{1}{2\pi} \int\limits_{{\mathbb R}} e^{\imath k x} \cdot \int\limits_{{\mathbb R}} z e^{-\imath k z} 1_{z \ge y} dz \cdot dk \\ &=&\int\limits_y^\infty z \cdot \frac{1}{2\pi} \int\limits_{{\mathbb R}} e^{\imath k (x-z)} dk \cdot dz \\ &=& \int\limits_y^\infty z \delta(z-x) dz \\ &=& x \cdot 1_{x > y} \end{eqnarray}
Second way: We split the integral in question into two parts.
\begin{eqnarray} I_1 &=& \int\limits_{{\mathbb R}} \frac{e^{\imath k (x-y)}}{(\imath k)^2} = |x-y| \cdot \pi \\ I_2 &=& (-\imath y) \int\limits_{{\mathbb R}} \frac{1}{k} e^{\imath k (x-y)} dk = y sign(x-y) \cdot \pi \\ &&\Longrightarrow \\ I &=& I_1 + I_2 = \frac{1}{2} \cdot \left\{ \begin{array}{lll} x & if \quad x>y \\ -x & if \quad x< y \end{array} \right. \end{eqnarray}
The results are different so there must be an error somewhere. Now, I know that the second way is correct because I tested it numerically using the following code:
M = 100; eps = 10^(-4); {x, y} =
RandomReal[{-1, 1}, 2, WorkingPrecision -> 50];
NIntegrate[1/k^2 Exp[I k (x - y)], {k, -M, -eps}] +
NIntegrate[1/k^2 Exp[I k (x - y)], {k, eps, M}] -
2 Cos[eps (x - y)]/eps
-Abs[x - y] Pi
NIntegrate[ (-I y)/k Exp[I k (x - y)], {k, -M, -eps}] +
NIntegrate[ (-I y)/k Exp[I k (x - y)], {k, eps, M}]
Pi Sign[x - y] y // N
NIntegrate[1/k^2 (-1 - I k y) Exp[I k (x - y)], {k, -M, -eps}] +
NIntegrate[1/k^2 (-1 - I k y) Exp[I k (x - y)], {k, eps, M}] +
2 Cos[eps (x - y)]/eps
Pi (Abs[x - y] + Sign[x - y] y)
This means that the first way is wrong. Indeed in the first way we expressed the result through a double integral and then swapped the order of integration. The assumptions of the Fubini's theorem are not satisfied because the double integral from the modulus of the integrand does not exist.
Update: I have checked numerically that in the first approach the step from the definition to the first line, i.e. the expression as the double integral over $z$ and $k$, is wrong. Other steps seem to be right.
In[1222]:= M = 20; eps = 10^(-4); {x, y} =
RandomReal[{0, 1}, 2, WorkingPrecision -> 50];
If[x < y, tmp = x; x = y; y = tmp;];
(NIntegrate[1/k^2 (-1 - I k y) Exp[I k (x - y)], {k, -M, -eps}] +
NIntegrate[1/k^2 (-1 - I k y) Exp[I k (x - y)], {k, eps, M}] +
2 Cos[eps (x - y)]/eps)/(2 Pi)
NIntegrate[ Exp[I k x] z Exp[-I k z], {k, -M, M}, {z, y, M}]/(2 Pi)
NIntegrate[z Sin[M (x - z)]/(Pi (x - z)) , {z, y, M}]
x If[x > y, 1, 0]
Out[1224]= 0.255426 - 1.41358*10^-16 I
Out[1225]= 0.481756 - 4.85956*10^-15 I
Out[1226]= 0.481756
Out[1227]= 0.46787318300471700236343631477297526124617114997482
This is strange because that step follows very easily from rotating the integration contour in the complex $z$-plane and from elementary integration.
Having said all that my question is how can we salvage the first approach in order to get the right result? The use of Dirac's delta function is very handy and this function is well established in mathematics for a long time. So again, is it possible to get the correct result using the first approach?