Let $X \sim \mathcal{N}\left( {{\mu \_x},\sigma \_x^2} \right)$ and $Y \sim \mathcal{N}\left( {{\mu \_y},\sigma \_y^2} \right)$ be two univariate and independent Gaussian/normal random variables and let $Z = XY$ be their product. 


It is well-known that the probability distribution of $Z$ can be approximated by the normal distribution $\mathcal{N}\left( {{\mu \_x}{\mu \_y},\mu \_x^2\sigma \_y^2 + \mu \_y^2\sigma \_x^2} \right)$ when $\frac{{{\mu \_x}}}{{{\sigma \_x}}}$ and $\frac{{{\mu \_y}}}{{{\sigma \_y}}}$ are large. See for instance [An approach to distribution of the product of two random variables](http://pldml.icm.edu.pl/pldml/element/bwmeta1.element.bwnjournal-article-doi-10_7151_dmps_1146/c/dmps.1146.pdf).


I’m looking for similar results in the multivariate case : $X \sim {\mathcal{N}\_d}\left( {{\mu \_x},{\Sigma \_x}} \right)$ , $Y \sim {\mathcal{N}\_d}\left( {{\mu \_y},{\Sigma \_y}} \right)$ and $Z$ is the pointwise/Hadamard/Schur product of $X$ and $Y$ , i.e. $Z = \left( {{x\_1}{y\_1},{x\_2}{y\_2},...,{x\_d}{y\_d}} \right)$.


More generally, I’m looking for references dealing with the properties of the pointwise/Hadamard/Schur products of multivariate random variables (as a generalization of the scalar product of univariate r.v.s).



Normal approximation to the pointwise/Hadamard/Schur product of two multivariate Gaussian/normal random variables
$\newcommand{\si}{\sigma}
\newcommand{\Si}{\Sigma}
\renewcommand{\c}{\circ}
\newcommand{\tr}{\operatorname{tr}}$


Let $X\_1:=X$, $X\_2:=Y$, $\mu\_1:=\mu\_x$, $\mu\_2:=\mu\_y$, $\Si\_1:=\Si\_x$, $\Si\_2:=\Si\_y$, $N:=\mathcal{N}\_d$. Let $Y\_i:=X\_i-\mu\_i$. Then $X\_i=\mu\_i+Y\_i$, $Y\_i\sim N(0,\Si\_i)$, $Y\_1$ and $Y\_2$ are independent, and 
\begin{equation\*}
 X\_1\c X\_2=\mu\_1\c\mu\_2+V+Y\_1\c Y\_2,\quad\text{where}\quad V:=\mu\_1\c Y\_2+\mu\_2\c Y\_1
\end{equation\*}
and $\c$ denotes the Hadamard product. Next, 
\begin{equation\*}
 V\sim N(0,\Si),\quad\text{where}\quad \Si:=D(\mu\_1)\Si\_2 D(\mu\_1)+D(\mu\_2)\Si\_1 D(\mu\_2)
\end{equation\*}
and $D(\mu)$ is the diagonal matrix whose diagonal entries are the entries of the column matrix $\mu$; we are assuming that the matrix $\Si$ is nonsingular. Hence, 
\begin{equation\*}
 \Si^{-1/2}(X\_1\c X\_2-\mu\_1\c\mu\_2)=Z+\delta,\quad\text{where}\quad \delta:=\Si^{-1/2}(Y\_1\c Y\_2)
\end{equation\*}
and $Z\sim N(0,I)$, a standard normal random vector in $\mathbb R^d$. Further, 
\begin{equation\*}
 E\|\delta\|^2\le\|\Si^{-1/2}\|^2 E\|Y\_1\c Y\_2\|^2=\|\Si^{-1}\| \tr(\Si\_1\c\Si\_2), 
\end{equation\*}
$\|v\|$ is the Euclidean norm of a vector $v$, and $\|M\|$ and $\tr M$ are the corresponding operator norm and the trace of a matrix $M$, respectively. 


Thus, we come to the following conclusion: if $\mu\_1$, $\mu\_2$, $\Si\_1$, $\Si\_2$ vary in any such way that the matrix $\Si$ is nonsingular and 
$$\|\Si^{-1}\| \tr(\Si\_1\c\Si\_2)\to0,$$ 
then 
$$X\_1\c X\_2\approx N(\mu\_1\c\mu\_2,\Si)$$ 
in the sense that $\Si^{-1/2}(X\_1\c X\_2-\mu\_1\c\mu\_2)$ converges to a standard normal random vector $Z$ in distribution. 


In the particular case when $\Si\_i=\si\_i^2 I$ for some scalars $\si\_i>0$, the condition $\|\Si^{-1}\| \tr(\Si\_1\c\Si\_2)\to0$ becomes 
\begin{equation\*}
 \min\Big(\frac{\mu\_1^2}{\si\_1^2}+\frac{\mu\_2^2}{\si\_2^2}\Big)\to\infty, 
\end{equation\*}
where $\min v$ denotes the minimum of the entries of a (say, column) matrix $v$. 


Further specializing this to the case $d=1$, we have 
\begin{equation\*}
X\_1 X\_2\approx N(\mu\_1\mu\_2,\mu\_1^2\si\_2^2+\mu\_2^2\si\_1^2)\quad\text{if}\quad \frac{\mu\_1^2}{\si\_1^2}+\frac{\mu\_2^2}{\si\_2^2}\to\infty. 
\end{equation\*}


Here, as in the paper [An approach to distribution of the product of two random variables](http://pldml.icm.edu.pl/pldml/element/bwmeta1.element.bwnjournal-article-doi-10_7151_dmps_1146/c/dmps.1146.pdf) you cited, the asymptotic variance is $\mu\_1^2\si\_2^2+\mu\_2^2\si\_2^2$, not $\mu\_1^2\si\_2^2+\mu\_2^2\si\_1^2+\si\_1^2\si\_2^2$ as you wrote, even though the two values are relatively close to each other under the condition $\frac{\mu\_1^2}{\si\_1^2}+\frac{\mu\_2^2}{\si\_2^2}\to\infty$. 


However, there are two differences between our conclusion for the special case $d=1$ and the corresponding result in the cited paper: (i) there, the statement is not rigorous, saying that something tends to a "limit" that is in fact varying with $\mu\_1$, $\mu\_2$, $\si\_1$, $\si\_2$ and (ii) the condition there for such a "convergence" appears to be that both $\frac{\mu\_1}{\si\_1}$ and $\frac{\mu\_2}{\si\_2}$ "increase", whereas here we show that it is enough to have $\frac{\mu\_1^2}{\si\_1^2}+\frac{\mu\_2^2}{\si\_2^2}\to\infty$ -- that is, it is enough that at least one of the ratios $\frac{|\mu\_1|}{\si\_1}$ or $\frac{|\mu\_2|}{\si\_2}$ go to $\infty$. 


