What is the distribution of a random variable that is the product of the two normal random variables ?

What is the distribution of a random variable that is the product of the two normal random variables ?

Let $X\sim N(\mu_1,\sigma_1), Y\sim N(\mu_2,\sigma_2)$ and $Z=XY$

That is, what is its probability density function, its expected value, and its variance ?

I'm kind of stuck and I can't find a satisfying answer on the web. If anybody knows the answer, or a reference or link, I would be really thankful...

I will assume $X$ and $Y$ are independent. By scaling, we may assume for simplicity that $\sigma_1 = \sigma_2 = 1$. You might then note that $XY = (X+Y)^2/4 - (X-Y)^2/4$ where $X+Y$ and $X-Y$ are independent normal random variables; $(X+Y)^2/2$ and $(X-Y)^2/2$ have noncentral chi-squared distributions with $1$ degree of freedom. If $f_1$ and $f_2$ are the densities for those, the PDF for $XY$ is $$ f_{XY}(z) = 2 \int_0^\infty f_1(t) f_2(2z+t)\ dt$$

For the special case that both Gaussian random variables $X$ and $Y$ have zero mean and unit variance, and are independent, the answer is that $Z=XY$ has the probability density $p_Z(z)={\rm K}_0(|z|)/\pi$. The brute force way to do this is via the transformation theorem: \begin{align} p_Z(z)&=\frac{1}{2\pi}\int_{-\infty}^\infty{\rm d}x\int_{-\infty}^\infty{\rm d}y\;{\rm e}^{-(x^2+y^2)/2}\delta(z-xy) \\ &= \frac{1}{\pi}\int_0^\infty\frac{{\rm d}x}{x}{\rm e}^{-(x^2+z^2/x^2)/2}\\ &= \frac{1}{\pi}{\rm K}_0(|z|) \ . \end{align}