Let $S_n$ be a unit n-sphere in $\mathbb{R}_n$: $$ \sum_{i=1}^n x_i^2 = 1 $$
How to generate a probability measure $\mu$ that is uniform on $S_n$?
Well, there is an answer long time ago: You essentially take $y_i$ to be independent random Gaussian variables (with the same mean and variance); in the end you normalize them $x_i = y_i / \sqrt{\sum_{i=1}^n y_i^2}$ to get $\mu$.
One may speculate that the Wick theorem for Gaussian distribution still holds for $\mu$, i.e. correlations can be factorized into pairs, for example
$$
\mathbb{E}[x_i x_j x_k x_l] = \mathbb{E}[x_i x_j ] \mathbb{E}[x_k x_l ]+ \mathbb{E}[x_i x_k ] \mathbb{E}[x_j x_l ] + \mathbb{E}[x_i x_l ] \mathbb{E}[x_j x_k ]
$$
The rationale is that we can first generate the independent normal random variable $y_i$. These $y_i$s satisfy Wick theorem, and we normalize correlators after the factorization.
Is this true? Can we prove that Wick theorem for $\mu$ using the rotational invariance?
Update: According to joriki, the Wick theorem above does not hold. However, using the calculation by ben, we have a modified Wick contraction: $$ \mathbb{E}_{\mathbb{S}^{n-1}}[x_i x_j x_k x_l] = \frac{1}{f(2,n)} \mathbb{E}_{\mathbb{R}^{n}}[x_i x_j x_k x_l]\\ = \frac{1}{f(2,n)} ( \mathbb{E}_{\mathbb{R}^{n}}[x_i x_j ] \mathbb{E}_{\mathbb{R}^{n}}[x_k x_l ]+ \mathbb{E}_{\mathbb{R}^{n}}[x_i x_k ] \mathbb{E}_{\mathbb{R}^{n}}[x_j x_l ] + \mathbb{E}_{\mathbb{R}^{n}}[x_i x_l ] \mathbb{E}_{\mathbb{R}^{n}}[x_j x_k ]) \\ = \frac{f(1,n)^2}{f(2,n)} ( \mathbb{E}_{\mathbb{S}^{n-1}}[x_i x_j ] \mathbb{E}_{\mathbb{S}^{n-1}}[x_k x_l ]+ \mathbb{E}_{\mathbb{S}^{n-1}}[x_i x_k ] \mathbb{E}_{\mathbb{S}^{n-1}}[x_j x_l ] + \mathbb{E}_{\mathbb{S}^{n-1}}[x_i x_l ] \mathbb{E}_{\mathbb{S}^{n-1}}[x_j x_k ]) $$ For a general 2k point correlator, we need a factor of $\frac{f(1,n)^k}{f(k,n)}$. The contraction rules remain the same as the Wick theorem.
Nice idea, but unfortunately not.
If we choose $x_i=x_j=y$ and $x_k=x_l=z$, the second and third terms vanish by symmetry, and we’re left with $\mathbb E\left[y^2z^2\right]=\mathbb E\left[y^2\right]\mathbb E\left[z^2\right]$, which says that $y^2$ and $z^2$ are uncorrelated. That’s indeed the case for the Gaussian distribution, but not for the uniform distribution on the sphere, where these are negatively correlated. This correlation is introduced by the normalization.