What's the densitiy of the product of two independent Gaussian random variables? Suppose that $X,Y$ are two scalar independent normal random variables, $X \sim N(\mu_X,\sigma_X^2)$, $Y \sim N(\mu_Y,\sigma_Y^2)$. I'm particularly interested about the case where we don't assume $\mu_X = \mu_Y = 0$.
I'm interested in the random variable $XY$. What can be said about its PDF?
There's an existing question where an answer explains that $XY$ is the difference of two chi-squared variables.
For the zero-mean case, we know that the PDF is the normal product distribution. Is there a non-zero-mean generalization of this? 
I know that there's a 1970 SIAM paper by Springer and Thompson, but I don't have access to this. Is the part which is relevant for my question publicly available somewhere?
To add to my confusion, I found a note by Bromiley, where it is argued that the product of two normal independent random variables is a normal variable again - which I thought was not the case. The argument in the linked document goes like this: Products of gaussian PDFs are gaussian. The PDF of a product of two independent RVs is their convolution. The Fourier transform of a convolution is the product of the fourier transforms. Gaussians are mapped to gaussians under the (inverse) Fourier transform. Am I misunderstanding something? Is there something wrong with the proof?
 A: $\newcommand{\+}{^{\dagger}}
 \newcommand{\angles}[1]{\left\langle\, #1 \,\right\rangle}
 \newcommand{\braces}[1]{\left\lbrace\, #1 \,\right\rbrace}
 \newcommand{\bracks}[1]{\left\lbrack\, #1 \,\right\rbrack}
 \newcommand{\ceil}[1]{\,\left\lceil\, #1 \,\right\rceil\,}
 \newcommand{\dd}{{\rm d}}
 \newcommand{\down}{\downarrow}
 \newcommand{\ds}[1]{\displaystyle{#1}}
 \newcommand{\expo}[1]{\,{\rm e}^{#1}\,}
 \newcommand{\fermi}{\,{\rm f}}
 \newcommand{\floor}[1]{\,\left\lfloor #1 \right\rfloor\,}
 \newcommand{\half}{{1 \over 2}}
 \newcommand{\ic}{{\rm i}}
 \newcommand{\iff}{\Longleftrightarrow}
 \newcommand{\imp}{\Longrightarrow}
 \newcommand{\isdiv}{\,\left.\right\vert\,}
 \newcommand{\ket}[1]{\left\vert #1\right\rangle}
 \newcommand{\ol}[1]{\overline{#1}}
 \newcommand{\pars}[1]{\left(\, #1 \,\right)}
 \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
 \newcommand{\pp}{{\cal P}}
 \newcommand{\root}[2][]{\,\sqrt[#1]{\vphantom{\large A}\,#2\,}\,}
 \newcommand{\sech}{\,{\rm sech}}
 \newcommand{\sgn}{\,{\rm sgn}}
 \newcommand{\totald}[3][]{\frac{{\rm d}^{#1} #2}{{\rm d} #3^{#1}}}
 \newcommand{\ul}[1]{\underline{#1}}
 \newcommand{\verts}[1]{\left\vert\, #1 \,\right\vert}
 \newcommand{\wt}[1]{\widetilde{#1}}$
\begin{align}
\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!{\rm P}\pars{\xi}&=
\int_{-\infty}^{\infty}\dd x\,{1 \over \root{2\pi}\sigma}\,
\exp\pars{-\,{x^{2} \over 2\sigma^{2}}}
\\ & \
\int_{-\infty}^{\infty}\dd y\,{1 \over \root{2\pi}\sigma}\,
\exp\pars{-\,{y^{2} \over 2\sigma^{2}}}\
\overbrace{\delta\pars{\xi  - xy}}^{\ds{\delta\pars{y - \xi/x} \over \verts{x}}}
\\[3mm]&={1 \over 2\pi\sigma^{2}}\int_{-\infty}^{\infty}
\exp\pars{-\,{1 \over 2\sigma^{2}}\bracks{x^{2} + {\xi^{2} \over x^{2}}}}\,
{\dd x \over \verts{x}}
\\[3mm]&={1 \over \pi\sigma^{2}}\int_{0}^{\infty}
\exp\pars{-\,{1 \over 2\sigma^{2}}\bracks{x^{2} + {\xi^{2} \over x^{2}}}}\,
{\dd x \over x}\tag{1}
\end{align}
$\ds{\delta\pars{x}}$ is the
Dirac Delta Function.

With $\ds{x \equiv A\expo{\theta/2}\,,\quad A > 0\,,\quad\theta \in {\mathbb R}}$:
\begin{align}
{\rm P}\pars{\xi}&={1 \over \pi\sigma^{2}}
\int_{-\infty}^{\infty}
\exp\pars{-\,{1 \over 2\sigma^{2}}
\bracks{A^{2}\expo{\theta} + {\xi^{2} \over A^{2}}\expo{-\theta}}}\,
\pars{A\expo{\theta/2}\,\dd\theta/2 \over A\expo{\theta/2}}
\end{align}
We can choose $\ds{A}$ such that
$\ds{A^{2} = {\xi^{2} \over A^{2}}\quad\imp\quad A = \verts{\xi}^{1/2}}$:
\begin{align}
{\rm P}\pars{\xi}&={1 \over 2\pi\sigma^{2}}
\int_{-\infty}^{\infty}
\exp\pars{-\,{\verts{\xi} \over \sigma^{2}}\cosh\pars{\theta}}\,\dd\theta
={1 \over \pi\sigma^{2}}
\int_{0}^{\infty}
\exp\pars{-\,{\verts{\xi} \over \sigma^{2}}\cosh\pars{\theta}}\,\dd\theta
\end{align}

$$\color{#00f}{\large%
{\rm P}\pars{\xi} = {1 \over \pi\sigma^{2}}\,
{\rm K}_{0}\pars{\verts{\xi} \over \phantom{2}\sigma^{2}}}
$$
where $\ds{{\rm K}_{\nu}\pars{z}}$ is a
Second Kind Bessel Function.

A: This should probably be a comment but I cannot comment because of rep limit.
However, we know for sure that the product is NOT normally distributed itself. You can try in MATLAB
x=100*(randn(100000,1)+5);
 y=60*(randn(100000,1)+4);
 hist(x.*y,100);
and see that the resulting distribution is skewed. 
However, if the means of the two distributions are far apart (say replace the mean of y from 4 to 1000) then the product starts to look much more bell-shaped. 
A: Let me address your confusion about the argument. The density of the product of two independent random variables is not their convolution of their densities. It is more complicated. And it does not, therefore, correspond to the additive framework after you take the Fourier transform. This is possible, but it requires the Mellin transform.
A: I find this paper may be useful. 
Let $(X, Y)$ a bivariate Normal distribution with independent variables and parameters: $ \mu_{x}, \mu_{y}, \sigma_{x}, \sigma_{y} $ and let $Z=XY$. 
The author argues that
$$f_{Z}(z)=\int_{-\infty}^{\infty} \frac{1}{2 \pi|y| \sigma_{x} \sigma_{y}} \exp \left\{-\frac{1}{2 \sigma_{x}^{2}}\left(\frac{z}{y}-\mu_{x}\right)^{2}-\frac{1}{2 \sigma_{y}^{2}}\left(y-\mu_{y}\right)^{2}\right\} d y.$$
 and that there is no analytic solution. The author approximates by numerical method.
The first three moments are determined by
\begin{array}{l}{E(Z)=\mu_{x} \mu_{y}} \\ {V(Z)=\mu_{y}^{2} \sigma_{x}^{2}+\left(\mu_{x}^{2}+\sigma_{x}^{2}\right) \sigma_{y}^{2}=\left(1+\delta_{x}^{2}+\delta_{y}^{2}\right) \sigma_{x}^{2} \sigma_{y}^{2}}\\
\alpha_{3}(Z)=\frac{6 \delta_{x} \delta_{y} \sigma_{x}^{3} \sigma_{y}^{3}}{\left(\left(1+\delta_{x}^{2}+\delta_{y}^{2}\right) \sigma_{x}^{2} \sigma_{y}^{2}\right)(3 / 2)}.\end{array}
