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