Let $A_d$ be the area of a triangle whose vertices are chosen uniformly at random from a unit sphere in $\mathbb{R}^d$.

I claim that for $d\ge 2$ and $m\ge 1$:

$$
E(A_d^{2m})=\frac{3}{4^m} \prod _{q=1}^{m-1} \frac{3 d+6 m-2 q-6}{d+2 m-2 q-2}\prod _{q=1}^m \frac{d+2 m-2 q-1}{d+2 m-2 q}\\
= \frac{3\ \Gamma \left(\frac{d}{2}\right)^2 \Gamma \left(\frac{d-1}{2}+m\right) \Gamma \left(\frac{3 d}{2}+3 m-3\right)}{4^m\  \Gamma \left(\frac{d-1}{2}\right) \Gamma \left(\frac{d}{2}+m-1\right)
   \Gamma \left(\frac{d}{2}+m\right) \Gamma \left(\frac{3 d}{2}+2 m-2\right)}
$$

For $d=2$ (and any $m\ge 1$) this simplifies to the established formula:

$$E(A_2^{2m}) = \frac{(3m)!}{16^m\ (m!)^3}$$

For $m=1$ (and any $d\ge 2$) it simplifies to:

$$
E(A_d^2)=\frac{3(d-1)}{4d}
$$

To prove the general formula, first note that the squared area of a triangle can be described in terms of a Grammian determinant:

$$
A_d^2 = \frac{1}{4} \det{\left(s_i \cdot s_j\right)}
$$

where the triangle has vertices $v_0, v_1, v_2$ and:

$$
s_i = v_i - v_0, \: i=1,2
$$

For $d\ge3$, we can always rotate a triangle with vertices on the sphere into the configuration:

$$\begin{array}{rcl}
v_0 & = & e_0 \\
v_1 & = & \cos(\theta_1)\, e_0 + \sin(\theta_1)\, e_1 \\
v_2 & = & \cos(\theta_2)\, e_0 + \sin(\theta_2)\cos(\phi_2)\, e_1 + \sin(\theta_2)\sin(\phi_2)\, e_2
\end{array}$$

The expectation values of the even moments can then be expressed as an integral over three coordinates of a suitably weighted version of the squared area raised to a power:

$$
E(A_d^{2m}) =
\frac{(d-2)\ \Gamma \left(\frac{d}{2}\right)}{2 \pi ^{3/2}\ \Gamma \left(\frac{d-1}{2}\right)}
\int_{0}^\pi
\int_{0}^\pi
\int_{0}^\pi
(A_d^2)^m \sin(\theta_1)^{d-2} \sin(\theta_2)^{d-2} \sin(\phi_2)^{d-3}
\,d\theta_1\,d\theta_2\,d\phi_2
$$

where the vertex $v_1$ is a representative of a $(d-2)$-sphere of radius $\sin(\theta_1)$ over which it can be rotated while keeping $v_0$ fixed, and $v_2$ is a representative of a $(d-3)$-sphere of radius $\sin(\theta_2)\sin(\phi_2)$ over which it can be rotated while keeping $v_0, v_1$ fixed, and the weights incorporate the measures of these spheres with respect to the whole $(d-1)$-sphere.

For $m=1$, the integrand expands as a sum of products of non-negative integer powers of sines, and the integral can be carried out explicitly to obtain:

$$
E(A_d^2)=\frac{3(d-1)}{4d}
$$

For $m\ge 2$, we can integrate by parts to obtain the recursion relation:

$$
E(A_d^{2(m+1)}) = \frac{(d-1) (3 d+4 m)}{4 d^2} E(A_{d+2}^{2m})
$$

The general formula then follows by induction.

Although we derived this formula for even moments, it also gives correct values for odd moments using half-integer values for $m$, including the average area if we set $m=1/2$:

$$
E(A_d) = \frac{3\ \Gamma \left(\frac{3 (d-1)}{2}\right) \Gamma \left(\frac{d}{2}\right)^3}{2\ \Gamma \left(\frac{d-1}{2}\right)^2 \Gamma \left(\frac{d+1}{2}\right) \Gamma \left(\frac{3 d}{2}-1\right)}
$$

For example:

$$\begin{array}{rcl}
E(A_2) & = & \frac{3}{2\pi} \\
E(A_3) & = & \frac{\pi}{5}
\end{array}$$