The vector $(X,Y)$ is gaussian and $(U,V)=(X+Y,X-Y)$ is a linear function of $(X,Y)$ hence $(U,V)$ is gaussian as well. In particular, the distribution of $(U,V)$ is characterized by its mean vector $M$ and covariance matrix $C$. By definition,
$$
M=\begin{pmatrix}\mathbb E(U) \\ \mathbb E(V)\end{pmatrix},\qquad
C=\begin{pmatrix}\mathrm{var}(U) & \mathrm{cov}(U,V)\\ \mathrm{cov}(U,V)&\mathrm{var}(V) \end{pmatrix},
$$
with
- $\mathbb E(U)=\mu_X+\mu_Y$, $\mathbb E(V)=\mu_X-\mu_Y$,
- $\mathrm{var}(U)=\mathrm{var}(X)+\mathrm{var}(Y)+2\mathrm{cov}(X,Y)=\sigma_X^2+\sigma_Y^2+2\varrho$,
- $\mathrm{var}(V)=\mathrm{var}(X)+\mathrm{var}(Y)-2\mathrm{cov}(X,Y)=\sigma_X^2+\sigma_Y^2-2\varrho$,
- $\mathrm{cov}(U,V)=\mathrm{var}(X)-\mathrm{var}(Y)=\sigma_X^2-\sigma_Y^2$.
Except when $\varrho^2=\sigma_X^2\sigma_Y^2$, the distribution of $(U,V)$ has a density $f_{U,V}$, defined by
$$
f_{U,V}(u,v)=\frac1{2\pi\sqrt{\det C}}\exp\left(-\frac12\left(u-\mathbb E(U),v-\mathbb E(V)\right)^*C^{-1}\left(u-\mathbb E(U),v-\mathbb E(V)\right)\right).
$$
If $f(x)\,dx$ is a probability distribution with expected value $0$ and variance $1$, and the distribution of $X_i$ is
$$f\left(\dfrac{x-\mu_i}{\sigma_i}\right)\cdot\dfrac{dx}{\sigma_i},$$
and $X_i$ are independent, then certainly the distribution of $X_1+\cdots+X_n$ has expected value $\mu_1+\cdots+\mu_n$ and variance $\sigma_1^2+\cdots+\sigma_n^2$. Also, the higher cumulants would add together in the same way. (The fourth cumulant, for example, is $\mathbb E((X-\mu)^4) - 3(\mathbb E((X-\mu)^2))^2$, and the coefficient $3$ is the only number that makes this functional additive in the sense that the fourth cumulant of a sum of independent random variables is the sum of their fourth cumulants.)
We've tacitly assumed $\sigma_i<\infty$. I think if $\sum_{i=1}^\infty\sigma_i^2=\infty$, then as $n$ grows, the distribution would approach a normal distribution (I'm not recalling the appropriate generalization of the central limit theorem clearly enough to state it precisely.) But what happens for small $n$ is another question, and the answer would depend on what function $f$ is.
I said above that $f(x)\,dx$ has expectation $0$ and variance $1$. But one can also have perfectly good location-scale families in which the expectation, and a fortiori, the variance, do not exist. The most well-known case is the Cauchy distribution. The simplest result there is that $(X_1+\cdots+X_n)/n$ actually has the same Cauchy distribution as $X_1$ if these $n$ variables are i.i.d. It doesn't get narrower. So a lot depends on which function $f$ is.
Best Answer
That's the way to goal.
Let $X= \mu_1+\sigma_1~U~,~ Y= \mu_2+\sigma_2~(\rho~U+\sqrt{1-\rho^2}~V)$ where $U,V\mathop{\sim}\limits^\textrm{iid}\mathcal N(0,1)$ and $\rho=c/\sigma_1\sigma_2$
$$\begin{align}\mathsf {Cov}(X^2, Y^2) ~=~& \mathsf {Cov}\Big(\big(\mu_1+\sigma_1U\big)^2, \big(\mu_2+\sigma_2(\rho~U+\sqrt{1-\rho^2}~V)\big)^2\Big) \\[1ex] ~=~& \mathsf{Cov}\Big(\mu_1^2+2\mu_1\sigma_1U+U^2~,~ \mu_2^2+2\mu_2\sigma_2(\rho~U+\sqrt{1-\rho^2}~V)+\sigma_2^2\big(\rho^2U^2+2\rho\sqrt{1-\rho^2}~UV+(1-\rho^2)V^2\big)\Big) \end{align}$$
Simplify with the bilinearity of covariance, making use of: $\mathsf {Cov}(U,U) = 1, \mathsf {Cov}(U,V)=0, \mathsf{Cov}(U, U^2)=0,\mathsf {Cov}(U,UV)=0$ etc. (...why?)