I know that the product of two independent zero mean Normal random variables follows the normal product distribution.
I want to talk about the random variable $$z = xy$$ where $x \sim \mathcal{N}(0,1)$ and $y\sim \mathcal{N}(\mu,1)$ and they are independent.
In case of $\mu \gg 1$, simulations show that the distribution of $z$ can be approximated by $\mathcal{N}(0,\mu^2)$. I want to prove this fact.
My tries : $z = xy = x(\bar{y} + \mu)=x\bar{y}+x\mu$ where $\bar{y} \sim \mathcal{N}(0,1)$.
It is clear that $x\mu \sim \mathcal{N}(0,\mu^2)$ while $\mathbb{E}[x\bar{y}]=0$ and $\textrm{var}(x\bar{y})=\textrm{var}(x)\textrm{var}(\bar{y})=1$. Thus $z \approx x\bar{y} \sim \mathcal{N}(0,\mu^2)$.
Could anyone provide more rigorous proof ?
Edit to tighten up the idea sketched below:
Consider the scaled random variable $W_{\mu} := XY/\mu$ for $\mu \neq 0$. Then we have that the characteristic function is given by \begin{align} \varphi_{W_\mu}(t) &= \mathbb{E}\left[e^{itW_{\mu}} \right] \\ &= \mathbb{E}\left[e^{i(t/\mu)XY }\right] \\ &= \frac{\exp\left[-\frac{t^2}{2(1+t^2/\mu^2)}\right]}{\sqrt{1+t^2/\mu^2}}. \end{align} Therefore as $\mu \rightarrow \infty$ you have pointwise convergence $\varphi_{W_{\mu}}(t) \rightarrow e^{-\frac{t^2}{2}}$ of the characteristic functions so that $$ \frac{XY}{\mu} \stackrel{\mathscr{D}}{\rightarrow} \mathcal{N}(0, 1). $$
There are details for you to tighten up and prove exactly, but the following should get you well on the way.
If $X \sim \mathcal{N}(0, 1)$ and $Y \sim \mathcal{N}(\mu,1)$ then the random variable $Z=XY$ has moment generating function $$ \begin{align} M_Z(t) &= \mathbb{E}\left[ e^{tXY}\right]\\ &= \frac{\exp\left[ \frac{\mu^2 t^2}{2(1-t^2)}\right]}{\sqrt{1-t^2}}\\ &= 1 + \frac{1}{2}(\mu^2 + 1)t^2 +\frac{1}{8}(\mu^4 + 6\mu^2+3)t^4 +\mathcal{O}(t^6), \end{align} $$ now if $\mu$ is large all of these polynomials, in $\mu$, are well approximated by only their leading degree monomials leading to $$ M_Z(t)\approx 1 + \frac{\mu^2 t^2}{2}+\frac{1}{2}\left(\frac{\mu^2 t^2}{2}\right)^2 + \cdots $$ or $$ M_Z(t) \approx \exp\left(\frac{\mu^2 t^2}{2} \right) $$ which is the moment generating function of a normal random variable with density $\mathcal{N}(0, \mu^2)$.