Normal approximation to the pointwise/Hadamard/Schur product of two multivariate Gaussian/normal random variables 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.
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).
 A: $\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 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$. 
