Unexpected distribution of ab-cd where a,b,c,d are independent and N(0,1) distributed I have discovered an unexpected curiosity that if
a,b,c,d are independent random variables and are N(0,1) distributed, then
|ab-cd| is exponentially distributed Exp(1). It's easy to verify that numerically.
Is there a proof of that ? or does it stem from some deeper fact ?
I noticed that it when I was considering cross products of random Gaussian distributed vectors and noticed that their components have \exp(-|x|/a) distribution.
Demonstration with python
N = 1000000
x = np.random.normal(size=(N,4))
z = np.abs(x[:,0] * x[:,1] - x[:,2] * x[:,3])
plt.hist(scipy.stats.expon.cdf(z), bins=100)


 A: Writing
$$ab - cd = \left(\left[\left(\frac{a+b}{\sqrt 2}\right)^2 + \left(\frac{c+d}{\sqrt 2}\right)^2 \right] - \left[\left(\frac{a-b}{\sqrt 2}\right)^2  + \left(\frac{c-d}{\sqrt 2}\right)^2 \right]\right)/2$$
and noting that $(a+b, c+d, a-b, c-d)/\sqrt{2}$ has a standard Normal distribution, it is immediate (from the definitions of chi-square distributions arising from sums of squares of independent standard Normal variables) that the distribution of $ab -cd$ is twice that of the difference of two independent chi-squared(2) variables.
Since chi-squared(2) variables have Exponential$(1/2)$ distributions, $ab-cd$ is distributed as the difference of two Exponential$(1)$ distributions.  Since the characteristic function of an Exponential$(1)$ distribution is $\phi(t) = 1/(1 - it),$ the cf. of  the difference is
$$\phi(t)\phi(-t) = \frac{1}{1+t^2}.$$
The Fourier Transform of that is proportional to $f(\omega)=\exp(-|\omega|),$ showing that the absolute value $|ab-cd|$ must have a density proportional to $f(x)$ for $x\ge 0:$ that's the Exponential$(1)$ distribution.
The thread at Probability function for difference between two i.i.d. Exponential r.v.s gives a direct demonstration of this latter result via integration of the convolution.
Reference
(Courtesy Glen B: see comments to the question.)
The symmetry of standard Normal distributions around zero and the independence of $d$ from $(a,b,c)$ imply $(a,b,c,-d)$ has the same distribution as $(a,b,c,d),$ whence the distribution of $|ab-cd|$ is the same as that of $|ab-c(-d)|=|ab+cd|.$
A brief analysis of the question in this form appears in Stuart & Ord, Kendall's Advanced Theory of Statistics Volume I (3rd Ed. 1987) in the first half of exercise 11.21:

$x_r,y_r,$ $r=1,2,\ldots,$ are independent standardized normal variates.  Show that $z= x_1y_1 + x_2y_2$ is distributed exactly as $w_1-w_2,$ where the $w_j$ are independent with f.f. [probability density] $e^{-w_j}$ [NB: only for $w_j\ge 0$] and hence or otherwise that $z$ is distributed in the Laplace form $g(z)=\frac{1}{2}\exp\left(-|z|\right),$ and that $|z|$ is again exponentially distributed.

A: There is an extensive theory on the properties of random matrices, including the distribution of their determinants.
From the answer here for example you can see that if you form the matrix $W=AA^T$, where the elements of $A$ are your $\mathcal N (0,1)$ random variables, then
$$ \det W \sim \chi^2_2 \chi^2_1 $$
namely the determinant of $W$ has the distribution of a product of two independent $\chi^2$ random variables with 2 and 1 degrees of freedom.
The distribution of the product can be derived using the equivalent distribution of a product of Gamma random variables. In the general case this is expressed in terms of a modified Bessel functions, but here it reduces to a simpler expression, and using the transformation $|\det A|=\sqrt{\det W}$ it can be worked out quite easily that indeed $|\det A|$ has an exponential distribution with mean 1.
There might be a more direct and maybe simpler way of driving this property that i'm missing at the moment - maybe someone else can find it.
