Let $\varepsilon \sim N\left ( 0,I_{k} \right )$, consider the following function of $\varepsilon$,
$y=\left ( A+B\varepsilon \varepsilon {}'B{}' \right )^{^{\frac{1}{2}}}\varepsilon $,
where $A$ is a positive semi-definite symmetric $k\times k$ matrix, $B$ is a $k\times k$ matrix. I want to find the distribution or the approximate distribution of $y$, or some moments of $y$, such as $E\left ( yy{}' \right )$, $E\left ( yy{}'\otimes yy{}' \right )$.
let me calculate $\mathbb{E}[y^\top y]=\text{tr}\,\mathbb{E}[yy^\top]$ (using Isserlis theorem):
$$\mathbb{E}[y^\top y]=\sum_{ij}\mathbb{E}\bigl[A_{ij}\epsilon_j\epsilon_i\bigr]+\sum_{ijkl}\mathbb{E}\bigl[B_{ij} \epsilon_j \epsilon_l B_{kl} \epsilon_k\epsilon_i\bigr]$$ $$\qquad=\sum_{ij}A_{ij}\delta_{ij}+\sum_{ijkl}B_{ij}B_{kl}(\delta_{jl}\delta_{ki}+\delta_{ji}\delta_{lk}+\delta_{jk}\delta_{li})$$ $$\qquad=\text{tr}\,A+(\text{tr}\,B)^2+\text{tr}\,B^2+\text{tr}\,BB^\top.$$
