Solution (source code)

= Solution

Taking the <transpose> gives $M^T=(X^T+X)/\sqrt2=M$, so $M$ is a <symmetric matrix>. For $i<j$, $M_{ij}=(X_{ij}+X_{ji})/\sqrt2$ has the <standard normal distribution>, whereas $M_{ii}=\sqrt2X_{ii}$ has the <normal distribution> $N(0,2)$. Distinct upper-triangular entries depend on disjoint sets of the original <independent random variables>, so they remain <independent>.

Consequently the <joint probability density function>, relative to $\prod_{i\leq j}dm_{ij}$, is
$$
\boxed{
f((m_{ij})_{i\leq j})
=(4\pi)^{-n/2}(2\pi)^{-n(n-1)/4}
\exp\!\left(-\frac14\sum_i m_{ii}^2-\frac12\sum_{i<j}m_{ij}^2\right)
=\frac{\exp[-\operatorname{tr}(m^2)/4]}
{(4\pi)^{n/2}(2\pi)^{n(n-1)/4}}.}
$$
Here $\operatorname{tr}(m^2)=\sum_i m_{ii}^2+2\sum_{i<j}m_{ij}^2$ because $m$ is a <symmetric matrix>. This is the \b[<Gaussian orthogonal ensemble>], with off-diagonal <variance> one and diagonal <variance> two; the rescaled <matrix> $M/\sqrt n$ has the corresponding $1/n$ off-diagonal <variance>.

We will also use its entry <covariance>. Expanding in the <independent> entries of $X$, for which $\mathbb E[X_{ij}X_{kl}]=\delta_{ik}\delta_{jl}$, gives
$$
\boxed{\mathbb E[M_{ij}M_{kl}]
=\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}.}
$$
The symbols $\delta_{ab}$ are <Kronecker deltas>. This formula includes the diagonal case, where the two terms both contribute.