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.
Best Answer
As suggested by Jim B, The $n$th moment of a circular distribution is defined as, $$ m_n = \int_\phi^{\phi+2\pi} \rho(\theta) e^{in\theta}\, \mathrm{d\theta}, $$ where $\phi$ is any angle and the probability density $\rho$ is periodic with period $2\pi$. It is illuminating to see how these moments vary with coordinate changes. Rotating $\rho$ by an angle $\varphi$ we have, $$ m_n' = \int_\phi^{\phi+2\pi} \rho(\theta-\varphi) e^{in\theta}\, \mathrm{d\theta} = e^{in\varphi}\int_{\phi-\varphi}^{\phi-\varphi+2\pi} \rho(\theta') e^{in\theta'}\, \mathrm{d\theta'} = e^{in\varphi}m_n $$ Where in the last step we used the fact that periodic functions have the same integral over any range that is exactly as long as their period. It is clear that $m_0$ is real and equal to 1 for well-behaved distributions. $m_1$ has a complex argument equal to the average angle, and by $e^{i \varphi}$ transforms as we would expect with rotating coordinate systems. For $n>1$, we can conclude that the complex argument of $m_{n>1}$ cannot be meaningful as a property of $\rho$, as it is not invariant with respect to rotations.
To see how we might recover the "usual" moments from the circular moments, suppose that $\rho(\theta)$ is concentrated in a small region that we may without loss put near $\theta=0$. Expanding $e^{in\theta}$ around that point, we arrive at, $$ m_n \approx \sum^\infty_{j=0}\frac{(in)^j}{j!}\int_{-\pi}^{\pi} \theta^j \rho(\theta) \, \mathrm{d\theta} $$ It is clear that the "usual" moments are all hidden together within $m_n$, summed with weights that depend on $n$. They may be isolated with polynomials of $m_n$ for many $n$, although this is not usually done.