Suppose that $(X,Y)$ is distributed ranodm 2-vectors having the normal distribution with $\mathbb{E}X=\mathbb{E}Y=0,\mathbb{Var}(X)=\mathbb{Var}(Y)=1,$ and $\text{Cov}(X,Y)=\theta\in(-1,1).$ Show the value of $\mathbb{E}(X^2Y^2).$
The joint pdf of $(X,Y)$ is $$p(x,y)=\frac{1}{2\pi\sqrt{1-\theta^2}}\exp\left\{-\frac{1}{2(1-\theta^2)} \left [ x^2-2\theta xy+y^2\right ]\right\}.$$
Then $$\mathbb{E}(X^2Y^2)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}x^2y^2\cdot p(x,y)\mathrm{d}x\mathrm{d}y.$$ Let $$\begin{cases} u=\frac{x-y}{\sqrt{2}} \\ v=\frac{x+y}{\sqrt{2}} \end{cases}\Rightarrow \mathbb{E}(X^2Y^2)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\frac{(u^2-v^2)^2}{2}\cdot p\left (\frac{u+v}{\sqrt{2}},\frac{u-v}{\sqrt{2}}\right )\mathrm{d}u\mathrm{d}v.$$ But I don't know how to deal with this integral
$$\begin{align*} &\Delta:=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\frac{(u^2-v^2)^2}{2}\cdot \exp\left \{-\frac{1}{2(1-\theta^2)}\left [ (1-\theta)u^2+(1+\theta)v^2\right ]\right \}\mathrm{d}u\mathrm{d}v\\ &\quad=4\int_{0}^{\infty}\int_{0}^{\infty}\frac{(u^2-v^2)^2}{2}\cdot \exp\left \{-\frac{1}{2(1-\theta^2)}\left [ (1-\theta)u^2+(1+\theta)v^2\right ]\right \}\mathrm{d}u\mathrm{d}v,\quad \theta\in(-1,1). \end{align*}$$