We know the product of two independent Normal random variables has a normal product distribution, or Variance Gamma distribution if they are correlated.
But, what if there are three Normal random variables?
So, here is the question: Suppose $x,y,z$ are three independent normal random variables ($x, y, z\sim N(0,1)$), what's the distribution of $xy+xz+yz$?
On
The probability density function $f$ is given by $$f(x)=\begin{cases} \dfrac1{\surd3}\mathrm e^x & \text{if $x<0$}, \\ \dfrac2{\surd3}\mathrm e^x[1-\Phi(\sqrt{3x})] & \text{if $x\geqslant0$}, \end{cases}$$where $\Phi(x):=\dfrac1{\surd(2\pi)}\int_{-\infty}^x\exp\dfrac{-t^2}2\mathrm dt$ is the standard-normal cumulative distribution function.
It is convenient to work at first with $2(YZ+ZX+XY)$, since this is$$(X+Y+Z)^2-(X^2+Y^2+Z^2).$$This may be written as $R^2(H^2-1)$, where $R:=\sqrt{X^2+Y^2+Z^2}$, $H:=0$ when $R=0$, and $$ H:=\frac{|X+Y+Z|}{|R|}\quad(R\neq0).$$Since $X$, $Y$, an $Z$ are independent, their joint distribution function is the product of their individual distribution functions. From the definition of the standard-normal distribution, and by expressing the product of exponentials as the exponential of their sum, it is easy to see that the density in $(x,y,z)$ space then depends only on (the square of) $r:=\surd(x^2+y^2+z^2)$. So, if we divide the space into spherical shells centered on the origin, the density is constant on each shell. By choosing the equatorial plane of the shell of radius $r$ to be $x+y+z=0$, and slicing into thin, equally spaced, rings parallel to the equator, it is easy to show that the “mass” of each ring is the same. The plane of each ring has the form $x+y+z=hr$, where $0\leqslant |h|\leqslant\surd3$. Thus the random variable $H$ has a uniform distribution with support $[0\,\pmb,\, \surd3]$. It is clear also that $H$ and $R$ are independent.
It remains to find the distributions of $R^2$ and $H^2-1$. The former is $\chi^2_3$ by a standard result:$$f_{R^2}(x)=\chi^2_3(x)=\frac1{\surd(2\pi)}x^{1/2}\mathrm e^{-x/2}\quad(x\geqslant0).$$ The distribution of $H^2-1$ is obtained by first considering the cumulative distribution function of $H$, getting that of $H^2$ from it, differentiating, and then shifting by $1$: $$f_{H^2-1}(x)=\frac{\pmb1\{-1<x\leqslant2\}}{2\surd3\surd(x+1)}.$$ Next we use the formula for the distribution $f_{12}$ of the product of two independent random variables with distribution functions $f_1$ and $f_2$: $$f_{12}(x)=\int_{-\infty}^\infty f_1(t)f_2\left(\frac xt\right)\frac{\mathrm dt}{|t|}.$$Thus, substituting $f_1=f_{H^2-1}$ and $f_2=\chi_3^2$ , and discarding the range of integration ($t<0$) outside the support of $f_2=\chi_3^2$ , yields$$f_{12}(x)=\int_0^\infty\frac{t^{1/2}\mathrm e^{-t/2}}{\surd(2\pi)}\frac{\pmb1\{-1<x/t<2\}}{2\surd3\surd(x/t+1)}\frac{\mathrm dt}{t}=\int_0^\infty\frac{\pmb1\{-t<x<2t\}\mathrm e^{-t/2}}{2\surd(6\pi)\surd(x+t)}\mathrm dt.$$After the substitution $u^2=x+t$, with $\mathrm dt=2u\,\mathrm du$, this may be written$$f_{12}(x)=\frac{\mathrm e^{x/2}}{2\surd(6\pi)}\int_0^\infty\pmb1\{u^2>3x/2\}\mathrm e^{-u^2/2}\,\mathrm du.$$This reduces to $f_{12}(x)=\dfrac{\mathrm e^{x/2}}{2\surd3}$ when $x<0$, and $f_{12}(x)=\dfrac{\mathrm e^{x/2}}{\surd3}\left[1-\Phi\left(\sqrt{3x/2}\right)\right]$ when $x\geqslant0$. Finally, the distribution is scaled according to the formula $$f_{\alpha S}(x)=\frac1\alpha f_S\left(\frac x\alpha\right)$$for any random variable $S$ and scaling constant $\alpha$, which in this case is $\frac12$.
Remark: This approach may readily be generalized to any number of independent standard-normal random variables $X_i\;(i=1,...,n)$ to obtain the distribution of $\sum_{i<j}X_iX_j$: As in the above case ($n=3$), the distribution can be boiled down to that of the product of just two independent random variables—one having a $\chi_n^2$ distribution while the other's distribution is derived, as above, from a uniform distribution.
On
Alternative solution:
Remark: I use Maple to calculate the integral. I get the same result as @John Bentin's one.
Denote $w = [x, y, z]^\mathsf{T}$. We have $xy + yz + zx = \frac{1}{2} w^\mathsf{T} A w$ where $$A = \left( \begin{array}{ccc} 0 & 1 & 1 \\ 1 & 0 & 1 \\ 1 & 1 & 0 \\ \end{array} \right). $$ Let $A = U\mathrm{diag}(-1, -1, 2)U^\mathsf{T}$ be the eigendecomposition of $A$ where $U$ is an orthogonal matrix. Then, we have $xy + yz + zx = \frac{1}{2} w^\mathsf{T} U\mathrm{diag}(-1, -1, 2)U^\mathsf{T}w$. Let $v = [v_1, v_2, v_3]^\mathsf{T} = U^\mathsf{T}w$. Then, $v_1, v_2, v_3 \sim \mathcal{N}(0, 1)$ are independent. We have $xy + yz + zx = - \frac{1}{2} (v_1^2 + v_2^2) + v_3^2$. Let $u_1 = v_1^2 + v_2^2 $ and $u_2 = v_3^2$. Then, $u_1 \sim \chi^2(2)$ and $u_2 \sim \chi^2(1)$ are independent. We have $xy + yz + zx = - \frac{1}{2}u_1 + u_2$.
Denote $T = xy + yz + zx$.
For $t> 0$, we have \begin{align} F_T(t) &= \mathrm{Pr}(xy + yz + zx \le t)\\ &= \mathrm{Pr}(- \tfrac{1}{2}u_1 + u_2 \le t)\\ &= \int_0^\infty \frac{1}{2}\mathrm{e}^{-u_1/2}\left(\int_0^{t + \tfrac{1}{2}u_1} \frac{1}{\sqrt{2\pi u_2}}\mathrm{e}^{-u_2/2} du_2\right) \mathrm{d} u_1\\ &= \mathrm{erf}(\tfrac{\sqrt{2t}}{2}) - \tfrac{1}{\sqrt{3}}\mathrm{e}^t \mathrm{erf}(\tfrac{\sqrt{6t}}{2}) + \tfrac{1}{\sqrt{3}}\mathrm{e}^t \end{align} where $\mathrm{erf}(u) = \frac{2}{\sqrt{\pi}}\int_0^u \mathrm{e}^{-v^2} \mathrm{d} v$ is the error function.
For $t\le 0$, we have \begin{align} F_T(t) &= \mathrm{Pr}(xy + yz + zx \le t)\\ &= \mathrm{Pr}(- \tfrac{1}{2}u_1 + u_2 \le t)\\ &= \int_0^\infty \frac{1}{\sqrt{2\pi u_2}}\mathrm{e}^{-u_2/2} \left(\int_{2u_2 - 2t}^\infty \frac{1}{2}\mathrm{e}^{-u_1/2} \mathrm{d}u_1\right) \mathrm{d} u_2\\ &= \tfrac{1}{\sqrt{3}}\mathrm{e}^t. \end{align}
Thus, we have $$F_T(t) = \left\{\begin{array}{ll} \tfrac{1}{\sqrt{3}}\mathrm{e}^t & t\le 0 \\ \mathrm{erf}(\tfrac{\sqrt{2t}}{2}) - \tfrac{1}{\sqrt{3}}\mathrm{e}^t \mathrm{erf}(\tfrac{\sqrt{6t}}{2}) + \tfrac{1}{\sqrt{3}}\mathrm{e}^t & t > 0. \end{array} \right.$$ Then, the PDF of $T$ is given by $$f_T(t) = \left\{\begin{array}{ll} \tfrac{1}{\sqrt{3}}\mathrm{e}^t & t\le 0 \\ \tfrac{1}{\sqrt{3}}\mathrm{e}^t\left(1 - \mathrm{erf}(\tfrac{\sqrt{6t}}{2})\right) & t > 0. \end{array} \right.$$
Not a full answer, but I wanted to be able to provide an image. I did $10^8$ simulations of this distribution and plotted a density histogram. Just looking at the shape suggests the PDF is not simple; the asymmetry was surprising to me at first but it makes sense upon some reflection.
Frankly, I was surprised at how quickly Mathematica performed the simulation. Its implementation of standard normal variates must undoubtedly be extremely efficient.