I have a question about the bivariate normal r.v.'s
Given $X, Y \sim \operatorname{Normal}(0,1)$ with correlation coefficient $\rho$. Let $Z=\max(X,Y)$. Show that $\operatorname E Z^2=1$.
My attempt:
I found that $\operatorname E Z=\frac{1}{2}[\operatorname E(X+Y+|X-Y|)]=\frac{1}{2} \operatorname E(|X-Y|)=\sqrt{\frac{1-\rho}{\pi}}$. I found this by obtaining the distribution of $U=X-Y$, which is Normal as well, $U\sim \operatorname{Normal}(0,2(1-\rho))$.
Now, $\operatorname E Z^2=\operatorname{Var} Z+(\operatorname E Z)^2 = \frac{1}{4} \operatorname{Var}(X+Y+|U|) + \frac{1-\rho}{\pi}$.
Then, $\operatorname{Var}(X+Y+|U|)=\operatorname{Var} X+\operatorname{Var} Y + \operatorname{Var}|U| + 2\operatorname{Cov}(X,Y)+2\operatorname{Cov}(X,|U|) + 2\operatorname{Cov}(Y,|U|)$.
$$\operatorname{Var}(X+Y+|U|) = 1 + 1 + \operatorname{Var}|U| + 2\rho + 2\operatorname{Cov}(X,|U|) + 2\operatorname{Cov}(Y,|U|).$$
My question is, is there any way to calculate these three covariances because I couldn't find a way to tackle them without doing messy integration? Any other approaches are highly appreciated.