Covariance of squared bivariate normal random vector $\mathbf{Y}=(Y_1, Y_2)^T$ is a bivariate Gaussian random vector with mean zero.
I want to prove that $Cov[Y_1^2, Y_2^2]=2\times(Cov[Y_1, Y_2])^2$.
I know that
$$Cov[Y_1^2, Y_2^2]=E[Y_1^2Y_2^2]-E[Y_1^2]E[Y_2^2]=E[Y_1^2Y_2^2]-Var[Y_1]Var[Y_2]$$
and
$$2(Cov[Y_1, Y_2])^2=2(E[Y_1Y_2]-E[Y_1]E[Y_2])^2=2(E[Y_1Y_2])^2$$ as $\mu_{Y_1}$,$\mu_{Y_2}$ are 0.
What is the next step I should proceed?

Some facts: 

  
*
  
*If $(Y_1,Y_2)$ is jointly normal (with mean 0) with parameters, $\rho, \sigma_1$ and $\sigma_2$, then $$Y_1 | Y_2 \sim N\left(\frac{\rho\sigma_1}{\sigma_2}Y_2, (1-\rho^2)\sigma_1^2\right)$$ 
  
*The fourth central moment of a normal random variable is $\mathbb{E}[Y] = 3\sigma^4$. 

Using these facts, we can write, 
$$\begin{eqnarray}
\mathbb{E}[Y_1^2Y_2^2] &=& \mathbb{E}[\mathbb{E}[Y_1^2Y_2^2 | Y_2]]\\
&=& \mathbb{E}[Y_2^2\mathbb{E}[Y_1^2|Y_2]]\\
&=& (1-\rho^2)\sigma_1^2\mathbb{E}[Y_2^2] + \frac{\rho^2\sigma_1^2}{\sigma_2^2}\mathbb{E}[Y_2^4]\\
&=& (1-\rho^2)\sigma_1^2\sigma_2^2 + 3\rho^2\sigma_1^2\sigma_2^2
\end{eqnarray}$$
Substitute this into your expression for $Cov(Y_1^2,Y_2^2)$ to get $2\rho^2\sigma_1^2\sigma_2^2$ which is identically $2(Cov(Y_1,Y_2))^2$.
