Analogue of Wick formula for orthogonal polynomials n-point correlations of Gaussian random variables can be simplified with Wick expansion.
$$ \langle x_{i_1} x_{i_2} \dots x_{i_{2n-1}} x_{i_{2n}} \rangle  = \int_{\mathbb{R}^n} x_{i_1} \dots x_{i_{2n}} e^{- \frac{1}{2}\sum x_i^2}  = \sum_{\sigma \in \text{matchings}} \; \prod_{\{i,j\}\in \sigma} \langle x_i x_j \rangle$$
A matching is a partition of $\{ 1,2,\dots ,2n\}$ into two-element sets, e.g. $\{ \{1,2\},\dots, \{2n-1,2n\}\}$.
The Chebyshev polynomials are orthogonal with respect to $\mu(dx) = \sqrt{4-x^2} dx$  but there probably isn't a way to simplify n-point correlations,
\[ \int_{S^n}  x_{i_1}\dots x_{i_{2n}} \sqrt{1 - (x_1^2+\dots + x_{n}^2)}  \; dx \]
Is there an analogue of Wick's formula for other orthogonal polynomials?  It seems rather unlikely since the measure can be arbitrary.
Maybe look at Wick expansion as the sum over matchings.  Is there a version of the Wick expansion for non-crossing matchings or for Schröder paths?

 A: The integral you want to compute is related to multivariable beta functions, the Dirichlet distribution and Dirichlet integrals. In particular one has
$$\int\int\cdots\int f(x_1+\cdots+x_n) x_1^{k_1}\cdots x_n^{k_n} dx_1\cdots dx_n$$
$$=\frac{\Gamma(k_1+1)\cdots \Gamma(k_n+1)}{\Gamma(k_1+\cdots+k_n+n)}\int_0^1 f(t)t^{k_1+\cdots+k_n+n-1}dt$$
where the first integral is over the simplex $\sum x_i\le 1$, and in your case $f(t)=\sqrt{1-t}$. Since you are considering moments the $k_i$'s will end up being half-integers and therefore the result will be some sort of generalized Catalan number.
There is also the "Free probability" point of view, which I think is closer to what you want. In fact the joint distribution of free Gaussian random variables has an expansion similar to the Wick expansion but over noncrossing matchings. Of course in free probability the role of the normal distribution is played by the Wigner semicircle distribution.
From the perspective of orthogonal polynomials I believe the correspondence goes like this: Hermite polynomials correspond to the normal distribution and complete matchings. Chebyshev polynomials of second kind correspond to the semicircle distribution and to noncrossing matchings. Charlier polynomials correspond to set partitions, Laguerre polynomials correspond to permutations etc.
