Solution 1:

Let $X_0, X_1, \dots, X_n$ be i.i.d. standard normal vectors in $\mathbb{R}^n$ (so each $X_i \sim \mathcal{N}(0, I_n)$). Writing $Y_i = X_i - X_0$ for $i = 1, \dots, n$, we have that the $n$-volume of the $n$-simplex with vertices $X_0, X_1, \dots, X_n$ is equal to $$\frac{1}{n!} |\det(Y_1, \dots, Y_n)|$$ where we consider $Y_1, \dots, Y_n$ as column vectors.

Define $(W_1, \dots, W_n) = (Y_1, \dots, Y_n)^T$, i.e. $W_{i, j} = X_{j, i} - X_{0, i}$, so $W_1, \dots, W_n$ are independent, and $W_i \sim \mathcal{N}(0, \Sigma)$, where the covariance matrix $\Sigma$ has $2$'s on the diagonal and $1$'s off the diagonal. Note that $J_n$ (the matrix of ones) has eigenvalues $n, 0, \dots, 0$, hence since $\Sigma = I_n + J_n$, $\Sigma$ has eigenvalues $n+1, 1, \dots, 1$ and thus $\det \Sigma = n+1$. Now, defining $Z_i = \Sigma^{-1/2} W_i$ for $i = 1, \dots, n$, we have that $Z_1, \dots, Z_n$ are independent with each $Z_i \sim \mathcal{N}(0, I_n)$, and also that $$\det(Y_1, \dots, Y_n) = \det(W_1, \dots, W_n) = \det(\Sigma^{1/2}Z_1, \dots, \Sigma^{1/2}Z_n) = \det \Sigma^{1/2} \cdot \det(Z_1, \dots, Z_n).$$ It follows that the desired expected volume is $$\frac{\sqrt{n+1}}{n!} \mathbb{E}[|\det(Z_1, \dots, Z_n)|]$$ for independent $Z_1, \dots, Z_n \sim \mathcal{N}(0, I_n)$. To finish, we compute $\mathbb{E}[|\det(Z_1, \dots, Z_n)|]$.

Let $Z_1', \dots, Z_n'$ be the result of performing the Gram-Schmidt process to $Z_1, \dots, Z_n$ without normalizing, so for each $k$, we have $\mathrm{span}(Z_1', \dots, Z_k') = \mathrm{span}(Z_1, \dots, Z_k)$, and we inductively define $Z_k' = Z_k - P_kZ_k$ (with $Z_1' = Z_1$), where $P_k$ is the orthogonal projection onto $\mathrm{span}(Z_1', \dots, Z_{k-1}')$. Notably, these are all elementary column operations, so $\det(Z_1', \dots, Z_n') = \det(Z_1, \dots, Z_n)$, and $Z_1', \dots, Z_n'$ are orthogonal, so $|\det(Z_1', \dots, Z_n')| = \prod_{k=1}^n |Z_k'|$. Equivalently, we have $Z_k' = P_k' Z_k$, where $P_k'$ is the orthogonal projection onto the orthogonal complement of $\mathrm{span}(Z_1', \dots, Z_{k-1}')$, so $Z_k'$ can be seen as a standard normal vector on this $(n-k+1)$-dimensional space. This means that conditioning on $Z_1', \dots, Z_{k-1}'$, $|Z_k'|$ has the chi distribution with $n-k+1$ degrees of freedom, so in fact $|Z_k'|$ is independent of $Z_1', \dots, Z_{k-1}'$ with $$\mathbb{E}[|Z_k'|] = \sqrt{2} \frac{\Gamma((n-k+2)/2)}{\Gamma((n-k+1)/2)}.$$ It follows that all $|Z_k'|$ are independent, giving \begin{align*} \mathbb{E}[|\det(Z_1, \dots, Z_n)|] &= \prod_{k=1}^n \mathbb{E}[|Z_k'|]\\ &= \prod_{k=1}^n \sqrt{2} \frac{\Gamma((n-k+2)/2)}{\Gamma((n-k+1)/2)} \\ &= \prod_{k=1}^n \sqrt{2} \frac{\Gamma((k+1)/2)}{\Gamma(k/2)} \\ &= 2^{n/2} \frac{\Gamma((n+1)/2)}{\Gamma(1/2)} \end{align*} so the expected volume is $2^{n/2} \frac{\Gamma((n+1)/2) \sqrt{n+1} }{\Gamma(1/2) n!}$. At $n = 3$ (the given case), this is $\frac{2}{3} \sqrt{\frac{2}{\pi}}$.

Higher moments can be computed in the same way, using the corresponding higher moments of the chi distribution.