What's the densitiy of the product of two independent Gaussian random variables?

$\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.

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.

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.