Multiple Fourier Integrals involving Heaviside Theta Function I want to evaluate the integral:
$$I=\int_{-\infty}^{\infty}dx_1 \int_{-\infty}^{\infty}dx_2 \  \Theta(x_1-x_2) \ e^{i(ax_1+bx_2)}$$
where $\Theta(x)$ is the Heaviside function.
What I was doing now was taking the relation for $\Theta$: $\Theta (x)=-\frac{1}{2\pi i}\int_{-\infty}^{\infty}d\tau \frac{1}{\tau + i\epsilon} e^{-ix\tau}$ and I got:
$$I= -\frac{1}{2\pi i}\int_{-\infty}^{\infty}dx_1 \int_{-\infty}^{\infty}dx_2\int_{-\infty}^{\infty}d\tau \ \frac{1}{\tau + i\epsilon}e^{-i(\tau-a)x_1}e^{-i(\tau+b)x_2}\\=2\pi i\int_{-\infty}^{\infty}d\tau\ \frac{1}{\tau + i\epsilon}\delta(\tau-a)\delta(\tau+b) =2\pi i\frac{\delta(a+b)}{a + i\epsilon}$$
I didn't know if the integral was convergent and I could simply interchange the integrals, so I tried it in a different form with $X=x_1+x_2$ and $x=x_1-x_2$ :
$$I=\frac{1}{2}\int_{-\infty}^{\infty}dX \int_{-\infty}^{\infty}dx \ \Theta(x) \ e^{ia\frac{x+X}{2}+ib\frac{X-x}{2}} \\ =\pi \int_{-\infty}^{\infty}dx \ \Theta(x) e^{-2\pi i\frac{b-a}{4\pi}}\delta(\frac{b+a}{2})$$
With the Fourier Transform of the Heaviside function $\int_{-\infty}^{\infty}dk\ \Theta(k)e^{-2\pi i kx}=\frac{1}{2}(\delta(x)-\frac{i}{\pi k}) $ I get
$$I=\pi \left(2\pi\delta(b-a)-\frac{4 i}{b-a}\right)\delta(a+b)=2\pi^2\delta(a)\delta(b)+2\pi i\frac{\delta(a+b)}{a}$$
I don't know yet where the $\delta(a)\delta(b)$ should come from in the first method. When I want to check that now and integrate $I$ over $a$ and $b$ I get from the first line:
$$\int_{-\infty}^{\infty}da \int_{-\infty}^{\infty}db \ I = \int_{-\infty}^{\infty}dx_1\int_{-\infty}^{\infty}dx_2 \Theta(x_1-x_2) \delta(x_1)\delta(x_2)  \\ = \Theta(0)=\frac{1}{2}$$ and from the second result:
$$\int_{-\infty}^{\infty}da \int_{-\infty}^{\infty}db \ I=4\pi^2-\int_{-\infty}^{\infty}db \frac{1}{b}=-\infty$$
Where did it go wrong? Is the integral correct?
EDIT: corrected mistake in derivation because of comment.
 A: Products of distributions is a well-known delicate topic, see e.g. this mathoverflow post. Unless if each distribution factor depend on different variables, such as, e.g. $\delta^3(\vec{r})=\delta(x)\delta(y)\delta(x)$. So the safest thing is to write OP's integral in two manifestly independent variables. E.g. 
$$ I(a_1,a_2) ~:=~ \iint_{\mathbb{R}^2} \! \mathrm{d}x^1 \mathrm{d}x^2~ 
e^{i (a_1 x^1+a_2 x^2)}  \theta(x^1-x^2) $$
$$\tag{1}~=~ \iint_{\mathbb{R}^2} \! \mathrm{d}y^1 \mathrm{d}y^2 ~ 
e^{i (b_1 y^1 + b_2 y^2)}  \theta(y^1)~=~I_1(b_1)~I_2(b_2) .$$
Here we have defined a linear coordinate transformation (with unit determinant to avoid a Jacobian factor for simplicity)
$$ \tag{2} \begin{bmatrix}y^1 \\ y^2 \end{bmatrix}
~=~\begin{bmatrix} \frac{1}{2} &-\frac{1}{2} \\ 1-t & t \end{bmatrix} \begin{bmatrix}x^1 \\ x^2 \end{bmatrix} $$
in such a way that the new variable $y^1$ becomes the argument of the Heaviside step function. Here $t\in\mathbb{R}$ is a free parameter. It is interesting to trace how the integral (1) stays independent of the value of this free parameter $t$. The inverse transformation is
$$ \tag{3} \begin{bmatrix}x^1 \\ x^2 \end{bmatrix}
~=~\begin{bmatrix}  t&\frac{1}{2} \\ t-1 & \frac{1}{2} \end{bmatrix} \begin{bmatrix}y^1 \\ y^2 \end{bmatrix}. $$
In order for the argument $a_1 x^1+a_2 x^2=b_1 y^1 + b_2 y^2$ of the exponential (1) to stay form invariant, we define new variables
$$ \tag{4} \begin{bmatrix}b_1 \\ b_2 \end{bmatrix}
~=~\begin{bmatrix}  t&  t-1 \\ \frac{1}{2}& \frac{1}{2} \end{bmatrix} \begin{bmatrix}a_1 \\ a_2 \end{bmatrix}. $$
The first integral factor is
$$\tag{5} I_1(b_1)~:= \int_{\mathbb{R}} \! \mathrm{d}y^1~e^{i b_1 y^1}\theta(y^1)
~=~\frac{i}{b_1+i0^+}~=~P\frac{i}{b_1}+\pi\delta(b_1),$$
cf. the Sokhotski–Plemelj formula. The second integral factor is
$$\tag{6} I_2(b_2)~:= \int_{\mathbb{R}} \! \mathrm{d}y^2~e^{i b_2 y^2}
~=~2\pi ~\delta(b_2).$$
While the two factors (5) and (6) clearly depend on $t$, their product
[=OP's integral (1)] is independent of $t$:
$$ I(a_1,a_2)~=~I_1(b_1)~I_2(b_2)~=~\frac{2\pi i}{a_1+i0^+}\delta(a_1+a_2)$$
$$\tag{7} ~=~2\pi i ~P\frac{1}{a_1}\delta(a_1+a_2)+2\pi^2 \delta(a_1)\delta(a_2). $$
